Another Look at the Hartman-Watson Distributions

The article deals with the Hartman-Watson distributions and presents a new approach to them by investigating a special function u. The function u is strictly related to the distribution of the exponential functional of Brownian motion appearing in the mathematical finance framework. The study of the latter leads to new explicit representations for the function u. One of them is through a new parabolic PDE. Integral relations of convolution type between Hartman-Watson distributions and modified Bessel functions are presented. It turns out that u can be represented as an integral convolution of itself and the modified Bessel function K0. Finally, excursion theory and a subordinator connected to the hyperbolic cosine of Brownian motion are involved in order to obtain yet another representation for u. Possible applications of the resulting explicit formulas are discussed, among others Monte Carlo evaluations of u.


Introduction
The recent paper of Lyasoff [17] has shed a new light on the long studied subject of finding the distribution of an important additive functional of geometric Brownian motion, that is A t = t 0 e 2B u du, where (B t , t ≥ 0) is a real typical Brownian motion, and also on Hartman-Watson distributions. A t has been the subject of the deep studies of Marc Yor and co-authors in several articles; see e.g. [19,22,23] and many others. Such functionals, besides their intrinsic interest, are important for mathematical finance, in particular for pricing financial instruments including the so-called Asian options (see e.g. [7,18,23]). They belong to the Maciej Wiśniewolski m.wisniewolski@mimuw.edu.pl Jacek Jakubowski jakub@mimuw.edu.pl 1 first generation of exotic options and are especially useful for hedging the risk of oscillations of the asset price at the expiration date of the option (see Wystup [27,Sect. 1.5.4]). The difficulties in pricing Asian options come from the payoff function. The latter is an integral functional of the asset price; for instance, for a call option with exercise price K > 0 and expiration date t it is given by (A t − K) + , where A t = 1 t t 0 S u du and S u denotes the price of the basic financial asset at u. Clearly, the task of pricing Asian options reduces to the description of the probability distribution of A t , which is a well known and difficult problem in financial mathematics (see Geman and Yor [9]). In particular, there are known connections between the distribution of A t and the family of Hartman-Watson probability distributions. This paper is devoted to the study of new connections between these objects.
Yor [28] established the distribution of A t by applying the so-called Lamperti relation, which connects geometric Brownian motion, the additive functional A and an appropriate (squared) Bessel process. This resulted in finding the density of the joint distribution of the vector (e 2B t , A t ). Another way of obtaining this distribution was presented by Matsumoto and Yor [22], who applied the method of Feynman-Kac. It has turned out that the distribution of (e 2B t , A t ) is described by the Hartman-Watson (HW) distribution. A random variable has a HW distribution with parameter x > 0 if it has the density η(t, x) = u(t,x) I 0 (x) , t > 0, where u(·, x) is a non-negative function such that for each λ > 0, where I λ denotes the modified Bessel function of index λ ∈ R. HW distributions are not only enigmatic and mathematically interesting distributions, but also a very important object for applications. They appear in the description of the so-called von Mises-Fisher distributions, which are for instance important for statistical inference (see Hartman and Watson [12]). Recall that X r follows the von Mises-Fisher distribution with parameter r > 0 if it has density g X r (x) = 1 2πI 0 (r) e r cos(x) , x ∈ (−π, π) (see e.g. [4]). If X r is such a random variable, then for a two-dimensional Brownian motion W and an independent (of W ) random variable r following the HW distribution with parameter r we have (cos(X r ), sin(X r )) law = W r .
Recall the deep connection of HW distributions with the skew-product representation of planar Brownian motion. Suppose φ t is a continuous determination of the angular component of a complex Brownian motion Z, that is, Z t = |Z t |e iφ t , Z 0 = z = 0. Recalling that η denotes the density of HW distribution and using Eq. 1 we have from [22,Remark 2] for λ ∈ R. Actually, HW distributions were discovered in [12] by combining the distribution of a Brownian diffusion on the n-dimensional sphere with the von Mises-Fisher distribution. However, the partial differential equation (PDE) describing HW distributions, appearing in [12], requires a correction: two terms of the equation should be multiplied by 1 2 . Another interesting object involving HW distributions is the time of explosion of a diffusion given by the stochastic differential equation (SDE) As observed by Kent [14], the time of explosion of X is finite a.s. and follows the HW distribution with parameter x. It is well known that HW distributions are infinitely divisible (Hartman [11]) and belong to the Bondesson class, a special class of Bernstein functions (see Schilling, Song and Vondracek [26,Chapter 9]). Lyasoff [17] has showed that the integral representation of the density of the vector (e 2B t , A t ) obtained by Yor is one of the many possible integral representations (see [17, (5.10a), (5.10b) and Remark 5.1]). In particular, the density of A t conditioned on B t may be interpreted as the contour integral, so the problem of finding the joint distribution of (e 2B t , A t ) can be reduced to rewriting the integral in a tractable form that makes the calculation of residues possible (see [17,Sect. 5]). The idea of Lyasoff was based on the simple form of the Fokker-Planck equation obtained from an application of a random affine equation (see also Bertoin and Yor [5, (17)]). The most important ingredient of HW distributions and of the distribution of (e 2B t , A t ) is a special function u (see Eq. 8 below for its definition). As observed by Barrieu, Rouault and Yor [2], the formula for u obtained by Yor [28] is not numerically tractable and yields significant numerical errors (see Bernhart and Mai [4,Section 3]). This is due to the oscillating nature of the function t → u(t, x) for small values of t. In order to overcome numerical problems, in [4] the Gaver-Stehfest algorithm is applied to obtain u(·, x) via inversion of the Laplace transform. In Section 5 numerical methods for finding HW distributions are presented. Although complicated, the numerical methods are claimed to be efficient (see [4]; for another numerical treatment of HW distributions see also Gerhold [10]).
Our goal is to study the probabilistic nature of HW distributions, to find their connections with modified Bessel functions I and K and finally to provide a new numerically tractable description of HW distributions. The most important results of this work are two new and useful representations of the function u. In what follows we now outline our results obtained. First, we show that and then, we find the following interesting connection of the density g(t, ·) of A t with the function u: Note that, if we assume (without loss of generality) that the interest rate on a financial market is 0 and the asset price is modeled by a geometric Brownian motion, then the price at time 0 of an Asian type instrument with payoff f (A t ) for some positive measurable function f , takes the form x e −xy 2 /2 u(t, y) y dxdy where we have used Fubini's Theorem and defined This is a motivation for developing a computational procedure for integrals of the type ∞ 0 h(y) u(t,y) y dy. We present such a procedure for h belonging to some special class. Possible choices of h will be presented in Example 3.5. However, the procedure is not suitable for finding g(t, ·).
Note that the known formulas for the density of A t , including the recent results of Lyasoff [17], are complicated and difficult to analyze numerically (see also [22,Sect. 4]). We present a new formula for g(t, ·), more suitable for numerical applications, and a new explicit formula for u. First, we show that Although the formula for t (ix) is well known (see Eqs. 12 and 13 below), it seems hardly possible to retrieve a closed formula for t (x) from t (ix). One of the crucial steps in our considerations is showing that which yields Using this and the result of Alili and Gruet [1] we arrive at a new representation of the function u: We will discuss possible applications of the last formula. In particular, comparing the two different representations of u, we prove that (H 1 (t, y)) 4 3 (H 2 (t, y)) and B is a standard Brownian motion. With the use of these representations, both H 1 and H 2 may now be evaluated independently by Monte Carlo methods.
In Section 4, we use the diffusion (Z t ) t = (A t /e B t ) t to derive the PDE satisfied by u. Recall that Z was the subject of extensive studies by Matsumoto and Yor [20,21]. We start from the observation that for a fixed t > 0, P(B t ∈ dx|Z t = z) depends only on z and not on t. A special representation of Z t (Corollary 4.3) enables us to deduce PDEs for the function u and the density of A t (Theorem 4.4). In particular, we retrieve the correct form of the parabolic PDE for u, namely Note that, the parabolic PDE satisfied by u was presented in [12, Thm. 6.1]. However, the result of Hartman and Watson differs by a multiplicative constant from ours. Thus, our method is completely different from that of [12]. We discover that g(t, ·) satisfies the following parabolic PDE: In Section 5 we investigate connections between the modified Bessel function K 0 and the function u. We discover in Theorem 5.1 that, for every t > 0, Next, through the theory of excursions of Brownian motion, we arrive at a time-space convolution identity for u: (that is Theorem 5.2). This yields another implicit description of u (see Theorem 5.3 and Corollary 5.4). For a fixed t > 0, the function φ t : R + → R + given by defines the Lévy measure of a subordinator, and for each λ > 0, Since t → u(t, x) solves an Abel integral equation of convolution type, we can describe u in terms of φ by the formula Next, we prove that the function u defined by is another solution of the parabolic PDE satisfied by u above, and for this solution we are able to produce some identities, analogous to these obtained earlier for u.
We show yet another application of the time-space convolution identity. Using it and the theory of Volterra integral equations we obtain a new description of the modified Bessel functions I n for n ∈ N (Theorem 5.10). It turns out that all I n depend on I 0 , I 0 and on some well defined generating sequences.

Preliminaries
On a complete probability space ( , F, P) endowed with a filtration (F t ) t∈[0,∞) satisfying the usual conditions, consider a standard Brownian motion (B t ) t≥0 . A geometric Brownian motion (GBM) is then defined by Y t = e B t . The classical additive functional A is defined by A t = t 0 e 2B u du. The acclaimed result of Matsumoto and Yor [22,Thm. 4.1] is We denote by g(t, x) the density of A t for fixed t > 0, that is, and we put g(0, x) = 0, for all x.
Recall that the density of a HW distribution is defined, for fixed x > 0, by where u is the function given by To investigate the properties of u we use the process Z defined by Z t = A t /Y t , t ≥ 0. It is well known that Z is a diffusion, and the generator of Z is given by Matsumoto Yor [20,21]). It turns out that the distribution of Z t for any fixed t > 0 can be described in terms of u, namely (see [20,Prop. 4.4]). To investigate the properties of u we apply the properties of Z t and the known fact (see [21,Prop. 1.7]) The following transform of Z t can be obtained from Eq. 10: Indeed, (see also [21,Corr. 12.4]). Some other special identities for HW distributions are known (see [22,Section 4]), and in particular it is well known that for any x ∈ R and t > 0, and Moreover, as mentioned in the Introduction, the basic identity for u is for any λ > 0, where I λ denotes the modified Bessel function with index λ.

The Formula for the Density of A t and a New Representation of u
The goal of this section is to provide new representations for the density of A t and for u. We start from a useful lemma connecting the density of A t and u.

Lemma 3.1 For any t, y > 0 we have
Proof To obtain (15) we integrate (5) with respect to x and we change the variables z = e x /y.
From this lemma we obtain an integral formula which relates the function u, the modified Bessel function K 0 and the hyperbolic sine of the standard Brownian motion. Theorem 3.2 For any t > 0 and α > 0, Moreover, for any t > 0, Proof We know from [22, (3.2)] that for each λ > 0, Taking λ = α 2 /2 we write the last identity in the form where B is a standard Brownian motion. Inserting Eq. 15 for the pdf of A t , denoted by Hence, using Fubini's Theorem and the formula for the density of the GIG distribution (see e.g. [21, (9.1)]), we obtain and Eq. 16 follows. Since α → K 0 ( α 2 + y 2 ) is strictly monotone on (0, ∞) (see [15, (1.11)]), using the monotone convergence Theorem we get Obviously E cos α sinh(B t ) converges to 1 as α ↓ 0, and the theorem is proved.
The next theorem provides a universal method for computing integrals of the form for t > 0 and for f belonging to a broad class described below. Observe that, in the Black-Scholes model the price of an Asian type option with payoff f (A t ), for some function f , is of the form Eq. 19 (see Eq. 3). This follows easily from Lemma 3.1 and Fubini's Theorem. As we will further see (Theorem 3.3 and Corollary 3.4), for f from some special class D, a closed formula for the price of such an option can be obtained. However, the classical payoffs f c (x) = (x − K) + and f p (x) = (K − x) + do not belong to D and another approach has to be worked out for them.
For a measurable function f on (0, ∞) we denote by L −1 f the inverse Laplace transform of f , i.e.,

Theorem 3.3 Let f be a Borel function on [0, ∞) and define g(y)
= e y f (y). Assume that If ϕ is integrable, then Proof Using formula (12), Fubini's Theorem and the fact that ϕ is integrable, we conclude that In order to obtain (21) we need to prove that By the assumption on g we have Recalling that g(y) = f (y)e y we have so to prove (22) it is enough to show that for w ≥ 0, For this, we observe that for fixed , Eq. 23 is a Fredholm integral equation of the first kind, and its unique solution is given by (see [8, 3.3.16]). This finishes the proof.
Now we define a class D of functions f satisfying the assumptions of the last theorem. We say that a function f belongs to D if it is Borel and where for some integrable function ϕ on R + . (21) holds for every f ∈ D.

Corollary 3.4 Formula
Proof Let where ϕ is defined by Eq. 25. By changing variables, for any y ≥ 0, The assertion now follows from Theorem 3.3.
Now we present a few examples of functions f ∈ D.
In general, we can use [3] or [13] to obtain other examples of functions in D.
Next, we provide a new formula for the density of A t . For this, we consider the Fourier transform of the function y → u(t,y) y , integrable for fixed t > 0. We denote it by t : R → C, so recall that Observe that t (x) = ψ t (−ix), for x ≥ 0. Using Eqs. 12 and 13 we can give the closed form of ψ t : Theorem 3. 6 The density of A t for fixed t > 0 is given on (0, ∞) by Re t (sinh z) cosh(z)dz.
Proof Combining (17) and (18) we obtain, for λ > 0, From Fubini's Theorem and the fact that for any x, y > 0, [3, (27), p. 17]), we may write (28) in the form For simplicity of presentation, denote t = Re t . Since we may rewrite (29) as Using a standard trick involving a gamma random variable with parameter 1/2, denoted by Changing the variables w = cosh 2 (z)/(2v) and using Fubini's Theorem yield so the assertion follows.
It turns out that formula (27) leads us to the new representation of the function u which is the most important result of this Section. This corresponds to the discovery of Lyasoff that the form of the density of (e 2B t , A t ), and hence the form of the density of the HW distributions, depends on the choice of parametrization, that is, on the choice of the special contour used for the integral representation of the density of (e 2B t , A t ). In what follows we present another representation of u using a completely different approach compared to the one used in Lyasoff [17].
Proof Fix x > 0. We compare formula (27) for g with formula (4) obtained by Alili and Gruet. As a result, we get Re t (sinh y) cosh(y)dy for each x > 0. Observe that the integral on the RHS can be rewritten as Hence, we conclude that for y > 0, Thus, from the unique trigonometric representation of a complex number, , it is clear that the characteristic function of y → u(t,y) y , i.e. t , is integrable. Thus, from the Fourier inversion Theorem, we obtain, for t, y > 0, dx .
Observe now that the middle integral in the last expression is equal to which finishes the proof.
Proof This follows from Theorem 3.6 and Eq. 32.

Corollary 3.9
The function u has the representation where B is a standard Brownian motion.
Proof The substitution Arc sinh x = z in Eq. 30 and some simple algebra allow one to rewrite u in the form (34).
We have the following upper bound for u. Proposition 3.10 Fix t > 0 and take > 0 such that 1 2 e 8t − 1 ≤ 2 and K = e π 2 /(8t) . Then, for every y > 0, Proof Define Using [13], it can be verified that u(t, y) = 0 for every t > 0 and y ≥ 0. So, Eq. 34 and the elementary identity cos x − cos y = −2 sin( x+y Since by Bougerol's identity, sinh(B 4t ) law = W A 4t , where W and B are independent Brownian motions, we have and the assertion follows from the fact that E A 4t = 1 2 e 8t − 1 .
Comparing two representations of u, namely (8) and (34), yields yet another formula for u. For this let Observe that both H 1 and H 2 can be evaluated in a standard way by Monte Carlo methods.
The following proposition shows how independent evaluations of H 1 and H 2 by Monte Carlo methods can be used to evaluate the function u.
Inserting this directly into Eq. 36 yields the assertion.

Partial Differential Equations
In what follows we apply the distribution of Z t for fixed t (see Preliminaries), to obtain a parabolic PDE satisfied by the function u. To this end, define for fixed t, Lemma 4.1 For fixed t > 0 and any z > 0, Proof Suppose α > 0. We have from Eq. 10, changing the variables w = e x , that where we have used the formula for the density of a GIG distribution (see [21, (9.1)]).

Lemma 4.2 Let W be a standard Brownian motion independent of B and f, g be given by
Eq. 37. For every t ≥ 0 and a positive Borel function h on [0, ∞) we have Proof It is clear that for any λ > 0 and t ≥ 0 we have so if we replace λ with λ/e B t , we see from the definition of Z that Thus, by the properties of Brownian motion and independence of W and B we deduce that Invoking Fubini's Theorem, conditioning on Z u on the LHS, conditioning on Z t on the RHS and using the definitions of f and g, we have Let G t,u (z)dz = P(Z u e W t−u ∈ dz). By Fubini's Theorem, Eq. 38 can now be written as So, from the equality of Laplace transforms we conclude that for any positive Borel function h on [0, ∞), which is our statement.

Corollary 4.3 Assume
where W is a standard Brownian motion independent of the standard Brownian motion B.
Proof It follows from Lemma 4.2 with h = H and the simple fact that E g(Z t ) = e t/2 .
We are now ready to produce a parabolic PDE satisfied by the function u. Before, let us comment on the methodology applied. The original derivation of such a PDE was based on the fact that I λ , which is one of the two basic solutions of the Bessel ODE equation, is the Laplace transform of t → u(t, x), for every fixed x > 0. Using this fact, the PDE for u was found in [12, Thm. 6.1], except that one side of the result in [12] should be multiplied by a constant. Our methodology for proving a PDE satisfied by u is new and is interesting by itself. It is based on the occurrence of u in the density of the diffusion Z and Lemmas 4.1 and 4.2. The second element of the Theorem below is the parabolic PDE for the density g of A t .

Theorem 4.4
The function u given by Eq. 8 and g given by Eq. 6 satisfy the following parabolic PDEs: and Proof We first argue that u is a C 1,2 function on (0, ∞) × (0, ∞). Write where (x, t, k) = e π 2 /(2t) √ 2π 3 t e − k 2 2t sin(π k/t). Fix > 0. We see, by direct computation, that there exists c such that for all t ≥ we have ∂ (x,t,k) ∂t < c , for all x, k ≥ 0. For x ≥ 0, define It is clear that |θ(x)| ≤ K 1 (x), for every x ≥ 0 (see [24, 10.32.9]). Therefore, since is chosen arbitrarily, we conclude from Eq. 42 that ∂u ∂t exists and is continuous. Let be a function such that u(t, x) = ∞ 0 (x, t, k)dk. Fix δ > 0. We see that, for any fixed t > 0, there exists a constant c δ,t such that for i = 1, 2 and all k ≥ 0 and x ≥ δ > 0. Now, we conclude that u belongs to C 1,2 ((0, ∞)× (0, ∞)) and we can differentiate in x under the integral sign. Now, we prove that u satisfies (39) on (0, ∞) × (0, ∞). For this suppose α > 1 and define H (z) = z α , z > 0. Recall that for fixed t > 0 the random variable Z t has the density given by Eq. 9. Using this, Corollary 4.3 and Lemma 4.1 we obtain which after reordering becomes By arguments analogous to those used at the beginning of the proof we may differentiate both sides with respect to t to obtain After eliminating e − (α−1) 2 2 t , this can be rewritten as Observe that Eq. 43 amounts to I t = I I t + I I I t , where we have which after the substitution v = xz and by Fubini's Theorem becomes On the other hand Inserting (45) and (46) and the formula for I I t into Eq. 44 and using the Mellin transform argument shows that a.e. we have We now rewrite this as v and take the derivative with respect to v. We obtain

Rewriting this as
and taking the derivative with respect to v we obtain Since K 1 (x) > 0 for x > 0, the substitution x = 1/v yields By the known identities (x −1 K 1 (x)) = −x −1 K 2 (x) and (x 2 K 2 (x)) = −x 2 K 1 (x) (see [6, p. 638]), the first bracket on the RHS of Eq. 47 is 0, and to simplify the expression in the second bracket we use the formulas As a result we obtain (39). It follows directly from Eq. 8 that u(t, 0) = 0 for every t ≥ 0. The fact that u(0, x) = 0 for every x ≥ 0 is proved in [2, Thm. 1]. Thus Eq. 40 holds.
x . We conclude from Eq. 5 that for any positive From this identity and Fubini's Theorem we obtain Hence, for any λ > 0, On the other hand, so using Fubini's Theorem and a Laplace transform argument we conclude that Taking the derivative with respect to x (here we use the integrability of the function z → u(s, z)z) we obtain It is clear from Eqs. 15, 12 and 13 that we can differentiate both sides with respect to t > 0 to obtain Now we retrieve both integrals on the RHS of Eq. 48 to obtain the desired result (41). Towards this end we use the fact that for any fixed x > 0 and n ∈ N we have ∞ 0 u(t, z)z n e − z 2 2 x dz < ∞, so we can take the first derivative with respect to y in Eq. 15. As a result, for x > 0 we have Using again (49) we can take the x-derivative of Eq. 50 to conclude that Putting together (48), (50) and (51), after some tedious algebra we get (41). This finishes the proof of the theorem.
It turns out that we can also derive a PDE satisfied by the function η given by Eq. 7. Recall that η, for fixed x > 0, is the density of a HW distribution. Observe that the RHS of the next equation is the generator of the diffusion X given by Eq. 2.

Corollary 4.5 The function η satisfies on
Proof Since by definition η(t, x) = u(t, x)/I 0 (x), the corollary follows from Theorem 4.4, the identities (see [6, p. 638]) and some algebra, analogous to the end of the proof of Theorem 4.4.
To end this section we give a simple proof of the following well known result.

Corollary 4.6 The Laplace transform of u is a solution of the Bessel PDE.
Proof We use Theorem 4.4. For λ > 0 we have Integration by parts and u(0, x) = 0 enable writing the LHS in the form while the arguments in the first part of the proof of Eq. 39 of Theorem 4.4 show that the RHS is Denoting U(x) = t 0 e − λ 2 2 t u(t, x)dx and comparing the RHS and LHS we finally obtain

Convolutions
For α ≥ 1 define a Banach space G α of measurable functions on R + by Consider the linear operator T t : G α → G α given by where K 0 is the modified Bessel function. By identity (12) we conclude that for any fixed t > 0 the function x → u(t, x) belongs to G α . We show that T t u(t, x) = u(t, x) for any t > 0, i.e. u(t, ·) is a fixed point of T t .
Proof Using Eq. 12 we see that the derivative of x → ∞ 0 e −y cosh(x) u(t,y) y dy exists and is equal to − sinh(x) ∞ 0 e −y cosh(x) u(t, y)dy. Again by Eq. 12, we conclude that Integration by parts yields If cosh(x) = α (so x = ln(α + √ α 2 − 1)), then using Eq. 12 we see that the last equality takes the form As (see the table of Laplace transforms in [8]), u(t, 0) = 0 and 1 α = ∞ 0 e −αy dy, a Laplace transform argument yields as stated.
Due to the theory of excursions of Brownian motion we present a new implicit description of u in terms of a Lévy measure of some subordinator. This also yields a new recurrent description of the modified Bessel functions I n , n ∈ N. We start with an important result on a time-space convolution identity for u.
Proof Take a positive Borel function f on R. For λ > 0, t > 0 and a Brownian motion B under P, using Eq. 10, we obtain For a positive, bounded and measurable function f on R define a functional F λ,f on U δ by For each element ω of the canonical space there exists an element i 0 (ω) ∈ U such that i 0 (ω)(t) = ω(t) for t < ζ(ω), and i 0 (ω)(t) = 0 for t ≥ ζ(ω). We can now write It is well known that for L the local time at 0 of a Brownian motion, Since (59) and (60) determine the RHS of Eq. 58, and the product of two Laplace transforms is the Laplace transform of the convolution of the functions involved, we finally obtain By the Markov property we have It is clear that I I → 0 as t → ∞. Now, comparing (61) with Eq. 56 and using Fubini's Theorem, and a Laplace transform argument along with the fact that f was chosen arbitrary, we conclude that for any y = 0 and v > 0, By Eq. 9 we rewrite the last equality for y > 0 as which after the substitution x = 1/z on the LHS becomes If γ = cosh(y) and y > 0, then y = ln(γ + γ 2 − 1), and after using Eq. 12 and Fubini's Theorem we obtain Using Eq. 54 and another identity of this kind, (see the table of Laplace transforms [8]), we conclude from Eq. 62 and from a Laplace transform argument that for any x > 0 and v > 0, Hence, from Theorem 5.1 we finally obtain which is Eq. 55.
It turns out that Theorem 5.2 leads to a new implicit description of the function u in terms of the functions and φ 0 ≡ 0. Recall that Eq. 63 is an Abel integral equation of convolution type with respect to u. Therefore the function u is given by Moreover, the function x → φ t (x) is the density of a probability distribution and the Lévy measure of a subordinator.
Proof To prove (65) we first observe that from Eq. 55 we have, for any α > 0, Hence by Eq. 12 and the formula for the Laplace transform of I 0 we have This implies that Letting α tend to ∞ we easily see that c = 1. Now fix an arbitrary λ > 0. Taking α such that λ = cosh(α) − 1 we rewrite the last identity in the form which proves the first statement of the theorem.
The second statement follows from the identity (65). Fix t. If λ goes to ∞ we see that φ t (·) defines the density of a probability distribution. Moreover, the function λ → P(cosh(B t ) ≤ 1 + λ) is a Bernstein function (see [26,Ch. III]) and φ t (·) defines the Lévy measure of some subordinator (see [26, Ch. III, Comments 3.12] ). Theorem 5.3 and properties of exponential random variables easily yield: For the next implicit description of the function u define, for α > 0, two functions on (0, ∞): Observe that for fixed t > 0 the function α is the Mellin transform of x → e −x u(t, x). We have So, again we obtain the well known Abel integral equation of convolution type with respect to α , which is solved uniquely by Eq. 67 (see [8, Volterra integral equations 1.1.5]).
For the next few results denote The function u is the convolution of y → u(t, y)/y and the modified Bessel function K 0 . In Theorem 5.1 we prove that u is the convolution of y → u(t, y)/y and the modified Bessel function I 0 , so u and u have similar structures. From Theorem 4.4 we know that u is a solution of the PDE (39). It turns out that u is another solution of Eq. 39.

Theorem 5.6
The function u defined by Eq. 68 is a solution of Eq. 39.
which is our statement. Now, our goal is to derive for the function u an identity similar to Eq. 12. Let ζ t , for fixed t > 0, be an arcsine random variable on [0, t], that is, a random variable with density P(ζ t ∈ ds) = 1 [0,t] (s) 1 π 1 √ s(t − s) ds.

Proposition 5.7
For any α ∈ R and t > 0, where in the third equality we have used Fubini's Theorem and Eq. 12.
The assertion follows now from the usual Laplace transform arguments.
From Theorem 5.2 and Eq. 14 we have Corollary 5.9 For any λ > 0,