{'fraction_non_alphanumeric': 0.056730038022813685, 'fraction_numerical': 0.021191381495564005, 'mean_word_length': 3.7617380808690406, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the distribution of the time-integral of the geometric Brownian motion

We study the numerical evaluation of several functions appearing in the small time expansion of the distribution of the time-integral of the geometric Brownian motion as well as its joint distribution with the terminal value of the underlying Brownian motion. A precise evaluation of these distributions is relevant for the simulation of stochastic volatility models with log-normally distributed volatility, and Asian option pricing in the Black-Scholes model. We derive series expansions for these distributions, which can be used for numerical evaluations. Using tools from complex analysis, we determine the convergence radius and large order asymptotics of the coefficients in these expansions. We construct an efficient numerical approximation of the joint distribution of the time-integral of the gBM and its terminal value, and illustrate its application to Asian option pricing in the Black-Scholes model.


Introduction
The distributional properties of the time integral of the geometric Brownian motion (gBM) appear in many problems of applied probability, actuarial science and mathematical finance.
The probability distribution function is not known in closed form, although several integral representations have been given in the literature, as we describe next.
A closed form result for the joint distribution of the time integral of a gBM and its terminal value was obtained by Yor [34]. Denoting A (µ) t = t 0 e 2(Ws+µs) ds, where W t is a standard Brownian motion, this result reads where θ r (t) is the Hartman-Watson integral defined by ∞ 0 e − ξ 2 2t e −r cosh ξ sinh ξ sin πξ t dξ .
and r is a positive parameter. This is related to the Hartman-Watson distribution [14] which is defined by its probability density function f r (t) = θr(t) θr(s)ds . Integrating (1) over x gives an expression for the distribution of the time-average of the gBM as a double integral. For certain values of µ this can be reduced to a single integral [4]. Chapter 4 in [21] gives a detailed overview of the cases for which such simplifications are available.
The numerical evaluation of the integral (2) is challenging for small values of t [1], which motivated the search for alternative numerical and analytical approximations. We summarize below limit theorems in the literature for the joint distribution of the time integral of the gBM and its terminal value.
An asymptotic expansion for the Hartman-Watson integral θ ρ/t (t) as t → 0 at fixed ρ was derived in [28]. The leading order of this expansion is and depends on two real functions F (ρ), G(ρ). The exact definition of F and G is technical and so we postpone it to Section 4.2, see formulas (53) and (54).
Upon substitution into (1), the expansion (3) gives the leading t → 0 asymptotics of the joint distribution of the time-integral of a gBM and its terminal value The function I(a, v) appears also in the asymptotic expansion of the log-normal SABR stochastic volatility model discretized in time, as the number of time steps n → ∞ under appropriate rescaling of the model parameters, see Sec. 3.2 in [29]. The application of Laplace asymptotic methods gives an expansion for the distribution function of the time integral of the gBM, see Proposition 4 in [28] (6) where g(a, µ) is given in Proposition 4 of [28] and J(a) = inf v≥0 I(a, v).
1.1. Literature review on the small time asymptotics of the distribution of (A (µ) t , W t ). In order to put these results into context, we summarize briefly the known short maturity asymptotics of the continuous time average of the geometric Brownian motion a T := 1 T T 0 e 2(Ws+µs) ds. Dufresne [5] proved convergence in law of this quantity to a normal distribution, see Theorem 2.1 in [5] This can be formulated equivalently as a log-normal limiting distribution in the same T → 0 limit, see Theorem 2.2 in [5]. The log-normal approximation for the time-integral of the gBM is widely used in financial practice, see [5] for a literature review. When applied to the pricing of Asian options, this is known as the Lévy approximation [17].
The log-normal approximation is also used for the conditional distribution of A (µ) t at given terminal value of the gBM, as discussed for example in [22,23]. This conditional distribution appears in conditional Monte Carlo simulations of stochastic volatility models [33], see [2] for an overview and recent results.
The result (7) can be immediately extended to the joint distribution of (a T , z T ) with z T = e W T +µT . We state it without proof which can be given along the lines of the proof of Theorem 2.1 in [5].
Proposition 1.1. Denotex T = (a T , z T ) with z T := e W T +µT . As T → 0 we have: where1 = (1, 1) and The covariance matrixΣ is computed as (with W sT = √ T B s and 0 ≤ s ≤ 1) In other words, as T → 0, the joint distribution of (a T , z T ) approaches a bivariate normal distribution with correlation ρ az = √ 3 2 . In analogy with Theorem 2.2 in [5], this result can be formulated equivalently as a bivariate log-normal limiting distribution.
Another asymptotic result in the T → 0 limit is the large deviations property of the timeaverage a T . Using large deviations theory methods it was showed in [26] that, for a wide class of diffusions which includes the geometric Brownian motion, the time average of the diffusion satisfies a large deviations property in the small time limit, see Theorem 2 in [26].
For the particular case of a gBM this result reads, in our notations with J BS (a) a rate function given in closed form in Proposition 12 of [26]. For convenience the result is reproduced below in Eq. (20). In contrast to the fluctuations result of (7) which holds in a region of width a − 1 ∼ √ T around the central value a = 1, the large deviations result (13) holds for a − 1 ∼ O(1).

1.2.
Complete leading small time asymptotics. The explicit asymptotic result (4) for the small time asymptotics of the joint distribution of (A (µ) t , W t ) satisfies all these limiting results and improves on them by including all contributions up to terms of O(t). It is instructive to see explicitly the connection to the asymptotic results presented above.
The exponential factor in (6) reproduces the large deviations result (13) by identifying J(a) = 1 4 J BS (a). The emergence of the log-normal limit implied by Proposition 1.1 can be also seen explicitly by expanding the exponent in (4) around a = v = 1 using the results of Propositions 2 and 3 in [28]. Keeping the leading order terms in the expansion of F (ρ) = π 2 2 −1 −log ρ+ log 2 ρ+ O(log 3 ρ) in powers of log ρ gives the quadratic approximation The joint distribution (4) is thus approximated in a neighborhood of (v, a) = (1, 1) as where the (· · · ) in the exponent are the error terms in (14) and we defined Up to the factor C LN (a, v) and the error terms in the exponent, the expression (15) coincides with a log-normal bivariate distribution in variables (a, v). Recall the bivariate normal distribution of two random variables N(X, Y ; σ X , σ Y , ρ) with volatility and correlation parameters σ X , σ Y , ρ The density (15) corresponds to parameters (18) These parameters correspond precisely to the covariance matrixΣ appearing in the fluctuations result of Proposition 1.1.
Improving on the log-normal approximation (15) requires adding higher order terms in the coefficient C LN (a, v) and in the exponent. All these corrections are expressed in terms of F (ρ), G(ρ) which can be evaluated by expansion around ρ = 1. In this paper we present series expansions of F (ρ), G(ρ) in powers of log ρ that can be used for such evaluation, and we study in detail their properties. We determine their convergence radii and coefficient asymptotics in Proposition 4.2. We note that alternative series expansions which improve on the log-normal approximation have been considered in the literature (in the context of the expansion of the density of A (µ) t ), based on Gram-Charlier series [6]. Although closed form results are available for the functions F (ρ), G(ρ), they are inefficient for practical use because of the need to solve one non-linear equation for each evaluation. The series expansions for F (ρ), G(ρ) allow much faster and reliable numerical evaluation, inside their convergence domain. In Section 5 we construct efficient numerical approximations for F (ρ), G(ρ), based on the series expansions derived here within the convergence domain, which are easy to use in numerical applications. This gives a numerical approximation for the joint distribution of the time integral of the gBM and its terminal value, which appears in many problems of mathematical finance, such as the simulation of the SABR model [3] and option pricing in the Hull-White model [15]. As an application we discuss the numerical pricing of Asian options in the Black-Scholes model, and demonstrate good agreement with standard benchmark cases used in the literature.

Preliminaries
The mathematics of the series expansions for the functions J(a), F (ρ), G(ρ) involves the study of the analyticity domain and of the singularities of a function defined through the inverse of a complex entire function. Such functions appear in several problems of applied probability, for example in applications of large deviations theory to the short maturity asymptotics of option prices in stochastic volatility models [9] and Asian options in the local volatility model [26,27]. A similar problem is encountered in the evaluation of functions defined by the Legendre transform of a function (19) F where θ * is the solution of the equation f ′ (θ) = x. If the optimizer θ * can be found exactly at one particular point x 0 , then expansion of f (θ) around θ * gives a series expansion for F (x) in x − x 0 which can be used for numerical evaluation of the Legendre transform.
For definiteness, consider the evaluation of the rate function J BS appearing in the short maturity asymptotics of the distribution of the time-average of the gBM. This determines also the leading short maturity asymptotics of Asian options in the Black-Scholes model 1 . The rate function can be computed explicitly [26] and is given by J BS (x) : (0, ∞) → R + , defined by 1 More precisely, in the short maturity limit T → 0 the prices of out-of-the-money Asian options approach zero at an exponential rate C(K, T ) ∼ e − 1 T JBS(K/S0) , determined by the rate function J BS (x).
where ξ is the solution of the equation and ζ is the solution in [0, π] of the equation See Proposition 12 in [26]. The function J BS (x) vanishes at x = 1 and can be expanded in a series of powers of x − 1 and log x, see [26] 2 . In practical applications one is interested in the range and rate of convergence of these series expansions. We answer these questions in Proposition 4.1.
Section 3 studies the analyticity properties of the complex inverse of the function determining ξ, ζ in (20). These properties are used in Section 4 to study the singularity structure of the rate function J BS (x) and of the functions F (ρ), G(ρ) appearing in Eq. (3). The nature of the singularity on the circle of convergence determines the large-order asymptotics of the coefficients of the series expansions of a complex function around a regular point. This is made precise by the transfer results (see Flajolet and Odlyzko [7,8]) which relate the expansions around the singular point to the large-order asymptotics of the expansion coefficients. We apply these results in Section 4 to derive the leading large-order asymptotics of the expansion coefficients for the functions J BS (x) and F (ρ), G(ρ) around x = ρ = 1, see Proposition 4.2.
Besides the rigorous proofs, numerical tests also confirm a general decreasing trend of the error terms of the leading order asymptotics for the expansion coefficients.
3. On the complex inverse of x(ξ) = sinh ξ ξ An important role in the study of the expansions considered in this paper is played by the singularities of a function defined through the inverse of another function, as in the rate function J BS (x), see Eq. (20). The main issues appear already in the study of the functions ξ(x), ζ(x). For this reason, we discuss the analytic structure of these functions in details in this section.
Clearly, f is an entire function with Taylor series around zero given by We want to make sense of the inverse of f , restricted to some complex domains. First observe that the equation sinh ξ = ξ has a unique real solution ξ = 0 and infinitely many complex solutions. This shows that f is not invertible on the entire complex plane. In this section we discuss how to make sense of an inverse. Although the results will not surprise experts in complex analysis, this question has not been discussed in the present context before.

Local analysis.
Let us now restrict f to a small neighborhood of the origin. Then the range of f is in a small neighborhood of 1. Clearly, f is not locally invertible as f is an even function, that is f (ξ) = f (−ξ), and thus is not injective. A natural way to make the function at least locally invertible is to define In light of (24), we have Here, and in the sequel, unless stated otherwise, the notation √ z will only be used in the argument of even functions so as both values of √ z give the same result. Recall the following simple fact: if g ′ (z) = 0, then there is a small neighborhood D of z so that g : D → g(D) is biholomorphic (see for example Theorem 3.4.1 in [30]). Since g ′ (0) = 0, g is locally invertible at zero. We also note that g is an entire function as its Taylor series converges everywhere. We will show next that g : C → C is surjective. One says that an entire function f has finite order ρ if constants A, B can be found such that |f (z)| ≤ Ae B|z| ρ for all z ∈ C. See for example Chapter 5 in [31] and Section 1.3 in [16] for a detailed discussion of entire functions of finite order. The order of the entire function g(z) = ∞ n=0 a n z n can be expressed in terms of the coefficients of its Taylor expansion as, see Theorem 2 in Sec. 1.3 of [16] ρ = lim sup n log n − log |a n | .
Using the expansion (25) it follows that the function g has order ρ = 1 2 . Any entire function with order < 1 is surjective, see Section 5.1 of [16]. Since in our case ρ = 1/2, g : C → C is surjective.
We study next the critical points of f and g which are given by the zeros of f ′ and g ′ , respectively. To identify the zeros of g ′ , let us write 0 = η 0 < η 1 < η 2 < ... for the nonnegative solutions of the equation tan η = η (and note that all real solutions are given by Table 1 for the numerical values for k = 1, ..., 5). Now we have The zero sets of f ′ and g ′ are given by Proof. If ξ = 0, then the equation f ′ (ξ) = 0 reduces to ξ = tanh(ξ), whose solutions are exactly the numbers ξ k , k = 0. In case ξ = 0, we see by (24) that f ′ (0) = 0. The first statement follows. For the second statement, observe that g ′ (z) = f ′ ( √ z)/(2 √ z) and g ′ (0) = 1/6 = 0.
Proof. Since g maps R to R, in this proof we will denote by √ z for the square root function . Now we are ready to state the analyticity of the local inverse of g on the region that is most important for our applications.
Proof. Fix an arbitrary ε > 0. It suffices to show that h can be defined analytically on (1). First observe that the real interval I = [z 1 +ε − , ε + ] is being mapped, under the function g, to the real interval g(I) = [ω 1 + ε, 1 + ε] for suitable ε ± . Furthermore, the derivative of g is non-zero on I. Thus restricting the domain of g to a small neighborhood of I, the restricted function provides a bijection between a neighborhood of I and a neighborhood of g(I). More precisely, since the derivative of g on I is non-zero, for every ω ∈ g(I) there is some δ ω so that g can be analytically inverted on B δω (ω). By compactness of g(I), we can choose a finite cover of g(I) by the sets B δω (ω) showing that the inverse function h is well defined on a small neighborhood B δ ′ (g(I)) of g(I).
Next we extend the function h analytically to the ball B. This can be done by a standard procedure. Recall that a point a ∈ C is an asymptotic value of g if there is a curve Γ tending to ∞ so that g(z) → a az z → ∞ along Γ. Furthermore, if g ′ (c) = 0, then g(c) is a critical value. The singularity set is defined as the union of asymptotic values and critical values. It is known that for any open set D that is disjoint to the singularity set, the map g : g −1 (D) → D is a covering (see [25]).
In our case, we have already identified the set of critical values {ω k , k ≥ 1}. An elementary computation shows that the only asymptotic value is ω = 0. Now let D = B \ B δ ′ /2 (g(I)) (whereB is the closure of B). Then D is open and so g : g −1 (D) → D is a covering. This implies that along any simple polygonal arc γ starting from B δ ′ /2 (I) and staying inside B, h has an analytic extension. Indeed, by the covering map property, every ω ∈ γ has a neighborhood U ω so that g −1 (U ω ) is a union of disjoint open sets, each of which is mapped homeomorphically onto U by g. Let us now choose a finite cover of γ by neighborhoods U ω 1 , ..., U ωn so that U ω 1 ∩ B δ ′ /2 (I) = ∅ and U ω k−1 ∩ U ω k = ∅ and by induction on k we can extend h analytically to ∪U ω k . Observing that D is simply connected, the result now follows from the monodromy theorem.
3.2. Global analysis. Theorem 3.1 tells us that the inverse h(ω) can be defined with h(1) = 0. Furthermore, h is analytic and the Taylor expansion has a radius of convergence of at least 1 − ω 1 . We note that the inverse g −1 can be defined in the sense of a sheaf (also known as global analytic function, or Riemann surface). In this general sense, g −1 has a second order algebraic branch point at ω 1 = g(z 1 ) by Lemma 3.2.
Thus the function h cannot be analytically extended to ω 1 and so the power series (27) has radius of convergence exactly 1 − ω 1 . We will be mostly interested in this expansion (see Section 3.4). However, for completeness we first discuss how to visualize the sheaf g −1 .
First, observe that h can be analytically extended to C \ Γ 1 , where Γ 1 is any curve starting at ω 1 and tending to ∞ (for example Γ 1 = ω 1 + iR ≥0 ; however as we will see later, sometimes other choices are more convenient). Indeed, replacing D by in the last step of the proof of Theorem 3.1 and noting that this domain is simply connected, analyticity of h follows (of course the expansion (27) is not valid on the larger domain). The function h can be identified with a chart of the Riemann structure of the sheaf g −1 (in this identification, its domain is sometimes referred to as Riemann sheet).
Next, we can define another analytic function Since the function z → g(z) − ω 1 has a zero of order 2 at z 1 , we can extend h analytically from {ω 1 + re iθ : r < ε, θ ∈ (π/2, 3π/4)} to {ω 1 + re iθ : r < ε, θ / ∈ [3π/4, π]}. Table 1. Numerical values of the first few numbers of interest as defined in (26). Let the resulting function be denoted byh. Clearly,h takes values in a small neighborhood of z 1 and since there is a branch point at Since the branch point is of order 2, we conclude thath(ω 1 + t) for t as above has to be real, and slightly bigger than z 1 . Thus we can extend and denote this restriced function by h 2 . Then h 2 can be extended analytically to C\(Γ 1 ∪Γ 2 ) by the monodromy theorem as in the proof of Theorem 3.1. Furthermore, we can "glue together" the domains of h 1 and h 2 along Γ since lim x→ω 1 ± h 1 (x+ yi) = lim x→ω 1 ∓ h 2 (x + yi) and the gluing is analytic (see the analyticity ofh and the fact that the branch point is algebraic of order 2). Since g preserves the real line and all branch points are of second order, this construction can be continued inductively, i.e. we can define the Riemann sheet k ≥ 2 by h k and glue together with the sheets k ± 1 along Γ k±1 . This gives a complete elementary description of the sheaf g −1 .
3.3. Visualizing the sheaf. For better visualization of the sheaf constructed in the previous section, we include some figures. We start with listing the numerical values of the first few critical points η k , η 2 k , ω k and |iπ + log |ω k || (the latter will be needed later) in Table 1. Given ω, there may be infinitely many z's so that g(z) = ω. A plot of some of these z's for ω ∈ (ω 1 , 1) are shown in Fig. 3.1 (note that both ω and z are real). The first branch corresponds to z ∈ (z 1 , ∞) when ω > ω 1 . On this branch h 1 (1) = 0. We call this Riemann sheet the principal branch. The second branch corresponds to z ∈ (z 2 , z 1 ), where ω ∈ (ω 2 , ω 1 ), and the k-th branch corresponds to z ∈ (z k , z k−1 ), where ω is in between ω k−1 and ω k .   Schematic representation of the first three Riemann sheets of the function g −1 (x). The blue dots show g −1 (1) on the various sheets. On the main sheet g −1 (1) = 0. This is joined to the second sheet at the branch cut Γ 1 starting at the branch point ω 1 ≈ −0.217. This is joined to the third sheet at the branch cut Γ 2 starting at the branch point ω 2 ≈ 0.128, and so on. There is an infinite number of Riemann sheets joined successively at the branch points ω k → 0 which accumulate towards zero as k → ∞.
former. The explanation for this fact is that the branch point at ω 2 is not "visible" on the main Riemann sheet, and thus the branch point at ω 1 is the dominant singularity around ω = 1.

3.4.
Properties of the Taylor series. Next, we study the coefficients c n of the Taylor expansion (27). The first few coefficients are easily computed by the Lagrange inversion theorem using (25) and Theorem 3.1. In particular, we obtain Writing e y = ω we can also expand h as a series in y = log ω (here, log ω means the principal value of the logarithm). In applications to Asian options pricing [26], the parameter ω is related to the option strike. In practice, it is usual to expand implied volatilities in logstrike, which corresponds to expanding in y = log ω. See for example [11].
Comparing the first few derivatives with (28), the expansion (28)  Explicitly computing high order terms of these expansions is difficult (cf. [24], pages 411-413). However, we can study the asymptotic properties of the sequences c n and d n .
Proof. (i) We first prove (30). The asymptotics of the coefficients c n for n → ∞ depend on the nature of the singularity on the circle of convergence. The leading asymptotics of h(ω) around this singularity can be translated into asymptotics for the coefficients of the series expansion using the "transfer results" (see Flajolet and Odlyzko [7,8]). Theorem VI.1 and VI.3 of [8] imply that if h is analytic on C \ Γ ′ where Γ ′ = R ≤ω 1 and furthermore as ω → ω 1 , 3 then (30) holds with We defined here h(ω 1 ) = z 1 . With this definition, h is continuous at ω 1 (of course there is no open neighborhood of ω 1 where h is defined). Indeed, Theorem VI.1 of [8] says that the Taylor coefficients a n of the function (1−ω) 1/2 = a n ω n have the large n asymptotics a n = n −3/2 Γ(−1/2) + O(n −5/2 ). The term C 2 (ω − ω 1 ) is a polynomial and so it has a finite Taylor expansion and by Theorem VI.3 of [8], the Taylor coefficients of (1 − ω)(1 − ω) 1/2 are of O(n −5/2 ). We obtain (33) by rescaling the singularity from 1 to ω 1 .
To finish the proof of (30), it remains to verify (32). We will prove that ω 1 is an algebraic branch point of order 2, meaning that near ω 1 , the function h "behaves like a branch of the square root function near 0" . This behavior results from the inversion of the function f ( √ z) around the critical point z 1 . It is known, see for example Theorem 3.5.1 in [30], that the inversion of a non-constant analytic function f (z) = w around a point z 0 with f ′ (z 0 ) = 0 gives a multi-valued function which can be expanded in a convergent series in fractional powers of w − w 0 with w 0 = f (z 0 ), called a Puiseux series. To make these statements mathematically precise, first note that h can be analytically defined on C\Γ ′ exactly as in Section 3.2 (replacing Γ 1 by Γ ′ . With a slight abuse of notation we also denote this function by h.) Then by Theorem 3.5.1 in [30] there is some ε > 0 and an analytic function φ on D ′ := B ε (0) \ R ≤0 so that φ 2 (x) = x (i.e. φ is one branch of the square root function) and another function ψ, analytic on B ε (0) with ψ(0) = z 1 , ψ ′ (0) = 0 so that h(ω − ω 1 ) = ψ(φ(ω)) for all ω ∈ D ′ . Furthermore, the function h has a Pusieux expansion on D ′ , that is where C 1 satisfies C 2 1 = −8z 1 /f ′′ (iξ 1 ). Thus we obtain (32) which completes the proof of (30).
(ii) First we prove that the Taylor expansion (29) has a radius of convergence ρ x := | log |ω 1 | + iπ|. First note that the right hand side of (29) converges to an analytic functioñ h if |ℑy| < π. Indeed, in this case, e y / ∈ R − and soh = h(e y ) is well defined as h is analytic on C \ Γ ′ . Now we claim thath can be analytically extended toD = B ρx (0). To prove this claim, letD + =D ∩ {y : ℑy > π} andD − =D ∩ {y : ℑy < −π}. Let us denote by H the upper half plane, that is H = {y ∈ C : ℑy > 0}. Now observe that the image ofD ∩ H under the exponential function is disjoint to Γ + := ω 1 + e i(2π−ǫ) R ≥0 for ǫ small enough. Thus defining h with the branch cut Γ + , h(e y ) is well defined and analytic forD ∩ H. Likewise, we can analytically extendh toD − using the branch cut Γ − := ω 1 + e iǫ R ≥0 in the definition of h. Clearly, there are singularities at y ± = ±iπ + log |ω 1 | as e y ± = ω 1 . Thus we have verified that the right hand side of (29) converges to an analytic functionh on B ρx (0).
Next, we prove (31). The new feature here is the presence of two singularities on the circle of convergence in the y complex plane, with radius ρ x . This situation is discussed in Chapter VI.5 of [8]. One feature specific to this situation is the oscillatory pattern of the leading asymptotics for the coefficients of the series expansion (29).
We present next some numerical tests of these asymptotic results. In Fig. 3.3 (left) we plot |c k | 1/k vs 1/k for the first 100 terms of the series (28) (left panel) and |d k | 1/k vs 1/k for the first 100 terms of the series (29) (right panel). The red lines in these plots are at 1/|1 − ω 1 | (left) and 1/| log |ω 1 | + iπ| (right). The general trend as k → ∞ is towards the limiting result shown by the red lines is consistent with Proposition 3.1.

Applications: Series expansions
We apply in this section the results of the previous section to two specific expansions appearing in mathematical finance which were discussed in the Introduction.
4.1. The rate function J BS (x). The rate function J BS appears in the small-time asymptotics for the density of the time average of a gBM, and determines the leading short maturity asymptotics for Asian options in the Black-Scholes model. The explicit form of this function for real positive argument x > 0 was given above in equation (20). Complexifying the function J BS (x) and denoting the complex variable by ω, J BS can be written as where h(ω) is the function introduced in Theorem 3.1, and we define Note that J is well defined as the function z → z tanh(z/2) is even.
The series expansion of the rate function J BS (ω) = ∞ n=0 c J,n (ω − 1) n has convergence radius R J = 1. The first few terms are given by and the n → ∞ asymptotics of the coefficients c J,n are given by (ii) The series expansion J BS (e y ) = ∞ n=0 d J,n y n around y = 0 converges for |y| < ρ x = |iπ + log |ω 1 || ≈ 3.49295. The first few terms are The large n asymptotics of the coefficients is with d J ≈ −23.4048 and θ x ≈ 2.02317.
Proof. (i) First, observe that the function J has isolated singularities (simple poles) along the real negative axis at z n = −(2n + 1) 2 π 2 for n ∈ Z, and is analytic everywhere else. The residues of J(z) at the simple poles z n are Res J (z = z n ) = −4z n . Recall that by the function J BS (ω) is analytic at ω. Noting that the only solution of the equation h(ω) = −(2n + 1) 2 π 2 is ω = 0, we conclude that J BS (x) is analytic on B 1−ω 1 (1) \ {0} and has a singularity at 0. Thus R J = 1.
Next, we prove (39). In order to apply the transfer results (cf. Flajolet and Odlyzko) we need to study the behavior of J BS (ω) around the singularity on the circle of convergence. This is a simple pole at ω = 0, so it is sufficient to compute its residue.
The Taylor expansion of h(ω) around ω = 0 is where h(0) = −π 2 and Substituting into (42) gives This means that the residue of J BS (ω) at ω = 0 is equal to 2.
To derive (39), we write the Laurent series of J BS (ω) on a punctured ball around 0 as 2 ω + ∞ n=0 a n ω n . With the notation φ(z) = ∞ n=0 a n ω n , we see that φ is analytic in a neighborhood of 0 and By uniqueness of analytic extension, φ is analytic and (46) holds on the domain B 1−ω 1 (1).
Transferring this result to the asymptotics of the coefficients c J,n yields (39).
(ii) To prove the second part of the proposition, let us write ψ(y) = J BS (e y ). By the first part of the proposition, ψ is analytic as long as h is analytic at e y . Indeed, since e y = 0 for all complex numbers y, the only singularities of J (h(e y )) are due to singularities of h.
The equation e y = ω 1 has infinitely many solutions of the form i(2k + 1)π + log |ω 1 |. Let us write y ± = ±iπ + log |ω 1 | and ρ x = |y ± | for the singularities which are closest to the origin y = 0. Then, exactly as in the proof of Proposition 3.1, we see that ψ is analytic on B ρx (0) and there are branch points at y ± . The picture is now similar to that of Proposition 2.1 (ii) and is different from point (i) discussed above as the singularities closest to the origin are branch points (and are not isolated).
It remains to verify (41). We again use the transfer results (cf. [7]). As in the proof of Proposition 3.1, we define h(ω 1 ) = z 1 . Next, we consider the Taylor expansion of J around z 1 := −η 2 1 (recall that J is analytic at z 1 ). First, we claim that J ′ (z 1 ) = 0. Indeed, for any η with η = tan η, we have η = 2 tan(η/2) 1−tan 2 (η 2 /2) and so Thus we have the Taylor expansion as z → z 1 . Substituting here z − z 1 → h(ω) − h(ω 1 ) and using the Puiseux expansion (32) gives an expansion of the form as ω → ω 1 . The coefficients in this expansion can be obtained in terms of the derivatives of J (z 1 ) and C i in (32) by expanding h(ω) for ω close to ω 1 .
The branch point singularity arises from the second term in (48) with coefficient Finally, changing the expansion variable from ω − ω 1 to y with ω = e y is done in the same way as in the proof of Proposition 3.1(ii). The transfer result now says that the Taylor coefficients A n of the function (1 − z) 3/2 = ∞ n=0 A n z n satisfy whence (41) follows by substituting y ± = ρ x e ±iθx .
The left plot in Figure 4.1 shows |d J,n | 1/n vs 1/n for the first 100 coefficients of the expansion in y n , which illustrates convergence to 1/ρ x , the inverse of the convergence radius in (40).
We study also the approximation error of the asymptotic expansion of the (signed) coefficients (52) ε J,n = d J,n /d asympt where we denoted d asympt J,n the leading asymptotic expression for the coefficients in Eq. (41). The error is formally of order O(n −1 ), so we expect that ε J,n vs 1/n approaches zero linearly.
The right plot in Figure 4.1 shows ε J,n vs 1/n, which confirms the expected general trend towards zero of the error terms. This is seen in more detail in the left plot in Figure 4.2 which zooms into the region close to the horizontal axis.
We note also the presence of large error outliers. For example for n = 97 the relative error is |ε J,97 | = 5.54 such that the subleading correction to the asymptotic result can be comparable with the leading order coefficient. The right plot in Fig. 4.2 shows the error ε n vs cos(θ X (n − 3 2 )) which shows that large errors are mostly associated with small values of the cos factor, which leads to accidental suppression in the leading order contribution.

Asymptotics of the Hartman-Watson distribution. A closed form result for the joint distribution of the time integral of the gBM and its terminal value was obtained by
Yor [34], see (1).
The numerical evaluation of the integral in (2) is challenging for small t ≪ 1. Alternative ways of numerical evaluation of this distribution have been studied in terms of analytic approximations [1,13].
An asymptotic expansion for the Hartman-Watson distribution θ ρ/t (t) as t → 0 at fixed ρ = rt was derived in [28]. The leading order result is expressed by (3) and the real functions    and λ ∈ (0, π) is the solution of the equation The constraints on κ and λ ensure that G(ρ) is positive for positive real ρ.
We will study in this section the series expansions of F and G around ρ = 1, which are relevant for the computation of the distribution (4)  2 ). This relation follows from Proposition 6 in [28] and the relation J(a) = 1 4 J BS (a). However it does not seem easy to use this relation to connect the expansions of these two functions, so we treat them separately.
The large n asymptotics of the coefficients is ii) The function G(e −τ ) has the expansion around τ = 0 as G(e −τ ) = ∞ n=0 d G,n τ n . The first few terms are The large n asymptotics of the coefficients is with d G ≈ 0.719253. Both series expansions converge for |τ | < ρ x = | log |ω 1 | + iπ| ≈ 3.49295.
Proof. (i) The function F (ρ) can be written as 2 . Now the situation (and hence the proof) is similar to Proposition 4.1.
The function F (z) has poles at z k = −k 2 π 2 , k ∈ Z and is analytic everywhere else.
In order to find the large n asymptotics of d F,n we study again the asymptotics of F (e −y ) around its branch points at y ± . This requires the Taylor expansion of F (z) around z 1 = −η 2 1 which is a regular point for this function. The analysis closely parallels that in Prop. 4.1(ii). It turns out that F ′ (z 1 ) = 0, just as for J (z) and so the expansion of F (z) around z 1 reads Substituting here the Puiseux expansion (32) gives that the leading singularity is a branch point at ω 1 arising from the second term in the expansion Changing variables to τ = log ω = − log ρ gives a similar Puiseux expansion in fractional powers of τ − τ ± . The asymptotics of the coefficients d F,n follows again by the transfer result and has the form (58) with (ii) The function G(ρ) can be written as G(ρ) = G(h(1/ρ)), where we defined Writing the numerator in this form ensures that G(ρ) is positive for real and positive ρ, which is the case required for the application to the asymptotics of the Hartman-Watson distribution (3). (This condition was ensured in the original expression by imposing the constraints κ ≥ 0 and λ ∈ (0, π).) This condition replaces the even property of F (z) which was used to resolve the sign ambiguity before. The function G(z) has branch points on the negative real axis at z k = −η 2 k where the denominator in (66) vanishes. The branch point closest to the origin z = 0 is at z 1 and will be relevant for our purpose as z 1 is mapped to ω 1 = 1/ρ, the branch point of h.
In the neighborhood of z 1 the function G(z) diverges as The asymptotics of G(h(1/ρ)) around 1/ρ = ω 1 are obtained by substituting the Puiseux series (32) into (67). Denoting x = 1/ρ for simplicity, we have as x → ω 1 This asymptotics can be transferred to the large n asymptotics of the coefficients d G,n as We chose the branches of the (y − y ± ) −1/4 functions such that they add up to a real function along the positive real axis. This introduces a minus sign between the two terms. This reproduces the stated result, with We present next some numerical tests of the asymptotic expansions for F (ρ) and G(ρ). is the first term on the right hand side of (58). The numerical results confirm a general decreasing trend of the approximation error with n as expected.

Application: Asian options pricing in the Black-Scholes model
Here, we provide a simple application of the expansions studied in this paper, namely the pricing of an Asian option in the Black-Scholes model. Assume that the asset price follows a geometric Brownian motion (71) dS t = rS t dt + σS t dW t started at S 0 , or equivalently, by Itô's formula, S t = S 0 e (r−σ 2 /2)t+σWt . Of course, in practice one has to consider the issue of mis-specification of the price process, which introduces model risk.
Asian option prices are given by expectations in the risk-neutral measure C and P stand for call and put options and the subscript A denotes Asian.
Using the self-similar property of the Brownian motion, the Asian option prices can be expressed in terms of expectations of the standardized integral of gBM A (µ) t = t 0 e 2(Bs+2µs) ds as . and the expectations in (75) are expressed in terms of reduced parameters see e.g. [12]. Denote the density of the standardized time-average of the gBM The leading t → 0 asymptotics (4) gives an approximation for this density where n(t) is a normalization integral, defined such as to ensure the normalization condition ii) outside the convergence domain we use the asymptotic expansions derived in points (i),(ii) Propositions 4 and 5 of [28] for F (ρ), G(ρ) as ρ → 0 and ρ → ∞, respectively. This gives the following approximation for the functions F (ρ), G(ρ), defined by  (49) and (58) of [28] and F R (ρ), G R (ρ) denote the ρ → ∞ asymptotic expansions given in equations (50), (59) of [28]. The first few coefficients d F,k , d G,k are tabulated in Table 2. The approximation F (ρ),Ḡ(ρ) can be made arbitrarily precise by increasing the truncation orders N F,G and by including higher order terms in the tail asymptotics.

Numerical examples.
We illustrate the application of the method on the seven test cases introduced in [10] and commonly used in the literature as benchmark tests for Asian options pricing. Table 3 shows the option prices following from the method proposed here for each of these scenarios, comparing with the precise evaluation in Linetsky [20] obtained using the spectral expansion method.
We used an approximation forF (ρ),Ḡ(ρ) keeping N F = N G = 6 terms in the series expansion. Adding more terms has no impact on the numerical values shown. The series Table 2. The first 10 coefficients of the series expansions of F (ρ), G(ρ) in powers of log(1/ρ). The leading terms of these series expansions are d F,0 = π 2 Table 3. Numerical tests for pricing Asian options in the Black-Scholes model on the seven scenarios commonly considered in the literature [10,20]. All scenarios have K = 2.0. The results from the method proposed here are shown in the column C A (K, T ) which are compared with the precise benchmarks from the spectral expansion method from Linetsky [20].  Table 3 is negligible.
For the numerical evaluation we changed the integration variable ρ to z = log ρ in both integrals in (79) and (80). The 2-dimensional integration was performed in Mathematica using NIntegrate with Method -> NewtonCotesRule and MaxIterations -> 100.
The results in Table 3 show good agreement with the precise benchmark values of [20], the difference being less than 1% in all cases. For most scenarios, the impact of including terms beyond quadratic order in the joint distribution (4) is larger than the impact of varying the volatility parameter by ∆σ ∼ 0.2%, which makes them relevant in practical applications where σ is known with precision of this order of magnitude.