The Hartman-Watson distribution was introduced in the context of directional statistics  and was studied further in relation to the first hitting time of certain diffusion processes . This distribution has received considerable attention due to its relation to the law of the time integral of the geometric Brownian motion, see Eq. (3) below .
The Hartman-Watson distribution is given in terms of the function defined as
The normalized function
defines the density of a random variable taking values along the positive real axis, distributed according to the Hartman-Watson law .
The function appears in the law of the additive functional of a standard Brownian motion
. This integral appears also in the distributional properties of stochastic volatility models with log-normally distributed volatility, such as thelog-normal SABR model  and the Hull-White model [13, 14].
An explicit result for the joint distribution ofwas given by Yor 
where the function is given by (1).
A precise evaluation of is required for the exact simulation of the time integral of the geometric Brownian motion , conditional on the terminal value of the Brownian motion . This problem appears for example in the simulation of the SABR model, see . The paper  proposed an exact simulation method by inverting the Laplace transform of .
The Yor formula yields also the density of by integration over . The usefulness of this approach is limited by the difficulty of evaluating numerically the integral in (1) for small , due to the rapidly oscillating factor [3, 5]. Alternative numerical approaches which avoid this issue were presented in [4, 16].
For this reason considerable effort has been devoted to applying analytical methods to simplify Yor’s formula. For particular cases simpler expressions as single integrals have been obtained for the density of by Comtet et al  and Dufresne . See the review article by Matsumoto and Yor  for an overview of the literature.
In the absence of simple exact analytical results, it is important to have analytical expansions of in the small- region. Such an expansion has been derived by Gerhold in , by a saddle point analysis of the Laplace inversion integral of the density .
In this paper we derive the asymptotics of the function at fixed using the saddle point method. This regime is important for the study of the small- asymptotics of the density of following from the Yor formula (3). The resulting expansion turns out to give also a good numerical approximation of for all . The expansion has the form
The leading order term proportional to is given in Proposition 1, and the subleading correction is given in the Appendix.
The paper is structured as follows. In Section 2 we present the leading order asymptotic expansion of the integral at fixed . The main result is Proposition 1. The leading term in this expansion is shown to give a reasonably good approximation for over the entire range of . Surprisingly, this approximation has also the correct asymptotics, with a coefficient which approaches the exact result in the limit. In Section 3 we compare this asymptotic result with existing results in the literature on the small expansion of the Hartman-Watson distribution [3, 12]. In Section 4 we present the application of this expansion to obtain the leading asymptotics of the density of the time averaged geometric Brownian motion . Assuming that the error term in the asymptotic expansion of is uniformly bounded in , we obtain the leading asymptotics of this density in the form . The exponential factor reproduces a known result obtained previously using Large Deviations theory [2, 21]. An Appendix gives an explicit result for the subleading correction.
2. Asymptotic expansion of as
We study here the asymptotics of as at fixed . This regime is different from that considered in , which studied the asymptotics of as at fixed . We derive this asymptotics in this section.
The asymptotics of the Hartman-Watson integral defined in (1) as is
The function is given by
and the function is given by
Here is the solution of the equation
and is the solution of the equation
The subleading correction is where is given in explicit form in the Appendix.
Proof of Proposition 1.
The function can be written as
These integrals have the form with .
We present in detail the asymptotic expansion for of the integral
where we denote for simplicity . The integral is treated analogously.
i) . The saddle points are given by the solution of the equation . There are three saddle points at with the solution of the equation . The second derivative of at these points is , and .
The contour of integration is deformed from the real axis to the contour shown in the left panel of Fig. 2, consisting of three arcs of curves of steepest descent passing through the three saddle points. Along these arcs we have . Denoting , the path is given by
where is the positive solution of the equation .
The integral is a sum of three integrals along each piece of the path. The real part of the integrand is odd and the imaginary part is even under. This follows from noting that we have and . This implies that i) the integral from B to A vanishes because the integrand is odd, and ii) the real parts of the integrals along and are equal and of opposite sign, and their imaginary parts are equal. This gives
Thus it is sufficient to evaluate only . This integral is written as
where we defined . This is expanded around as
with . Noting , this series is inverted as
The second factor in the integrand of (15) is also expanded in as
Substituting here the expansion (17) gives
The inequality implies that are imaginary, while are real.
By Watson’s lemma , the resulting expansion can be integrated term by term. The leading asymptotic contribution to is
The integral is . Substituting into (14) reproduces the quoted result by identifying . Since are real, the term in (20) does not contribute to , and the leading correction comes from the term. This is given in explicit form in the Appendix.
ii) . There are several saddle points along the imaginary axis. We are interested in the saddle point at with the solution of (9). At this point the second derivative of is .
Deform the integration contour into the curve in the right panel of Fig. 2. This is a steepest descent curve , given by , the positive solution of the equation . The contour passes through the saddle point at .
The integral . is the sum of the two integrals on the two halves of the contour. As in the previous case, the sum is imaginary since along the contour is real, and is even in . Noting that , the first term is odd in and cancels out, and only the second (imaginary) term gives a contribution. We have a result similar to (14)
As before, it is sufficient to evaluate only the integral, which is
where we introduced . This diference is expanded around as
which is inverted as (recall )
The integrand is also expanded in as
Substituting here (24) this can be expanded in powers of as
Again similar to the previous case, are imaginary, while are real, such that only the even order terms contribute.
By Watson’s lemma  we can exchange expansion and integration, and we get the leading asymptotic term
which is evaluated exactly as and gives . The subleading term is given in explicit form in the Appendix.
2.1. Properties of the functions and
Let us examine in more detail the properties of the functions and appearing in Proposition 1. We start by studying the behavior of the solutions of the equations (8) and (9) for respectively. For the solution of (8) approaches as and it increases to infinity as . For , the solution of (9) starts at for , and decreases to zero as .
The derivative can be computed exactly
This shows that the minimum of is reached for that value of for which . From (9) we get that this is , and at this point we have .
We can obtain also the asymptotics of for very small and very large arguments.
i) As the function has the asymptotics
ii) As the function has the asymptotics
iii) Around , the function has the expansion
i) As , the solution of the equation (8) approaches as
with . This follows from writing Eq. (8) as
Taking the logs of both sides gives
This can be iterated starting with the zero-th approximation , which gives (33).
Eliminating in terms of using (8) gives
Substituting (33) into this expression gives the quoted result.
ii) As , the solution of the equation (9) approaches . This equation is approximated as
which is inverted as
Substituting into gives the result quoted above.
iii) Consider first the case . As , we have . From (8) we get an expansion . Substituting into (36) gives the result (2). The same result result is obtained for the case, when is approached from above.
We give next the asymptotics of .
i) As the asymptotics of is
ii) As we have
iii) Around the function has the following expansion
i) The asymptotic expansion of for was already obtained in (33). Substitute this into
ii) Eliminating from the expression (7) for for gives
Using here the expansion of from (38) gives the result quoted.
iii) The proof is similar to that of point (iii) in Proposition 2. ∎
3. Numerical tests
The leading term of the expansion in Proposition 1 can be considered as an approximation for for all , by substituting . Consider the approximation for defined as
where is the solution of the equation
The correction to Eq. (45) is of order .
Table 1 shows the numerical evaluation of for and several values of , comparing also with the small expansion in Eq. (45). They agree well for sufficiently small but start to diverge for larger . In the last column of Table 1 we show also the result of a direct numerical evaluation of by numerical integration in Mathematica.
Figure 3 compares the approximation (black curves) against (blue curves) and direct numerical integration (red curves). We note the same pattern of good agreement between and at small , and increasing discrepancy for larger , which explodes to infinity for sufficiently large . The reason for this explosive behavior is the fact that the denominator in Eq. (45) approaches zero as approaches a certain maximum value. For larger than this value the denominator becomes negative and the approximation ceases to exist.
This maximum value is given by the inequality . For this condition imposes the upper bound .
We show in Figure 4 plots of vs at . The vertical scale of the plot is chosen as in Figure 1 (left) of Bernhart and Mai , which shows the plots of for the same parameters, obtained using a precise numerical inversion of the Laplace transform of . The shapes of the curves are very similar with those shown in .
3.1. Asymptotics for and of the approximation
As expected, the asymptotics matches precisely the exact asymptotics of . Although the approximation was derived in the small limit, it is somewhat surprising that it matches also the correct asymptotics of for , up to a multiplicative coefficient, which however becomes exact as .
Small asymptotics . The leading asymptotics of as was obtained in Barrieu et al 
The exponential factor agrees with the asymptotics of the exact result. Also the leading dependence of the multiplying factor reproduces the improved expansion following from .
Large asymptotics . The asymptotics of is given in Remark 3 in Barrieu et al  as
The dependence has the same form as the exact asymptotics (50).
The leading term in this expansion matches precisely the coefficient obtained from . We conclude that the right tail asymptotics of approaches the same form as the exact asymptotics in the limit .
3.2. Error estimates
We examine here the error of the approximation for obtained by keeping only the leading asymptotic term in Proposition 1. Define the normalized error
Using the subleading correction to the asymptotic expansion of Proposition 1, this error term has the leading term
where is given in Eq. (89). The plot of this function is given in the right panel of Figure 6. Using the asymptotics of from Proposition 3 and of from the Appendix, we get that the asymptotic error (54) vanishes as and . From the plot of this function in Figure 6 we see that it is bounded as for all .
Figure 5 shows the error term obtained by numerical evaluation of , for several values of . The results are in good qualitative agreement with the asymptotic result (54), which they approach as decreases. Similar results are obtained for ; we highlighted the region shown as the error is maximal in this region. Numerical stability issues in the evaluation of limited these tests to . The error shows a decreasing trend as decreases. Assuming that this trend continues all the way down to , these numerical tests suggest an error bound uniform in of the form
for all sufficiently small.
4. Application: The asymptotic distribution of
As an application of the result of Proposition 1, we give here the leading asymptotics of the distribution of the time average of the gBM .
The distribution of can be obtained from Yor’s formula (3) by integration over . Introducing a new variable this is expressed as the integral