On the Cumulants of the First Passage Time of the Inhomogeneous Geometric Brownian Motion

: We consider the problem of the ﬁrst passage time T of an inhomogeneous geometric Brownian motion through a constant threshold, for which only limited results are available in the literature. In the case of a strong positive drift, we get an approximation of the cumulants of T of any order using the algebra of formal power series applied to an asymptotic expansion of its Laplace transform. The interest in the cumulants is due to their connection with moments and the accounting of some statistical properties of the density of T like skewness and kurtosis. Some case studies coming from neuronal modeling with reversal potential and mean reversion models of ﬁnancial markets show the goodness of the approximation of the ﬁrst moment of T . However hints on the evaluation of higher order moments are also given, together with considerations on the numerical performance of the method.


Introduction
Stochastic diffusion processes with linear drift and multiplicative noise are often considered both in theory and applications because they constitute a good compromise between adherence to reality and mathematical tractability when used as models of real phenomena. Among these processes, the Geometric Brownian Motion plays a prominent role in particular in the context of financial modeling. Much is known about this process, although a lack of results emerges when dealing with its version with non-zero asymptotic mean, namely the Inhomogeneous Geometric Brownian Motion (IGBM). The IGBM belongs to the class of Pearson diffusions [1,2] but goes under different names according to the field of study. In the interest rates field, it is called the Brennan-Schwartz model [3,4], denoted as the GARCH model when used for stochastic volatility and for energy markets [5], as Lognormal diffusion process with exogenous factors when used for forecasting and analysis of growth [6,7], in real option literature, it goes under the names of Geometric Brownian motion with affine drift [8,9], Geometric Ornstein-Uhlenbeck [10] or meanreverting Geometric Brownian motion [11].
Here we want to address the related first-passage-time (FPT) problem that arises in many applications in which the stochastic process evolves in the presence of a threshold [12]. First-passage properties underlie a wide range of stochastic processes, spanning from diffusions with limited growth, to the dynamics of the spike generation in neurons or the triggering of stock options.
The mathematical study of the FPT problem rarely leads to a simple finding of the distribution of the random variable FPT, since it can be written in a closed-form expression only in a few cases [13]. Most of the time it consists in the finding of integral equations involving the FPT density or in the calculation of its Laplace transform ( see [14] for an extensive review). From these equations, information on the moments of the random variable or numerical evaluation of its density can be obtained.
Unlike its simpler version, the IGBM presents some difficulties when using these mathematical tools. The transition density of the IGBM does not have a practical closedform expression [8,15] and this complicates the use of Volterra integral equations involving the FPT probability density function (pdf). Even the calculation of the moments of FPT is unpractical, despite the closed-form formula for the Laplace transform of the FPT pdf is available in the literature [15]. In fact, the computation requires high order derivatives of a ratio of hypergeometric functions with respect to the first parameter, and to the best of our knowledge, only the expression of the first moment has been obtained with this strategy so far.
Moreover, an exact simulation scheme is not available. The classical Euler-Maruyama and Milstein methods exhibit non-zero probability to cross the entrance or exit boundaries even for arbitrary small time steps, that is the nature of the boundary is not preserved [16]. For this reason, analytical expressions are of primary importance, but a different strategy must be undertaken.
Recently an approach that uses cumulants to obtain FPT pdf and related statistics has been proposed [17,18]. In particular in [18], we propose an approximation of the FPT pdf of a square-root stochastic process by using cumulants and a Laguerre-Gamma polynomial approximation. The method is based on the closed-form expressions for cumulants recovered from the Laplace transform g * (z) of the FPT random variable T, thanks to the use of the algebra of formal power series of z. Cumulants encode many of the statistical properties of a pdf. For example, the third-order cumulant accounts for the skewness of the FPT pdf while the fourth-order cumulant involves the weight of tails in causing dispersion, that is the kurtosis. Moreover, from cumulants, moments of any order can be recovered.
To apply this method, the starting point is the availability of a closed-form expression of g * (z) in terms of power series. Unfortunately, this is not the case for the IGBM, but we exploit a property of the Tricomi function involving an asymptotic expansion. This expansion allows us to write g * (z) as a formal series in powers of z which results to be particularly suited when transformed via the logarithm function. Thus in the case of a strong positive drift, we are able to get an approximation of the cumulants of T of any order. In Section 4.2 we stress that a good approximation of the first cumulant is better achieved when the starting point of the dynamics is close to the threshold suggesting that, in this case, the approximated cumulants could be used to approximate the FPT pdf even for small times t.
We check the goodness of the approximation using realistic sets of parameters coming from two fields of application of the IGBM: neuronal modeling with reversal potential and mean reversion models of financial markets.
The IGBM has been used to model the changes in the membrane depolarization between two consecutive neuronal spikes in models in which the changes in depolarization are state-dependent and an inhibitory reversal potential is present [19][20][21][22][23]. Further, by solving the related FPT problem, the dynamics of the spike generation have been described especially through the instantaneous firing rate that is the reciprocal of the first FPT cumulant [24]. However, it is still undisclosed how the input statistics affect the moments of order higher than one of the neuronal output [25][26][27].
The second application deals with models of financial markets with mean reversion where the study of extreme events, such as defaults, can be studied using the problems of first-passage time or exit time from a region. In this context the IGBM is better know as the GARCH diffusion process or Brennan-Schwartz model and it is characterized by the properties that the changes in the short rate are state-dependent and unlimited excursions of the process are not allowed (see for instance [3,4,28]).

IGBM and the Related First-Passage-Time Problem
The inhomogeneous geometric Brownian motion (IGBM) is the diffusion process with infinitesimal mean and variance: It is described by the stochastic differential equation with θ, σ > 0, µ, v ∈ R and (W t ) t≥0 is a standard Wiener process. Equation (2) is a linear SDE and admits a unique strong solution.
According to the Feller classification of boundaries if µ > v/θ then v is an entrance boundary, i.e., it cannot be reached by X t in finite time, and there is no probability flow to the outside of the interval (c, +∞). We will always consider this case in the following.
We note that, for µ = 0, the solution of Equation (2) corresponds to the well-known Geometric Brownian motion [29]. Although for this last process a large literature is available, less is known about the IGBM. The transition density of the IGBM does not have a practical closed-form expression [8,30] and an exact simulation scheme is not available albeit very recently a splitting numerical method has been proposed [16]. Explicit expressions of the IGBM solution of Equation (2) can be found for instance in [5] or [31]. In particular in [31] the IGBM moments are written in terms of a transformed Brownian motion, obtained using a change of time method.
We want to address the problem of the evolution of the IGBM process X t in the presence of a constant threshold S. In particular, we are interested in the random time T the process X t starting in x 0 < S crosses S for the first time. This random variable goes under the name of first-passage time (FPT) and is defined as An analytical closed-form expression for the probability density function g(t) of T for the IGBM process is not available but its Laplace transform is known [15] where a = and Ψ is the confluent hypergeometric function of the second kind (or Tricomi's function [32]). When evaluated in z = 0, the Laplace transform (4) returns the probability of crossing S. In principle, moments of T (or equivalently cumulants) of any orders can be computed using higher derivatives of g * , when they exist. For the IGBM process, only the first moment is known to have an explicit analytical expression, due to the cumbersome expression of g * . Indeed, from Equation (4), in Ref. [23] the mean FPT has been written using hypergeometric functions p F q as follows (6) is intended as taking the following limit Asymptotic expansions of the mean first passage time around the starting position x 0 or around the boundary points are given in [15].
To the best of our knowledge, explicit expressions for the variance of T or higher-order moments are not available in the literature for the IGBM process. The derivation from the Laplace transform is impractical, involving at least second derivatives of the Tricomi functions with respect to the first parameter.
Often, where analytical methods are impractical, numerical solutions are sought. In terms of complexity, difficulties apply as well in the numerical evaluation of moments of T through the Siegert equation [33], although the stationary distribution W (x) := lim t→∞ f (x, t|x 0 , 0) of X t in the absence of a threshold is known to be a shifted inverse gamma distribution with the following shape, scale and location parameters [23] X ∞ ∼ Inv-Gamma 1 + 2 For all these reasons in the next section we propose a method based on the algebra of formal power series, to compute approximations of the cumulants of T using an asymptotic expansion of the Laplace transform g * . Indeed an important asymptotic expansion of with a n = a(a + 1) · · · (a + n − 1), n ∈ N the rising factorial. As lim |s|→∞ |ε N (a, b, s)| = 0 for fixed N, a, b, the rhs of (9) returns the hypergeometric series Note that the 2 F 0 series is divergent unless a or 1 + a − b are non-positive integers when it reduces to a polynomial. Nevertheless we might write Ψ(a, b, s) = s −a 2 F 0 a, 1 + a − b; − 1 s by using the Borel summation of the 2 F 0 series. In the following we deal with 2 F 0 as it was an exponential formal power series in a. Indeed Lemma 1 in the subsequent section ] the ring of formal power series in the indeterminate a with coefficients in R.

Formal Cumulants
Given a sequence {g k } k≥0 of real numbers, withg 0 = 1, the sequence of its formal cumulants {c k } k≥1 is such that Note that in the ring R[[z]] of formal power series in z, the (exponential) formal power series logg(z) is well defined [35] with Then {c k } k≥1 are also said formal cumulants ofg(z). A polynomial expression of formal cumulants {c k } k≥1 in terms of {g k } k≥1 involves the logarithmic (partition) polyno- where {B k,j } are the partial exponential Bell polynomials [36] B k,j (g 1 , . . . , with the sum taken over all sequences λ 1 , λ 2 , . . . , λ k−j+1 of non negative integers such that In addition, the sequence {g k } k≥1 might be recovered from its formal cumulants {c k } k≥1 by using the complete Bell (exponential) polynomials {Y k } with {B k,j } given in Equation (14) and {g k } k≥1 replaced by {c k } k≥1 .

Remark 1.
Both the k-th logarithmic polynomial (13) and the k-th complete Bell polynomial (15) are special cases of the k-th general partition polynomial [37] G k (y 1 , . . . , y k ; x 1 , . . . , In particular, for y j = (−1) j−1 (j − 1)!, x j =g j , j ≥ 1, we recover the k-th logarithmic polynomial as given in (13) whereas for y j = 1, x j =c j , j ≥ 1, we recover the the k-th complete Bell (exponential) polynomials as given in (15). Moreove if {y j } j≥1 and {x j } j≥1 are the coefficients of two exponential formal power series in z, let's say x j z j j! then the k-th general partition polynomial G k gives the k-th coefficient of the exponential formal power series composition f ( f (z) − 1) [35].
Wheng(z) is the Laplace transform of a pdf g(t), then the logarithmic and the complete Bell polynomials allow us to deal with cumulants and moments related to g(t). In particular, ifg(z) is the Laplace transform of the FPT pdf and the rhs of Equation (12) is its Taylor expansion about 0, theng 0 = 1,g and there exist cumulants of any order {c k (T)}, see for instance [38]. In particular, from Equations (11) and (13), we havẽ Approximations of FPT Cumulants so that, according to the definitions (5), we have 1 + a − b = −(a + 2B 1 ). From (9) and (10), in the range of parameters such that the asymptotic expression (9) holds, the ratio gives an approximation of the FPT Laplace transform g * (z) given in (4). The aim of the subsequent Theorem 1 is to prove that g(z) ∈ R[[z]]. Thus its logarithm has the exponential formal power series representation (11) whose coefficients {c k } k≥1 are the formal cumulants of g(z) in (19). From (17), the sequence {(−1) kc k } k≥1 will give an approximation of the cumulant sequence {c k (T)} k≥1 for the range of parameters such that the asymptotic expression (9) holds. In particular, an approximation of the mean of T is given in Corollary 1.
To this aim, let us start by proving that 2 , with A ∈ R. This result is given in Lemma 1 where [ n j ] denotes the unsigned Stirling numbers of first type and x denotes the integer part of x ∈ R + . Lemma 1. As exponential formal power series in a, the hypergeometric series 2 F 0 (a, −(a + 2B 1 ); A), A ∈ R has the following representation and D (n) Proof of Lemma 1. From (10), we have 2 F 0 (a, −(a + 2B 1 ); A) = ∑ n≥0 a n −(a + 2B 1 ) n A n n! .
To recover the expansion of 2 F 0 (a, −(a + 2B 1 ); A) as formal power series in a, we need to express both a n and −(a + 2B 1 ) n in powers of a. To this aim, recall that where [ n j ] are unsigned Stirling numbers of the first type. Similarly, By expanding the binomial (a + 2B 1 ) j and grouping with respect to the powers of a, the rising factorial −(a + 2B 1 ) n can be rewritten as Note that both a n and −(a + 2B 1 ) n are polynomials of degree n in a and, after some algebra, their product gives the following polynomial of degree 2n in a a n −(a + 2B 1 with the coefficients {D (n) j } given in (22). Note that 0 (B 1 ) = 1. Therefore, plugging (27) in (23), an expansion of 2 F 0 (a, −(a + 2B 1 ); A) in terms of powers of a is The expression of 2 F 0 (a, −(a + 2B 1 ); A) given in (20) follows after some algebra, exchanging the two summations in (28). (19) is the ratio of two formal power series in a = a(z). In the following theorem, we first show that g(z) ∈ R[[z]], then we expand log g(z) as given in (11) in order to get a closed form formula of its formal cumulants {c k } k≥1 . Theorem 1. The formal cumulants {c k } k≥1 of g(z) in (19) are such that where for A = A 1 , A 2 we have with P k the k-th logarithmic polynomial (13), (20), {B n,k } the partial exponential Bell polynomials given in (14) and Proof of Theorem 1. From Lemma 1, we Since a = a(z), let us first expand a in (19) as exponential formal power series in z.
Recalling that [39] (1 from (5) and (18) we have with Thus a has the following formal power series representation in z, namely a = ∑ n≥1 a n (B 1 ) z n n! , where a n (B 1 ) = By using (36), the next step is to compute the coefficients of 2 F 0 (a, −(a + 2B 1 ); A) as exponential formal power series in z. This follows by observing that 2 F 0 (a, −(a + 2B 1 ); A) as given in (20) results to be the composition of the two exponential formal power series f (z; A, B 1 ) = 1 + ∑ n≥1 f n (A, B 1 ) z n n! and f (z; B 1 ) = 1 + ∑ n≥1 a n (B 1 ) z n n! with coefficients given respectively in (20) and (36). Applying the Faà di Bruno's formula (16), we get where and {B n,k } are the partial exponential Bell polynomials (14). From (11) we recover where c * k (A; B 1 ) are given in (30). Thus the formal cumulants {c k } k≥1 in (29) are obtained by observing that log g(z) = log 2 F 0 (a, −(a + 2B 1 ); A 1 ) − log 2 F 0 (a, −(a + 2B 1 ); A 2 ), (40) using (39) with A replaced respectively by A 1 and A 2 and using the linearity property of formal power series. Finally, h n (A; B 1 ) given in (38) can be further simplified as follows.
As corollary, an approximation of the k-th FPT cumulant c k (T) for A 1 , A 2 → 0 isc k (T) such thatc k (T) = (−1) kc k (42) withc k the k-th formal cumulant as given in Theorem 1.

Corollary 1. For A 1 , A 2 → 0 the mean FPT of T is approximated bȳ
Proof of Corollary 1. The first formal cumulant of g(z) in (19) and from (22) Since then and the thesis follows.

Remark 2.
The approximation of the mean FPT can be rewritten in terms of generalized hypergeometric functions as follows: Corollary 2. For A 1 , A 2 → 0 the FPT variance is approximated bȳ Proof of Corollary 2. The second formal cumulant of g(z) in (19) is with h i (A; B 1 ) given in (31) for i = 1, 2. Thus (50) follows from (53) taking into account that Recall that h 1 (A; B 1 ) has been given in (45), from which (51) follows with h 1,j (B 1 ) as given in (52). From (31), we have (32) and B 2,1 (1, 1) = B 2,2 (1) = 1. As in the proof of Corollary 1, we have From (20) we recover where H j−1 are the generalized harmonic numbers. Therefore from (58) we have Replacing (58) and (55) in (54), after some algebra the result follows.

Applications
In this section we consider two case studies coming from different areas of research to show how the problem of the FPT for the IGBM can arise in applied mathematical modeling and how the proposed methodology can be implemented. We limit to the case of the first moment of T mainly for two reasons. The first one is that in this case, we can compare our results with the one existing in the literature. The second one is that this case is simpler from a numerical point of view and the main scope of the present paper is to propose a mathematical approach rather than perform a complete study of the computational properties of the numerical routine. However, we stress that cumulants of any order can be evaluated numerically by implementing a symbolic-numeric procedure using the R package kStatistics [41]. Indeed the computation of the k-th formal cumulant c * k in (30) involves the computation of the k-th logarithmic polynomial (13), which is a special case of the general partition polynomial (16). The routine GCBellPol in the package kStatistics generates general partition polynomials. Thus the k-th logarithmic polynomial as given in (30) can be recovered by using the routine eGCBellPol in the package kStatistics with where h j (A; B 1 ) and B 1 are given in (31) and (18) respectively. The computation of h j (A; B 1 ) in (59) can be traced back to general partition polynomials by setting with a * j (B 1 ) and f j (A; B 1 ) given in (32) and (21) respectively. As f j (A; B 1 ) is the sum of a numerical series, we have performed its computation by truncating the contribution of the remainder term when smaller of a fixed tolerance.

Stochastic Neuronal Model with Inhibitory Reversal Potential
Stochastic diffusion processes have been extensively used to model the evolution in time of the potential across the membrane of neuronal cells [42]. They describe the changes in the membrane depolarization between two consecutive neuronal spikes. The reason why they are so important is that it is generally agreed that neurons transmit information about their synaptic inputs through these spikes. A spike occurs if the membrane voltage reaches a certain threshold value. From a mathematical point of view, this dynamics is studied by solving the related first-passage-time problem [43]. The classical leaky integrate-and-fire model is described by a linear stochastic differential equation of Ornstein-Uhlenbeck type. However several modifications of the model were proposed to include more realistic features like the state-dependence of the the changes in depolarization and to avoid unrealistic properties like the unlimited state space of the process [19][20][21][22]. These alternative diffusion models maintain the linear drift while the additive noise is replaced by a multiplicative one. Among these diffusion processes, also the IGBM was proposed as a model of the neuronal activity to include the presence of an inhibitory reversal potential V I [23,44].
The IGBM leaky integrate-and-fire model is described by an Itô stochastic differential equation of the following type where µ characterizes the synaptic input, θ is called membrane time constant and takes into account the spontaneous voltage decay towards the resting potential (assumed equal to zero here) in the absence of input, W = {W t } t≥0 is a standard Wiener process and x 0 is the starting depolarization. The diffusion coefficient σ(X t − V I ), where V I represents the inhibitory reversal potential, determines the amplitude of the noise and guarantees at the same time that the changes in the membrane potential are smaller if X t approaches V I . It is assumed that neurons express information about their input mainly by means of the average frequency of spikes described usually by the reciprocal value of the interval between consecutive spikes [24,45]. For this reason, the study of the first moment of the FPT is of primary importance.
The range of parameters such that the proposed asymptotic approximation (43) holds can be interpreted in the neuronal modeling context as the presence of a strong excitatory input. This intuition is confirmed in the following figures where we plot the mean FPT for increasing values of µ.
In Figure 1 we compare the mean FPT of the IGBM process obtained with the asymptotic expression (43) (red dashed line) and that given by the exact expression (6) (blue solid line) for S = 1, σ = 0.8, θ = 1.1, x 0 = 0.2, V I = 0. We note that the proposed approximation slightly underestimates the value of E(T), but the agreement of the two curves improves for larger µ. Moreover for values of µ corresponding to strong external inputs the exact expression (6) is unstable, while the approximation is not.
In the following, we consider physiologically realistic parameter values as suggested in [19] and [46], with their unit measures. The resetting potential is chosen equal to zero, i.e., x 0 = 0 mV, the inhibitory reversal potential is V I = −10 mV, and the firing threshold S = 10 mV. The parameter of spontaneous decay towards the resting potential is chosen θ = 5 ms.
To check how sensitive is the asymptotic expression (43) to a change in µ or σ, we show the mean FPT as functions of the input for three different values of the noise amplitude in Figure 2. As expected the mean waiting time before the first spike decreases as µ increases. Moreover, for fixed µ, the mean FPT decreases for increasing σ (see the plot from the top downwards).
The measure units of the neuronal membrane potential are usually expressed in millivolt, while here we use centivolt. It is just a simple rescale of the SDE (60) and does not affect the dynamics. The change of variable is used here as a numerical trick. In fact the series (43) is more computationally stable for small parameter A i .
Finally, we observe that for the same computational accuracy, although non-optimal for small values of µ, the approximated expression (43) is way more stable than the exact expression (6). For instance, in Figure 3 we compare the mean FPT obtained via the two formulas for σ = 0.2 ms −1/2 and other parameters as in Figure 2. Using the software environment R, for equal time and accuracy of the computation we observe that the asymptotic expression (43), for µ sufficiently large, is a more reliable tool for a numerical study of the firing rate.

Mean Reversion Models of Financial Markets
Mean reverting stochastic processes are often used in models of financial markets (see for instance [3][4][5]8,10,11,15,28]). The mean reversion property implies that the process tends to drive the short rate value towards the long term average level favoring upward (downward) movements if the short rate is too low (high). Among these processes the IGBM plays a relevant role: as GARCH diffusion process or Brennan and Schwartz model it is adopted to describe the value of real options, to study interest rates and to describe the underlying asset in the study of option prices in stochastic volatility models.
Here we want to investigate the sensitivity of the approximated mean FPT to a change in the threshold value S and in the instantaneous volatility σ. To avoid cumbersome calculations involving derivatives of hypergeometric functions with respect to the parameters, we display the main features of the approximated FPT through figures, using realistic parameters coming from the option-pricing literature. We consider the model described by the following stochastic differential equation where µ > 0 represents the long-run mean, σ is the instantaneous volatility, W t is a standard Wiener process and λ is related to the speed of reversion towards the mean level. If λ is large, on average it takes less time for the process to move around the long-run mean and more time to move towards points far from it, i.e., more likely the process will not move far away from the long-run mean level. In the following, we consider realistic parameters values as suggested in [15] for Perpetual American Put Options: µ = 3.25, λ = 0.25, x 0 = 3, σ 2 ∈ [0.08, 0.2].
In Figure 4 we show for X t in Equation (61) the asymptotic first FPT cumulant given by Equation (43) as a function of S. The mean FPT for S = 3 is equal to zero being, in this case, S = x 0 , whereas, as expected, the mean time before the first crossing of S increases if the threshold moves away from x 0 . Moreover, we observe that the mean FPT decreases for increasing values of the instantaneous volatility σ. In particular, if the long-run mean level is below S and λ is large, σ is the main factor that can favor the threshold crossing. In Table 1 we show the relative error between the asymptotic expression (43) and the mean FPT given by (6) in the case σ 2 = 0.2 for five different choices of S. We observe that the asymptotic approximation underestimates slightly the waiting time before the first crossing. The error increases with S, although the relative error grows slowly. We note that the error is smaller for S close to x 0 . The reason is that our approximation performs better for A 1 , A 2 → 0 condition that, for fixed µ, is better fulfilled if x 0 S. Table 1. Relative error between the asymptotic expression (43) and the mean FPT given by (6) for the same choices of parameters of Figure 4 and σ 2 = 0.2.

Conclusions and Open Problems
We addressed the problem of the first passage time T of the IGBM through a constant threshold, giving approximations of the cumulants of T of any order obtained from an asymptotic expansion of the Laplace transform of T using the algebra of formal power series. Moreover from the expression of the cumulants, {c k (T)}, moments of T might be computed by using the recursion formula We have shown that the approximations work better if the IGBM has a large positive drift or if the starting position of the dynamics is close to the threshold. Due to the lack of results in the literature about moments of T of order higher than 1, even though they hold only in a certain parameters range the cumulants approximations constitute a novel result and a step forward in the comprehension of the FPT problem for this diffusion process. The first two approximated cumulants are given explicitly in Corollaries 1 and 2. The higher-order cumulants can be obtained from Theorem 1 analytically or by means of the R symbolic-numeric procedure sketched at the beginning of Section 4 and based on the package kStatistics. We stress that the computational cost of this procedure is lighter than the one for the classical approach with the integro-differential equations, although a more in-depth numerical analysis should be carried out.
These results can be used in many applications ranging from neuronal modeling with reversal potential to mean reversion models of financial markets both for understanding the underlying dynamics of extreme events and possibly for the simultaneous estimation of parameters involved in the model. The examples we have given show that this approximation works better when µ is large or x 0 is close to S.
These last observations suggest using the cumulants (29) and the Laguerre-Gamma polynomial approximation proposed in [18] to get an estimation of the FPT pdf g(t) for the range of parameters for which expression (9) holds. For large t, the asymptotic results from [47] can be applied to obtain information on g(t). Instead, for small-time t, the approximation obtained through the expansion with Laguerre-Gamma polynomials and the cumulants (29) should result to be competitive when compared with the Monte Carlo methods and when traced back to µ large or x 0 close to S. Indeed, the classical simulation schemes like Euler-Maruyama or Milstein method often provide uncorrected values of g(t) around zero and so the approximation of g(t) for small t turns out to be particularly relevant. This is on the agenda of our future research.
As the last remark, we stress that the methodology itself is not limited to the IGBM. As future work, we plan to extend this approach to other processes belonging to the class of Pearson's diffusions, since the expression of the Laplace transform of the FPT pdf for these processes is often written as a ratio of two hypergeometric functions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

IGBM Inhomogeneous Geometric Brownian Motion FPT
First Passage Time pdf probability density function