Sticky Feller Diﬀusions

We consider a Feller branching diﬀusion process X with drift c having 0 as a slowly reﬂecting (sticky) boundary point with a stickiness parameter 1 /µ ∈ (0 , ∞ ) . We show that (i) the process X can be characterised as a unique weak solution to the SDE system


Introduction
This paper is motivated by the question as to whether it is possible to find a closed-form expression for the transition density function of the Feller branching diffusion process X in [0, ∞) with drift c and diffusion coefficient a satisfying 0 < c < a when 0 is a slowly reflecting (sticky) boundary point for X with a stickiness parameter 1/µ ∈ (0, ∞) . The purpose of the paper is to provide an answer to this question.
Recall that the infinitesimal generator IL X of X acts in (0, ∞) by the rule where a > 0 and b, c ∈ IR are constants. We assume throughout that 0 < c < a in which case 0 is a regular (i.e. non-singular ) boundary point. Other cases are well analysed/understood (if c ≤ 0 then 0 is an exit boundary point, and if c ≥ a then 0 is an entrance boundary point) see e.g. [1, p. 134] combined with [10, p. 357] to obtain closed-form expressions for the transition density functions in these cases. The boundary classification in the present paper refers to the one due to Feller [7] (see e.g. [1, pp. 14-17] for a modern exposition).
In his original paper [6, p. 180] Feller found a closed-form expression for the transition density function of X when 0 is an absorbing/killing boundary point (i.e. when the motion of X terminates upon reaching 0 ). Feller's expression includes a modified Bessel function of the first kind having a positive index. Subsequently Molchanov [14, p. 312] found a closed-form expression for the transition density function of X (upon recalling its one-to-one connection to a Bessel process) when 0 is an instantaneously reflecting boundary point (i.e. when X spends no time at 0 with a strictly positive Lebesgue measure). Molchanov's expression includes a modified Bessel function of the first kind having a negative index. The same negative index expression appears in the paper by Cox et al [3, p. 391] where Feller's branching diffusion is used to model the motion of interest rates extending an earlier paper [21] by Vasicek where the Ornstein-Uhlenbeck process is used to this end.
The question of addressing a sticky boundary behaviour of X at zero (as another possible regular boundary behaviour) was raised by Longstaff [13]. Kabanov et al [11] present an answer in which an interest rate process gets 'stuck' at zero over a random interval of time. This boundary behaviour is different from a 'sticky' (slowly reflecting) boundary behaviour of X at zero (in the sense of Feller's boundary classification [7]) where the amount of time that X spends at zero has a strictly positive Lebesgue measure but contains no interval (thus being nowhere dense like a fat Cantor set). Renewed interest in the sticky boundary behaviour of stochastic interest rates in finance is documented in [12,Ch. 1] and [15] (see also [4] and the references therein) where the transition density function of the sticky reflecting Ornstein-Uhlenbeck process (Vasicek model) is expressed as an infinite series (spectral/eigenfunction expansion) involving confluent-hypergeometric/parabolic-cylinder functions, their derivatives, and related zeros. We will see below that this transition density function can be expressed in the closed form by means of a convolution integral involving a new special function which can be calculated using roots of an algebraic equation for the gamma function. This representation of the transition density function extends to all sticky Feller diffusions and to our knowledge has not been established before.
The paper is organised as follows. In Section 2 we use the Itô-McKean construction of sticky diffusions (cf. [8,Section 10]) and show that X can be characterised as a unique weak solution to the SDE system dX t = (bX t +c) I(X t > 0) dt+ 2aX t dB t (1.2) where B is a standard Brownian motion and 0 (X) is a diffusion local time process of X at 0 . In Section 3 we use the Green function of X (i.e. the Laplace transform of its transition density function in the time domain) to determine a (sticky) boundary condition at zero that characterises the transition density function of X as a unique solution to the Kolmogorov backward/forward equation of X . In Section 4 we calculate the limit at zero of the ratio between the first derivative of an increasing fundamental solution (eigenfunction) to the killed generator equation of X with respect to its scale function and the increasing fundamental solution itself. This limit coincides with the speed measure of X evaluated at the singleton consisting of zero. Making use of this expression in Section 5 we determine the Laplace transform of the transition density function of X in the time domain (i.e. the Green function of X ). The latter expression reveals a new convolution identity for sticky laws showing that the stickiness is entirely embodied in a single multiplication factor of the Laplace transform. In Section 6 we introduce a new special function (which reduces to the Mittag-Leffler function when b = 0 ) as the inverse Laplace transform of the multiplication factor. The function can be expressed by means of the complex inversion formula (Bromwich's integral) which in turn can be calculated by the residue theorem using roots of an algebraic equation for the gamma function. In Section 7 we apply Laplace inversion to the Green function of X from Section 5 and show that the transition density function of X can be expressed in the closed form by means of a convolution integral involving the special function from Section 6 and a modified Bessel function of the second kind. Letting µ ↓ 0 (absorption) and µ ↑ ∞ (instantaneous reflection) the closed-form expression for the transition density function of X reduces to the ones found by Feller [6] and Molchanov [14] respectively. In Section 8 we show that the results derived for sticky Feller diffusions translate over to yield closed-form expressions for the transition density functions of (a) sticky Cox-Ingersoll-Ross processes and (b) sticky reflecting Vasicek processes that can be used to model slowly reflecting interest rates.

Stochastic differential equations
In this section we use the Itô-McKean construction of sticky diffusions (cf. [8,Section 10]) and show that the Feller branching diffusion process X in [0, ∞) with drift c having 0 as a slowly reflecting (sticky) boundary point with a stickiness parameter 1/µ ∈ (0, ∞) can be characterised as a unique weak solution to the SDE system where b ∈ IR and 0 < c < a are given and fixed, B is a standard Brownian motion and 0 (X) is a diffusion local time process of X at 0 defined in (2.11) below. We assume in (2.1)+(2.2) that X starts at some x in [0, ∞) . The stochastic integral with respect to B in (2.1) is understood in Itô's sense. We refer to [5] for standard definitions of the weak/strong solutions to SDEs including their uniqueness that we will use throughout.
1. Recall that the infinitesimal generator IL X of X acts in (0, ∞) by the rule and the stickiness of X at 0 is characterised by the fact that the domain D( ) and the following (sticky) boundary condition at 0 holds [19, p. 310]) where s is the scale function of X given by , and m is the speed measure of X given by for α > 0 and z ∈ IR . The number ν := c/a−1 in (2.5) and (2.6) is referred to as the index of X . Note that ν ∈ (−1, 0) due to 0 < c < a . The stickiness of X at 0 is measured by is a given and fixed number. Letting µ ↓ 0 and µ ↑ ∞ we obtain the boundary behaviour of infinite stickiness (absorption) and zero stickiness (instantaneous reflection) of X at 0 respectively. Recall further that the diffusion local time process for t ≥ 0 . It is well known that the limit in (2.9) exists almost surely and that we have for all (bounded) measurable functions f : [0, ∞) → IR and all t ≥ 0 (cf. [9, pp. 174-175]). Moreover, the mapping (t, x) → x t (X) is continuous on IR + × [0, ∞) almost surely (this can be derived by noting that showing that the continuity claim reduces to the case of standard Brownian motion resolved in [20]). Hence from (2.10) we see that and the identity (2.2) is satisfied due to (2.8) above.
2. The problem as to whether a sticky Feller branching diffusion process X arising from (2.3)+(2.4) above can be obtained from an SDE system driven by a standard Brownian motion to our knowledge has not been considered in the literature. The main result of this section can now be stated as follows.
Proof. 1. We show that the system (2.1)+(2.2) has a weak solution. For this, let B 1 be a standard Brownian motion defined on a probability space (Ω 1 , F 1 , P 1 ) and let IF 1 be the natural filtration of B 1 . Consider the stochastic differential equation Then it is well known (see e.g. [10, p. 357]) that the stochastic differential equation (2.12) has a unique strong (non-negative) solution Z which moreover has 0 as an instantaneously reflecting boundary point due to 0 < c < a . These facts can be verified by exploiting the well-known (and easily verified) identity for t ≥ 0 where Y is a squared Bessel process of dimension δ := 2c/a solving (2.14) with Y 0 = y in [0, ∞) (taken equal to z above) and B 2 is a standard Brownian motion (with for t ≥ 0 ) upon recalling that the stochastic differential equation (2.14) has a unique strong (non-negative) solution Y (see e.g. [19, p. 439]) which moreover has 0 as an instantaneously reflecting boundary point due to δ ∈ (0, 2) (see e.g. [19, p. 442]). This means that s Z = s on [0, ∞) and m Z = m on (0, ∞) with m Z ({0}) = 0 in the notation of (2.5) and (2.6) above (recall (2.8) as well).
Following [8, p. 186 & pp. 200-201] consider the additive functional is well defined (finite) for all t ≥ 0 and satisfies the same properties itself. Moreover, since A = (A t ) t≥0 is adapted to IF 1 it follows that each T t is a stopping time with respect to IF 1 , so that T = (T t ) t≥0 defines a time change with respect to IF 1 . The fact that t → T t is continuous and strictly increasing with T t < ∞ for t ≥ 0 (or equivalently A t ↑ ∞ as t ↑ ∞ ) ensures that standard time change transformations are applicable to continuous semimartingales and their stochastic integrals without extra conditions on their sample paths (see e.g. [19, pp. 7-9 & pp. 179-181]) and they will be used in the sequel with no explicit mention. Consider the time-changed process

Ts
for t ≥ 0 . The identities (2.18)-(2.20) motivate the use of a variant of Doob's martingale representation theorem in order to achieve (2.1). For this, take another Brownian motion B 0 defined on a probability space (Ω 0 , F 0 , P 0 ) and let IF 0 denote the natural filtration of B 0 .
Set Ω = Ω 1 ×Ω 0 , F = F 1 ×F 0 , IF = IF 1 T ×IF 0 and P = P 1 ×P 0 . Then (Ω, F, IF, P) is a filtered probability space. Extend all random variables X 1 and X 0 defined on Ω 1 and Ω 0 respectively to Ω by setting X 1 (ω) := X 1 (ω 1 ) and X 0 (ω) := X 0 (ω 0 ) for ω = (ω 1 , ω 0 ) ∈ Ω . Then it is easily seen that B 1 T and B 0 remain (continuous) martingales with respect to IF and B 0 remains a standard Brownian motion on (Ω, F, P) as well (note that B 1 T and B 0 are independent). It follows therefore that the process B defined by for t ≥ 0 which shows that X and B solve (2.1) with for t ≥ 0 where in the second-to-last equality we use (2.19) above. From (2.24) we see that X satisfies (2.2) and this completes the proof of weak existence.
2. We show that uniqueness in law holds for the system (2.1)+ (2.2). For this, we will undo the time change from the previous part of the proof starting with the notation afresh. Suppose that X and B solve (2.1) subject to (2.2). As part of this hypothesis, we know that X and B are defined on a filtered probability space (Ω, F, IF, P) , both X and B are IF -adapted, and B is not only a standard Brownian motion with respect to P but also a martingale with respect to IF . Consider the additive functional Since t → T t is increasing and continuous it follows that its (right) inverse t → A t defined by is finite for all t ∈ [0, T ∞ ) . Note that t → A t is increasing and right-continuous on [0, T ∞ ) . Moreover, since T = (T t ) t≥0 is adapted to IF it follows that each A t is a stopping time with respect to IF , so that A t defines a time change with respect to IF for t ∈ [0, T ∞ ) . Consider the time-changed process is a continuous martingale with respect to IF given by and therefore the same is true for t → M t whenever s > 0 is given and fixed. It follows therefore that M At is a continuous martingale with respect to F At and we have . Using Lévy's characterisation theorem we can therefore conclude that Recalling that the stochastic differential equation (2.12) has a unique strong solution, this shows that Z t for t ∈ [0, T ∞ ) is a Feller branching diffusion process with drift c having 0 as an instantaneously reflecting boundary point. Moreover, using (2.2) we see that from which we find that Letting t ↑ T ∞ and using that the diffusion local time process 0 (Z) of the Feller branching diffusion process Z solving (2.31) is finite at every finite time, we see that This shows that T ∞ = ∞ almost surely and consequently the process Z solves (2.31) for all t ≥ 0 . From (2.35) we see that t → A t is strictly increasing (and continuous) and hence is the proper inverse for t ≥ 0 (implying also that t → T t is strictly increasing and continuous). It follows in particular that A T t = t so that Since Z is a unique strong solution to the stochastic differential equation (2.31), we see from (2.35)-(2.37) that X is a well-determined measurable functional of the standard Brownian motion W . This shows that the law of X solving (2.1)+(2.2) is uniquely determined and the proof of weak uniqueness is complete.
Remark 2. From (2.15)-(2.17) in the proof above we know that the process X from Theorem 1 is Markov (cf. [22]). Moreover, recalling that the process Z solving (2.12) is Feller (see e.g. [19, p. 440] and use (2.13) above) we know that X is strong Markov (cf. [22]). This fact combined with continuity of X , and the fact established in the proof above that X coincides with Z when away from 0 while t → T t = t 0 I(X s > 0) ds is strictly increasing on IR + , show that X is a regular diffusion process (in the sense of Itô and McKean [9]).

Remark 3.
In relation to definition (2.15) in the proof of Theorem 1 it is instructive to observe that the semimartingale local time L 0 (Z) of the process Z is identically equal to zero. Indeed, noting from (2.12) that d Z, Z t = 2aZ t dt we see that for all t ≥ 0 , where in the last equality we use the fact that 0 is an instantaneously reflecting boundary point for Z , implying the claim.

Sticky boundary condition
In this section we use the Green function of X (i.e. the Laplace transform of its transition density function in the time domain) to determine a (sticky) boundary condition at zero that characterises the transition density function of X as a unique solution to the Kolmogorov backward/forward equation of X .
1. Recall that a transition density function p of X with respect to the speed measure m of X is characterised by being valid for all measurable A ⊆ [0, ∞) whenever t > 0 and x ∈ [0, ∞) are given and fixed. It is well known that p can be chosen to be positive, jointly continuous in all the three variables, and symmetric in the two spatial variables (cf. [9, p. 149]).
2. In view of (2.6) we set m(x) = (x ν /a) e b a x for x > 0 and define for t > 0 and x, y ∈ (0, ∞) . From (3.1) we see that f is a transition density function of X with respect to Lebesque measure on (0, ∞) in the sense that for all measurable A ⊆ (0, ∞) whenever t > 0 and x ∈ (0, ∞) are given and fixed.

Recalling (2.3) we know that f solves the Kolmogorov backward equation
with y ∈ (0, ∞) given and fixed, and the Kolmogorov forward equation for t > 0 and y ∈ (0, ∞) where b(y) = by+c and σ(y) = √ 2ay with x ∈ (0, ∞) given and fixed. In (3.5) and (3.7) above δ z denotes a (formal) density function of the Dirac measure at z in (0, ∞) and the weak convergence is understood to hold for the corresponding probability distribution functions. The initial conditions (3.5) and (3.7) are insufficient to determine unique solutions to (3.4) and (3.6) respectively. In the reminder of this section we thus look for a slowly reflecting (sticky) boundary condition at zero (in the space domain) which when combined with a natural boundary condition at infinity (in the space domain) will accomplish this aim.
3. To exploit its symmetry in the two spatial variables we will perform our analysis in terms of the transition density function p satisfying (3.1) above. For this, recall that the Green function G of X is defined by where λ > 0 is given and fixed. It is well known that where the functions ϕ : [0, ∞) → IR and ψ : [0, ∞) → IR are uniquely determined (up to positive multiplicative constants) by and the Wronskian w is a constant (not dependent on x) defined by . Note that each of G , ϕ , ψ , w depends on λ but we will omit this dependence from the notation for simplicity. The condition (3.12) reflects the fact that 0 is a slowly reflecting (sticky) boundary point for X , and the condition (3.15) reflects the fact that ∞ is a natural boundary point for X . Setting τ c = inf { t ≥ 0 | X t = c } it is well known that the following probabilistic representation of the functions ϕ and ψ is valid where λ > 0 is given and fixed. From this representation it is easily seen that ϕ and ψ belong to the domain D(IL X ) of the infinitesimal generator IL X of X as used in (3.12) and (3.15) above respectively. We will see below that the solutions ϕ and ψ to (3.10)-(3.12) and (3.13)-(3.15) can be expressed in terms of modified Bessel functions ( when b = 0 ) or confluent hypergeometric functions ( when b = 0 ). 4. From (3.8) and (3.9) we see that where λ > 0 is given and fixed. Differentiating with respect to s in (3.18), letting x ↓ 0 and making use of (3.12) (with (2.8) above) and (3.10), we find that Similarly, we point out that ∂ 3 refers to a change of the third argument of p , and (∂ 3 p/∂s)(t; x, y) with (∂ 3 p/∂s)(t; x, 0+) are defined analogously for t > 0 and x, y ∈ (0, ∞) . Note that (∂ 2 p/∂s)(t; x, y) = ∂ x p(t; x, y)/s (x) and (∂ 3 p/∂s)(t; x, y) = ∂ y p(t; x, y)/s (y) for t > 0 and x, y ∈ (0, ∞) . These identities can be used to calculate the limits when either x or y from (0, ∞) tends to either 0 or ∞ as needed throughout.
Integrating by parts we find that for y ∈ (0, ∞) . Inserting this expression back into (3.19) we obtain for all λ > 0 . By the uniqueness theorem for the Laplace transform we can conclude that the following sticky boundary condition is satisfied for all t > 0 and all y ∈ (0, ∞) . Using exactly the same arguments, or exploiting the symmetry of p in the two spatial variables directly, we can conclude that the following sticky boundary condition is also satisfied for all t > 0 and all x ∈ (0, ∞) .
5. From (3.8) and (3.9) we see that where λ > 0 is given and fixed. Letting x → ∞ in (3.24) and using that the first limit value in (3.15) equals zero, we see that for all λ > 0 and y ∈ (0, ∞) . Similarly, differentiating with respect to s in (3.24), letting x → ∞ and using that the second limit value in (3.15) equals zero, we see that for all λ > 0 and y ∈ (0, ∞) . By the uniqueness theorem for the Laplace transform we can conclude from (3.25) and (3.26) that the following natural boundary condition is satisfied for all t > 0 and all y ∈ (0, ∞) . Using exactly the same arguments, or exploiting the symmetry of p in the two spatial variables directly, we can conclude that the following natural boundary condition is also satisfied for all t > 0 and all x ∈ (0, ∞) .
6. Recalling from (3.2) that for t > 0 and x, y ∈ (0, ∞) we see that the boundary conditions

Limit at zero
In this section we calculate the limit at zero of the ratio between the first derivative of an increasing fundamental solution (eigenfunction) to the killed generator equation of X with respect to its scale function and the increasing fundamental solution itself. This limit coincides with the speed measure of X evaluated at the singleton consisting of zero. It will enable us to calculate the Green function of X in the next section.
1. Consider the differential equation (3.10)+(3.13) rewritten as (4.1) axy (x) + (bx+c)y (x) − λy(x) = 0 for x > 0 where b ∈ IR and 0 < c < a are given and fixed. The case b = 0 is special and can be reduced to the case of Bessel processes studied in [17]. This is seen by letting b → 0 in (2.13) which yields Z t = Y at/2 for t ≥ 0 upon recalling (2.12) and (2.14) above. A quick analytic way to verify this connection is to observe that if z = z(x) solves (1/2)z (x)+(δ−1)/(2x)z (x)− λ z(x) = 0 where δ = 2c/a ∈ (0, 2) , then y defined by y(x) := z( 2x/a) for x > 0 solves (4.1). Using the general facts exposed following (4.2) in [17] we can therefore conclude that two linearly independent solutions to (4.1) are given by x → (2x/a) −ν/2 I ν (2 xλ/a) and x → (2x/a) −ν/2 K ν (2 xλ/a) where the first function is increasing and the second function is decreasing on (0, ∞) . Recall in these expressions that I ν and K ν denote the modified Bessel functions of the first and second kind that are respectively given by Recall also in all these expressions that ν = c/a−1 ∈ (−1, 0) . The two linearly independent solutions can then be used to find the solutions ϕ and ψ solving (3.10)-(3.12) and (3.13)-(3.15) respectively which in turn yield the Green function G of X given in (3.9) above. Since these calculations are very similar to those exposed in Sections 4 and 5 of [17] we will omit fuller details and instead derive the Green function G of X by first deriving it for b = 0 and then passing to the limit when b → 0 .
2. The case b = 0 is different because no reduction to Bessel processes is possible (notice that the space-time change (2.13) is no longer a simple constant multiple of time). Two linearly independent solutions to (4.1) are given by for a ∈ IR , b ∈ IR\{0} and x ∈ IR . This can be derived using the equation (108) combined with the table 15 on page 225 and the equation (70) on page 220 in [18]. The functions F 1 and F 2 are increasing on (0, ∞) with F 1 (0) = 1 and F 2 (0) = 0 as well as F 1 (∞−) = F 2 (∞−) = ∞ . The former two identities are evident while the latter two follow from the well-known (13.2.23) and (13.2.39) in [16]) to tackle (4.4) and (4.5) when b > 0 . The four identities, combined with the local convexity/concavity of each solution to (4.1) around any point at which its first derivative vanishes, then imply the increase of F 1 and F 2 as claimed. Note also that F 1 and F 2 are non-negative on (0, ∞) .
Having F 1 and F 2 at hand we can then distinguish two distinct cases as follows.
3. Consider first the case when b < 0 . Motivated by (4.4) and (4.5) we set for x > 0 and seek (positive) constants α and β such that ϕ satisfies (3.12) in addition to (3.10) and (3.11) above. Such constants have already be determined in (4.8) so that ψ satisfies (3.15) in addition to (3.13) and (3.14) above (note that ψ(0) > 0 and ψ(∞−) = 0 imply (3.14) by a similar local convexity/concavity argument as above) and the resulting function U is Tricomi's (confluent hypergeometric) function that is defined by (3.15) can be verified using the well-known asymptotic relation for M stated following (4.6) above combined with (2.5) above. 4. Consider next the case when b > 0 . Motivated by (4.4) and (4.5) we set for x > 0 and seek (positive) constants α and β such that ϕ satisfies (3.12) in addition to (3.10) and (3.11) above. Such constants have already be determined in (4.11) so that ψ satisfies (3.15) in addition to (3.13) and (3.14) above (note that ψ(0) > 0 and ψ(∞−) = 0 imply (3.14) by a similar local convexity/concavity argument as above) and the resulting functionŨ is no longer Tricomi's function (4.9) but a function defined by x) for a < 0 , b ∈ IR\Z and x < 0 . Note that (3.15) can be verified using Kummer's transformation combined with the well-known asymptotic relation for M stated following (4.6) above and (2.5) above. The reason for departure from Tricomi's function U to the modified Tricomi's functionŨ is that the latter asymptotic relation for M fails when x → −∞ (see (13.2.23) in [16] for details).

Laplace transform
In this section we make use of the limit at zero from the previous section and explicitly determine the Laplace transform of the transition density function of X in the time domain (i.e. the Green function of X ). The latter expression reveals a new convolution identity for sticky laws showing that the stickiness is entirely embodied in a single multiplication factor of the Laplace transform. We will focus on the multiplication factor itself in the next section.
Proposition 5 (Laplace transform in the sticky case). The Green function G of X from (3.8) can be explicitly expressed when 0 < m({0}) < ∞ as follows.
Adding and subtracting Γ( following F 1 (x) above and tiding up the resulting expression we end up with (5.1) above as claimed.
2. Consider next the case when b = 0 . We can then either pass to the limit by letting b → 0 in (5.1) or recall from (2.13) that Z t = Y at/2 for t ≥ 0 and use the closed-form expression for the Green function of the Bessel process X := √ Y given in (5.2) of [17] when 0 is a slowly reflecting (sticky) boundary point. Both methods yield (5.2) with (5.3) above as claimed. As this verification is somewhat lengthy and still straightforward we will omit fuller details. It needs to be noted in the second method that the stickiness parameter changes due to the space-time change applied. More specifically, if X F denotes the sticky Feller diffusion solving (2.1)+(2.2) with b = 0 in the present paper and X B denotes the sticky Bessel process whose square solves (2.1)+(2.2) in [17], and if µ F and µ B denote the (reciprocal of the) stickiness parameter of X F and X B respectively, then it is readily verified using the general identities (2.2) and (2.11) that µ F = aµ B . This identity explains why the constant c in (5.2) of [17] is (seemingly) different from the constant κ given in (5.3) above when a = 2 .
2. Looking at the Laplace transforms (5.1) and (5.4) we see that Laplace inversion needs to be applied to products of Kummer's functions and (modified) Tricomi's functions. Although Laplace inversions of such products do not appear to be tabled it turns out that one can calculate them by reduction to absorbing and instantaneously reflecting boundary behaviour as will be shown below. For this reason, and reversing the historical arrow, we now turn to deriving the Green functions of X when its boundary behaviour at 0 is described by m({0}) = 0 (instantaneous reflection) and m({0}) = ∞ (absorption).
Proposition 6 (Laplace transform in the instantaneously reflecting case). The Green function G of X from (3.8) can be explicitly expressed when m({0}) = 0 as follows.
(a) If b < 0 then G is given by Proof. If b = 0 then ϕ and ψ are given by (4.7)+(4.8) when b < 0 and (4.10)+(4.11) when b > 0 with β = 0 in both cases. This is formally seen from (4.18) since µ = ∞ implies that γ = ∞ so that β = 0 as claimed. The same arguments as in the first part of the proof of Proposition 5 above then show that (5.13) and (5.15) hold as claimed. If b = 0 then passing to the limit as µ → ∞ we see that (5.14) reduces to (5.2) with κ ↑ ∞ as claimed. This completes the proof. (a) If b < 0 then G is given by Proof. If b = 0 then ϕ and ψ are given by (4.7)+(4.8) when b < 0 and (4.10)+(4.11) when b > 0 with α = 0 in both cases. This is formally seen from (4.18) since µ = 0 implies that γ = 0 so that α = 0 as claimed. The same arguments as in the first part of the proof of Proposition 5 above then show that (5.16) and (5.18) hold as claimed. If b = 0 then passing to the limit as µ → 0 we see that (5.17) reduces to (5.2) with κ = 0 as claimed. This completes the proof. 3. We now disclose a revealing underlying structure of the sticky boundary behaviour in relation to the building blocks of absorbing and instantaneously reflecting boundary behaviour. For this, let G s denote the Green function of X when 0 is a slowly reflecting (sticky) boundary point, let G r denote the Green function of X when 0 is an instantaneously reflecting boundary point, and let G a denote the Green function of X when 0 is an absorbing boundary point. Note that Propositions 5-7 provide closed-form expressions for G s , G r and G a respectively. The following reshuffle shows how these closed-form expressions are interrelated (i.e. how G a and G r contribute to building G s ).
Corollary 8 (Convolution identity for sticky laws). In the notation specified above the following identity holds where the function L = L(λ) is defined by for λ > 0 and the constant κ is given by for b ∈ IR and 0 < c < a given and fixed with ν = c/a−1 ∈ (−1, 0) .

Special function
In this section we focus on the multiplication factor L from (5.20) that appears in the convolution identity (5.19) above. Applying Laplace inversion to L we obtain a new special function F when b = 0 that reduces to the Mittag-Leffler function when b = 0 . The function F can be expressed by means of the complex inversion formula (Bromwich's integral) which in turn can be calculated by the residue theorem using roots of an algebraic equation for the gamma function. Having F we will be able to apply Laplace inversion to the convolution identity (5.19) and disclose the transition density function of X in the next section.
3. Consider next the case when b = 0 . From (5.20) we see that the analogue of the Laplace transform equation (6.1) reads for λ > 0 . Recall that the Mittag-Leffler function is defined by x n Γ(αn+β) for x ∈ IR where α > 0 and β > 0 are given and fixed. Integrating term by term in (6.16), and summing up the resulting Laplace transforms of the power functions, one finds that for λ > 0 and c > 0 (see [17,Section 6]). Setting α = β = ν +1 and c = κ in (6.17) we see from (6.15) that F 0 is given explicitly by for t > 0 . This function is sometimes refereed to as the generalised Mittag-Leffler function (see [17,Remark 10] for further details).

Inverse Laplace transform
In this section we apply Laplace inversion to the Green function of X from Section 5 and show that the transition density function of X can be expressed in the closed form by means of a convolution integral involving the special function F from the previous section and a modified Bessel function of the second kind.
1. Motivated by the convolution identity (5.19) and following the historical arrow, we first turn to deriving the transition density function of X when its boundary behaviour at 0 is described by m({0}) = ∞ (absorption) and m({0}) = 0 (instantaneous reflection) respectively. The result of Corollary 9 below was originally derived by Feller in his equation (6.2) of [6, p. 180] (where the constant 4b 2 is erroneously included) and the result of Corollary 10 below was originally derived by Molchanov in his Remark 1 (with r = 1 & q = 0 ) of [14, p. 312] (up to the space-time change (2.13) above) both using different arguments. (a) If b = 0 then p is given by Proof. From (2.13) we see that for t ≥ 0 where R is a Bessel process of dimension δ = 2c/a ∈ (0, 2) . The Green function of R is given by the first term on the right-hand side of (5.2) in [17]. The transition density function p R of R (with respect to its speed measure) is given by the first term on the righthand side of (6.3) in [17]. This expression is derived by applying the tabled expression for the inverse Laplace transform (6.4) in [17]. From (7.3) we find by a direct calculation that for t > 0 and x, y ∈ [0, ∞) where f X and f R are transition density functions (with respect to Lebesgue measure) of X and R respectively. Recalling (3.2) we know that f X = p X m X and f R = p R m R . Inserting the known expressions for p R and m R to form f R in (7.4), we obtain a closed-form expression for f X which in turn yields a closed-form expression for p X upon recalling the known expression for m X (for m R and m X see (2.6) in [17] and (2.6) above respectively). It is a matter of routine to verify that the closed-form expression for p X obtained in this way coincides with the expression given in (7.1) when b = 0 . Letting b → 0 in (7.1) we get (7.2) and this completes the proof.
Corollary 10 (Inverse Laplace transform in the instantaneously reflecting case). The transition density function p of X when m({0}) = 0 can be expressed as follows.
(b) If b = 0 then p is given by Proof. This can be derived using exactly the same arguments as in the proof of Corollary 9 above if we replace −ν by ν throughout. Note that the tabled expression for the inverse Laplace transform (6.4) in [17] remains valid in this case and this completes the proof.
2. We next turn to deriving the transition density function of X when its boundary behaviour at 0 is described by 0 < m({0}) < ∞ (slow reflection). For this we let p a denote the transition density function of X from Corollary 9 and we let p r denote the transition density function of X from Corollary 10.
Proof. Both identities (7.7) and (7.8) are a direct consequence of the convolution identity (5.19) combined with the results of Corollaries 9 & 10 above and the results established in Section 6 above. This completes the proof.
analogue of (7.9)+(7.10) combined into a single identity for t > 0 and x ∈ [0, ∞) where κ and ν are the same as in (7.8) above. Plots of the transition density function p of X (with respect to its speed measure) and the transition density function f of X (with respect to Lebesque measure) are given in Figures 4 and 5 respectively. Note that unlike the limiting value of f (t; x, 0+) which equals ∞ , the limiting values of p(t; x, 0+) are both strictly positive and finite, whenever t > 0 and x ∈ [0, ∞) are given and fixed (see Figure 6).

Applications
In this section we show that the results derived for sticky Feller diffusions in the previous sections translate over to yield closed-form expressions for the transition density functions of (a) sticky Cox-Ingersoll-Ross processes and (b) sticky reflecting Vasicek processes that can be used to model slowly reflecting interest rates.
1. In the Cox-Ingersoll-Ross model (see [3, p. 390]) one assumes that the stochastic interest rate X is a Feller branching diffusion process with the infinitesimal generator given by (1.1) above. Focusing on the case when 0 < c < a and 0 is a slowly reflecting (sticky) boundary point for X with a stickiness parameter 1/µ ∈ (0, ∞) , we see by the result of Theorem 1 above that X can be characterised as a unique weak solution to the SDE system (2.1)+(2.2) above. Moreover, we see by the result of Theorem 11 above that the transition density function p of X with respect to its speed measure can be expressed in the closed form by (7.7) and (7.8) above depending on whether b in (1.1) is non-zero or zero respectively. Letting µ ↓ 0 (absorption) and µ ↑ ∞ (instantaneous reflection) the closed-form expression for the transition density function p of X reduces to the ones found by Feller [6] and Molchanov [14] as stated in Corollary 9 and Corollary 10 respectively. The remaining cases c ≤ 0 (exit) and c ≥ a (entrance) are addressed following (1.1) above.
2. In the Vasicek model (see [21,Section 5]) one assumes that the stochastic interest rate V is an Ornstein-Uhlenbeck process with the infinitesimal generator acting in IR by the rule (8.1) IL V = (α−βv)∂ v + σ 2 2 ∂ vv where α ∈ IR , β ∈ IR \ {0} and σ > 0 are given and fixed constants. The process V is known to be a unique strong solution to (8.2) dV t = (α−βV t ) dt + σ dB t with V 0 = v in IR . SettingṼ t := V t − α/β we see thatṼ solves dṼ t = −βṼ t dt + σ dB t so that without loss of generality we may and do assume that α = 0 in the sequel (otherwise replace V byṼ throughout). Define the process Z by setting Considering then the time-changed process (8.9) X t = Z T t defined as in (2.17) above for t ≥ 0 , and proceeding like in the rest of the proof of Theorem 1 above, we obtain a Feller branching diffusion process X having 0 as a slowly reflecting (sticky) boundary point with a stickiness parameter 1/µ ∈ (0, ∞) . The process X can be characterised as a unique weak solution to the SDE system (2.1)+(2.2) under (8.5) above. In parallel, recalling (8.6) and (8.9), we see that the root process (8.10) R t := X t = Z Tt = U Tt = |V Tt | for t ≥ 0 is a Vasicek process having 0 as a slowly reflecting (sticky) boundary point. Moreover, using that R and V have the same speed measure off 0 whose density function equals m R (x) = (2/σ 2 ) e −(β/σ 2 )x 2 for x > 0 , and recalling from (2.6) that m X (x) = (x ν /a) e (b/a)x = (1/2σ 2 ) e −(β/σ 2 )x / √ x for x > 0 , it can be easily verified that for t ≥ 0 . As in the proof of Theorem 1 above it can be further verified from (8.7) using (8.10) and (8.11) that R solves the SDE system with R 0 = |v| . The arguments above also show that R can be characterised as a unique weak solution to the SDE system (8.12)+(8.13). From (2.2) and (8.13) we see that the stickiness of R equals twice the stickiness of X . The arguments from the proof of Theorem 5 in [5] show that the SDE system (8.12)+(8.13) is equivalent to the single equation (8.14) dR t = −βR t I(R t > 0)dt + σI(R t > 0)dB t + µ 2 I(R t = 0)dt obtained by incorporating (8.13) into (8.12). Denoting the transition density function of R with respect to its speed measure by q and the transition density function of R with respect to Lebesque measure by g , and recalling that R = √ X , we find that q(t; x, y) = 1 m R (y) g(t; x, y) = 2y m R (y) f (t; x 2 , y 2 ) (8.15) = 2y m X (y 2 ) m R (y) p(t; x 2 , y 2 ) = 1 2 p(t; x 2 , y 2 ) for t > 0 and x, y ∈ [0, ∞) where the transition density function p of X with respect to its speed measure can be expressed in the closed form by (7.7) above with a = 2σ 2 , b = −2β = 0 , c = σ 2 and with x 2 & y 2 in place of x & y respectively.