Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion

We provide a comparative analysis of qualitative features of different numerical methods for the inhomogeneous geometric Brownian motion (IGBM). The limit distribution of the IGBM exists, its conditional and asymptotic mean and variance are known and the process can be characterised according to Feller’s boundary classification. We compare the frequently used Euler–Maruyama and Milstein methods, two Lie–Trotter and two Strang splitting schemes and two methods based on the ordinary differential equation (ODE) approach, namely the classical Wong–Zakai approximation and the recently proposed log-ODE scheme. First, we prove that, in contrast to the Euler–Maruyama and Milstein schemes, the splitting and ODE schemes preserve the boundary properties of the process, independently of the choice of the time discretisation step. Second, we prove that the limit distribution of the splitting and ODE methods exists for all stepsize values and parameters. Third, we derive closed-form expressions for the conditional and asymptotic means and variances of all considered schemes and analyse the resulting biases. While the Euler–Maruyama and Milstein schemes are the only methods which may have an asymptotically unbiased mean, the splitting and ODE schemes perform better in terms of variance preservation. The Strang schemes outperform the Lie– Trotter splittings, and the log-ODE scheme the classical ODE method. The mean and variance biases of the log-ODE scheme are very small for many relevant parameter settings. However, in some situations the two derived Strang splittings may be a better alternative, one of them requiring considerably less computational effort than the log-ODE method. The proposed analysis may be carried out in a similar fashion on other numerical methods and stochastic differential equations with comparable features.


Introduction
The inhomogeneous geometric Brownian motion (IGBM), described by the Itô stochastic differential equation (SDE) If the parameter µ = 0, the IGBM coincides with the well-known GBM [14,49], which has often been used as a test equation in the field of stochastic linear stability analysis in the mean-square or almost sure sense [50][51][52]. This theory has been introduced by Mitsui and Saito [52,53], based on stability theory in the sense of Lyapunov [54], and has been extended to systems of SDEs in [55][56][57][58]. Since the standard setting of this approach requires a constant equilibrium solution for which both the drift and diffusion components become zero, it cannot be applied to the IGBM. Nevertheless, known results for the Euler-Maruyama and Milstein schemes applied to the GBM are covered by our study as a special case. Thus, the results presented in this article are also related to stochastic stability analysis for SDEs with inhomogeneous drift coefficients. The paper is organised as follows. In Section 2, we introduce the IGBM and recall its properties. In Section 3, we provide a brief account of the numerical splitting and ODE approaches, introduce the considered numerical schemes for the IGBM and prove that the splitting schemes are mean-square convergent of order 1. In Section 4, we discuss the existence of the limit distribution of the investigated numerical solutions, provide closed-form expressions for their conditional and asymptotic means and variances, analyse the resulting biases and discuss the boundary preservation. In Section 5, we illustrate the theoretical results of Section 4 through a series of simulations, and investigate the required computational efforts of the different numerical methods. In addition, we analyse their ability to approximate the density of the limit distribution and study their behaviour at the lower boundary. Conclusions are reported in Section 6.

The IGBM and its properties
The IGBM is described by the Itô SDE where τ , σ > 0, µ ∈ R and W = (W (t)) t≥0 is a standard Wiener process defined on the probability space (Ω, F , P) with a filtration F = (F(t)) t≥0 generated by W . The initial value Y 0 is either a deterministic non-negative constant or an F (0)-measurable non-negative random variable with finite second moment. Since (1) is a linear and autonomous SDE, a unique strong solution process Y = (Y (t)) t≥0 exists [14,49]. The solution of the homogeneous SDE (if µ = 0) corresponds to the well-known GBM. The solution of the inhomogeneous equation can be expressed in terms of the embedded GBM. In particular, applying the variation of constants formula [49] to (1) yields ) . (2) Limit of the process. Since τ > 0, and thus 1/τ + σ 2 /2 > 0, the limit of the process (2) (for t → ∞), denoted by Y ∞ , exists almost surely and its distribution does not depend on Y 0 [2], see Section 5.3 for a description of the limit distribution when µ > 0. If µ = 0, the process Y ∞ degenerates to the origin 0 almost surely. In the context of linear stability analysis, this means that the equilibrium solution 0 of (1) is asymptotically stable [51].
Conditional and asymptotic mean and variance. Since Y 0 has finite second moment, the mean and variance of the process Y , conditioned on the initial value Y 0 , exist. They are explicitly known [2,3,9] and given by otherwise.
Since τ > 0, from (3), it follows that the asymptotic mean of Y exists. It is given by From (4), it follows that, under the condition σ 2 τ < 2, the asymptotic variance of Y exists. It is given by Boundary properties. Depending on the parameter µ, the IGBM possesses different properties at the boundary 0 according to Feller's boundary classification [41]. In particular, if µ = 0 and Y 0 > 0, the boundary 0 is unattainable and attracting, i.e., the process cannot reach 0 in finite time, but is attracted to it as time tends to infinity. This is in agreement with the fact that Y ∞ ≡ 0 almost surely for µ = 0. In the case that µ = Y 0 = 0, the process is absorbed at the boundary immediately. If µ > 0, then 0 is an entrance boundary, i.e., the process cannot reach the boundary in finite time if Y 0 > 0 or it immediately leaves 0 and stays above it if Y 0 = 0. If µ < 0, the boundary is of exit type, i.e., the process can reach the boundary in finite time and, as soon as it attains the boundary, it leaves [0, +∞) and cannot return into it. In many applications the process is stopped when it reaches an exit boundary, such that its state space is [0, +∞).
Feller's boundary classification is based on the idea of transforming the one-dimensional diffusion into a Wiener process, first by a change of space (through the scale density) and second by a change of time (through the speed density). The scale and speed densities are given by This implies the different types of boundary behaviour explained above, see Table 6.2 in [41]. According to this classification, we define the following properties, which are satisfied by the IGBM:

Numerical methods for the IGBM
. . , N, N ∈ N, t 0 = 0 and t N = t max . We denote byỸ (t i ) a numerical realisation of the process Y at the discrete time points t i = i∆, whereỸ (t 0 ) := Y 0 . Moreover, we denote by ξ i−1 := W (t i ) − W (t i−1 ) ∼ N(0, ∆), i = 1, . . . , N, the Wiener increments which are independent and identically distributed (iid) normal random variables with null mean and variance ∆. In the following, we recall different numerical methods used to generate valuesỸ (t i ) of the IGBM.

Itô -Taylor expansion approach
The most popular approach to derive numerical methods for SDEs is to use appropriate truncations of the Itô-Taylor series expansion [31,59].

Euler-Maruyama and Milstein schemes
Two of the most well-known methods in this class are the Euler-Maruyama and the Milstein schemes. The Euler-Maruyama method yields trajectories of the IGBM through the iteratioñ This method is mean-square convergent of order 1/2. This rate can be increased by taking into account additional terms of the Itô-Taylor expansion. In particular, the Milstein method yields trajectories of the IGBM viã and has a mean-square convergence rate of order 1.

Splitting approach
The second approach we focus on is based on splitting methods [19,24,26,32]. A brief account of their key ideas is provided in the following. Consider an Itô SDE of the form where the drift coefficient and the diffusion component can be expressed as Usually, there are several ways how to decompose the components F and G. The goal is to obtain subequations which can be solved explicitly. Once the explicit solutions are derived, they need to be composed. Two common procedures for doing this are the Lie-Trotter [22] and the Strang [23] approach. Let ϕ [l] t (Y 0 ) denote the exact flows (solutions) of the subequations in (10) at time t and starting from Y 0 . Then, the Lie-Trotter composition of flows and the Strang approach yield numerical methods for (9). The order of the evaluations of the exact flows can be changed, yielding different schemes within each approach.
The first equation, corresponding to the GBM, allows for an exact simulation of sample paths through The second equation is a simple ODE with its explicit solution given by The Lie-Trotter composition yields and the Strang approach results iñ with iid random variables ϕ i−1 , ψ i−1 ∼ N(0, ∆/2), representing independent Wiener increments of length ∆/2.
The proof of Proposition 1 is given in Appendix A, and it follows that of Lemma 2.1 in [60]. We also refer to [33], where this result has been proved for the second Lie-Trotter method (16) and for a more general class of SDEs. Note that, in contrast to the deterministic case [19], the mean-square convergence rate of Lie-Trotter splitting schemes for SDEs cannot be increased by using Strang compositions, i.e., compositions based on fractional ∆/2 steps [60].

ODE approach
An alternative approach to derive numerical solutions of SDEs is to solve properly derived ODEs, a methodology that we briefly recall in the following. Consider the Stratonovich version of (9) given by with G ′ (y) denoting the derivative of G with respect to y. Then, given a fixed time step ∆ > 0 and a Wiener increment ξ i−1 , a numerical solutionỸ (t i ) of SDE (19) can be obtained by defining it as the solution at u = 1 of the ODE This method has been observed to have a mean-square convergence rate of order 1, see, e.g., [38,40], and is called piecewise linear method, since it uses piecewise linear approximations of Brownian paths.
Recently, Foster et al. [40] proposed an extended variant of this approach, using polynomial approximations of Brownian motion. This yielded numerical schemes for SDEs with mean-square order 1.5. In particular, a numerical solutioñ (19) can be obtained by defining it as the solution at u = 1 of the ODE where [·, ·] denotes the standard Lie bracket of vector fields, and the du are rescaled space-time Lévy areas of the Wiener process over the intervals [t i−1 , t i ]. They are shown to have distribution ρ i−1 ∼ N(0, ∆/12) and to be independent of the Wiener increments ξ i−1 . Following the notion in [40], we call this method log-ODE scheme, and we refer to [40] for further details.

Properties of the numerical methods for the IGBM
We now examine the ability of the derived numerical methods to accurately preserve the properties of the process. First, we discuss the existence of the limit distribution of the numerical solutions. Second, we provide closed-form expressions for their conditional and asymptotic means and variances and analyse the resulting biases. Then, we show that the four splitting and the two ODE schemes preserve the boundary properties of the IGBM, while the Euler-Maruyama and Milstein schemes do not.
Formulas (26)-(31) closely resemble a discretised version of (2). Thanks to this, it is possible to show that the limit of the splitting and ODE methods exists for all time steps ∆ > 0 under the same condition as the IGBM. respectively. Since 1/τ + σ 2 /2 > 0, the limit of the numerical solution (for i → ∞), denoted byỸ ∞ , exists almost surely and does not depend on Y 0 , for any choice of the time step ∆ and for all model parameters.
The proof of Proposition 2 is given in Appendix B.

Investigation of the conditional moments
The previously derived relations (24)-(31) also allow for an investigation of the conditional means E[Ỹ (t i )|Y 0 ] and variances Var(Ỹ (t i )|Y 0 ) of the numerical solutions.

Closed-form expressions for the conditional means and variances
In Proposition 3, we provide closed-form expressions of the conditional mean and variance of a general random variable Z i that plays the role of a numerical solutionỸ (t i ) as in (24)-(31) for a fixed time t i . These expressions will allow for a straightforward derivation of the corresponding results for the numerical solutions of interest.

Proposition 3. Consider the real-valued random variable Z i defined by
. . , k, being iid with mean µ x ∈ R and second moment r > 0. The H k+1 , k = 0, . . . , I, are iid with mean µ h ∈ R and second moment r h > 0. Moreover, W k and H k+1 are independent and E[W l W k H k+1 ] = r k µ l−k x p, for k < l and p ∈ R. The mean of Z i conditioned on Z 0 is given by and the variance of Z i conditioned on Z 0 is given by The proof of Proposition 3 is given in Appendix C. Based on Proposition 3, we derive the conditional moments of the Euler-Maruyama, Milstein, splitting and ODE schemes.  (22) and (23), respectively, at time t i = i∆. Their means and variances conditioned on the initial value Y 0 are given by (33) and (34), respectively, with quantities µ x , µ h , r, r h , p, c 1 , c 2 , I, Z 0 and W 0 defined as reported in Table 1.
The proof of Corollary 1 is given in Appendix D.

Remark 3.
To make the results of Proposition 3 and Corollary 1 more approachable, the conditional means of the considered numerical methods are listed in closed-form as follows where L ∆,τ ,σ is defined as with erfi denoting the imaginary error function. The above expressions are obtained from (33) after calculating the geometric sums. Closed-form expressions of the conditional variances can be obtained analogously.
While the conditional means of the Euler-Maruyama and Milstein schemes are equal, their conditional variances are different. This results from the fact that the Milstein scheme takes into account an additional term that is related only to the diffusion coefficient of the SDE. Noting that µ∆ = µτ ∆/τ , it can be observed that the conditional means of the Euler-Maruyama, Milstein and splitting methods depend on ∆/τ and their conditional variances depend on ∆/τ and ∆σ 2 . Remarkably, only the conditional means of the ODE methods depend on σ , while this is not the case for the true conditional mean (3). If µ = 0, the conditional means and variances of the splitting schemes (15)- (18) and ODE schemes (22), (23) coincide with the true quantities (3) and (4), respectively, at time t i .

Remark 4.
Having closed-form expressions for the conditional moments of the numerical solutions allows for a direct control of the simulation accuracy through the choice of the time discretisation step ∆.

Conditional mean and variance biases
Corollary 1 implies that all methods yield conditional means and variances different from the true values. In the following, we study the introduced relative mean and variance biases defined by for each considered numerical method. These biases depend on the time step ∆, the time t i , the initial condition Y 0 and the parameters of the model. While the biases in the conditional means of the ODE methods depend on σ , that of the remaining methods are independent of σ . The biases in the conditional variance depend on all model parameters.
In the top left panel of Fig. 1, we report the relative mean bias (36) in percentage as a function of t i , for Y 0 = 10, ∆ = 0.1, µ = 1, τ = 5 and σ = 0.2. The relative mean biases (in absolute value) introduced by the Strang splitting schemes are significantly smaller than those of the Lie-Trotter splitting schemes and close to 0 for all t i under consideration, with the second Strang scheme performing slightly better than the first one (see the top right panel where we provide a zoom). Moreover, the piecewise linear method performs better than the Lie-Trotter methods, but worse than the Strang schemes.
For the chosen value of σ , the log-ODE method outperforms the Strang methods and produces a bias even closer to 0 for all times t i . However, this fact changes when σ is increased, as shown in the top right panel where we also consider σ = 1 and σ = 1.5. In particular, due to the dependence of the mean of the ODE schemes on σ , they may perform worse than all other methods in terms of preserving the mean when σ increases. Furthermore, it can be observed that at the beginning, the Strang and ODE methods clearly outperform the Euler-Maruyama and Milstein schemes. This changes with increasing time. In particular, the relative mean bias of the Euler-Maruyama and Milstein schemes approaches 0, suggesting an asymptotically unbiased mean (see Section 4.3).
In the bottom left panel of Fig. 1, we report the conditional variance biases (37) in percentage as a function of t i for the same values of Y 0 , ∆, µ, τ and σ = 0.2. All four splitting schemes and both ODE schemes yield better approximations of the conditional variance than the Euler-Maruyama and Milstein schemes for all t i under consideration. The log-ODE method yields again a bias close to 0 from the beginning, outperforming all other methods. This is also the case when σ is increased (figures not shown). Except for t i very small, the Strang schemes outperform the piecewise linear method, and also yield biases close to 0 from the beginning. Moreover, the relative variance biases (in absolute value) of the Lie-Trotter splitting schemes decrease in time and seem to coincide asymptotically with that of the first Strang scheme (see Section 4.3), as it can be observed in the bottom right panel of Fig. 1. Similar results are obtained for other parameter values, time steps and initial conditions. In Fig. 2, we report the relative biases of the conditional mean (36) (top panels) and variance (37) (bottom panels) in percentage as a function of the initial value Y 0 , for t i = 2 and the same parameters as before. All methods introduce larger biases for very small values of Y 0 . This may be explained by the fact that reproducing the features of the process near a boundary, i.e., near 0, is more difficult. For σ = 0.2, the log-ODE method outperforms the other methods, yielding relative biases close to 0 for any considered choice of the initial condition, not being strongly influenced by it. Similar to before, this changes when σ is increased, as illustrated in the top right panel where we also consider σ = 1 and σ = 1.5. The Strang methods (whose mean bias does not depend on σ ) do then introduce the smallest bias in the conditional mean. In general, the performance of the splitting and ODE schemes improves as Y 0 increases, while the Euler-Maruyama and Milstein schemes perform worse for large values of Y 0 . This is in agreement with the fact that Y 0 enters into the conditional means of the splitting and ODE schemes in the same way as in the true quantity, as evident when comparing

Investigation of the asymptotic moments
We now investigate the asymptotic mean and variance of the numerical solutions, i.e., comparing them with the true quantities (5) and (6), respectively.

Closed-form expressions for the asymptotic means and variances
In Proposition 4, we provide closed-form expressions of the asymptotic mean and variance of the random variable Z i introduced in Proposition 3. As before, these relations allow for a straightforward derivation of the corresponding results for the numerical schemes of interest, including necessary conditions that guarantee the existence of the asymptotic quantities. Proposition 4. Let the random variable Z i be defined as in Proposition 3. If |µ x | < 1, the asymptotic mean of Z i is given by If, in addition, r ∈ (0, 1), the asymptotic variance of Z i is given by ) .
The proof of Proposition 4 is given in Appendix E. Based on Proposition 4, we derive the asymptotic moments of the considered numerical schemes.

Corollary 2.
LetỸ (t i ) be the numerical solutions defined through (7), (8), (15)-(18), (22) and (23), respectively. The asymptotic means and variances of the Euler-Maruyama and Milstein schemes are given by The asymptotic means and variances of the splitting schemes are given by The asymptotic means and variances of the ODE schemes are given by where L ≡ L ∆,τ ,σ ,L ≡L ∆,τ ,σ ,L ≡L ∆,τ ,σ , K ≡ K ∆,τ ,σ ,K ≡K ∆,τ ,σ ,K ≡K ∆,τ ,σ are as in (35) (42) is more restrictive than that for the Euler-Maruyama method in (41), agreeing with similar results in the literature [50]. The asymptotic variances of the Lie-Trotter schemes and the first Strang scheme coincide, as previously hypothesised looking at Fig. 1. If µ = 0, the results for the Euler-Maruyama and Milstein methods in Corollary 2 are in agreement with those available in the linear stochastic stability literature for the GBM [51,52]. In particular, the conditions required in (41) and (42) are the same as those guaranteeing their mean-square stability. On the contrary, Corollary 2 implies that the splitting schemes (15)- (18) and the ODE methods (22), (23) are asymptotically first and second moment stable without needing extra conditions.

Asymptotic mean and variance biases
Corollary 2 implies that the derived schemes introduce asymptotic mean and variance biases. In the following, we analyse the resulting asymptotic relative biases with respect to the true quantities (5) and (6), for each considered numerical method. These biases depend on the time step ∆ and on the model parameters. All relative asymptotic biases do not depend on µ. In particular, except for the ODE methods, the asymptotic mean biases depend only on the ratio ∆/τ and their asymptotic variance biases depend on both ∆/τ and ∆σ 2 . As expected, all biases vanish as ∆ → 0, provided that the conditions of Corollary 2 are satisfied.
In the top left panel of Fig. 3, we report the relative biases of the asymptotic mean (53)   yield significantly smaller asymptotic mean biases (in absolute value) than the Lie-Trotter schemes, in agreement with the results reported in the previous section. Moreover, the mean bias (in absolute value) of the second Strang scheme is slightly smaller than that of the first Strang scheme as highlighted in the top right panel, where we provide a zoom. In addition, the log-ODE method introduces a smaller bias in the asymptotic mean than the piecewise linear method. This does not change when considering other values for τ and σ , see the top panels of Fig. 4 and Fig. 5 where we fix ∆ = 0.1 and consider (53) as a function of τ and σ , respectively. Furthermore, for small values of σ , the log-ODE method performs better than the Strang schemes in terms of preserving the asymptotic mean. However, this changes when σ is increased, see the top right panel of Fig. 3 and the top panels of Fig. 5.
In the bottom left panel of Fig. 3, we report the relative biases of the asymptotic variance (54) in percentage as a function of the time step ∆, for τ = 5 and σ = 0.5 fulfilling the conditions of Corollary 2. Note that, the Milstein scheme introduces a larger bias in the variance than the Euler-Maruyama method. Moreover, all splitting schemes yield significantly smaller asymptotic variance biases (in absolute value) than the Euler-Maruyama, Milstein and piecewise linear methods. The log-ODE method, however, outperforms the splitting methods. This fact holds true also for other values of τ and σ , see the bottom panels of Fig. 4 and Fig. 5, where we fix ∆ = 0.1 and plot (54) as a function of τ and σ , respectively. However, there exist combinations of τ and σ for which the condition σ 2 τ < 2 is satisfied and the first Strang and Lie-Trotter methods outperform the log-ODE method in terms of asymptotic variance preservation. This is illustrated in Fig. 6, where we provide a heatmap of for different values of τ and σ . In particular, the region within the white lines corresponds to combinations of τ and σ for which this ratio is smaller than 1, i.e., for which the first Strang and Lie-Trotter methods introduce a smaller relative asymptotic variance bias (in absolute value) than the log-ODE method. This can be also observed in the bottom right panel of Fig. 3, where we compare the log-ODE and splitting methods for σ = 0.5, 0.55, 0.6, with σ = 0.55 within the region marked by the white lines of Fig. 6.
Note also that all biases increase when τ is very small (cf. Fig. 4). This is because all biases depend on the ratio ∆/τ , requiring a small value of the time step ∆ to keep this ratio constant when τ is small. Moreover, for small values of τ the process gets closer to the boundary, where its behaviour is more difficult to preserve. The ODE methods, however, are less deterred by small values of τ , possibly due to their dependence on L ∆,τ ,σ (35),L ∆,τ ,σ (D.2) andL ∆,τ ,σ (D.3). Moreover, the relative biases (in absolute value) of the splitting and ODE methods decrease for large values of τ , while the relative variance biases (in absolute value) of the Euler-Maruyama and Milstein schemes initially decrease and then increase.
Interestingly, while the relative asymptotic variance biases (in absolute value) of the Euler-Maruyama, Milstein and ODE schemes increase in σ (bottom panels of Fig. 5), that of the second Strang method decreases as σ increases (bottom right panel), and that of the first Strang and Lie-Trotter schemes first decreases, and then increases again. The latter scenario is related to the fact that there do exist parameter values for which the first Strang and Lie-Trotter methods introduce a smaller relative variance bias (in absolute value) than the log-ODE method, see Fig. 6.
In general, if both τ and σ are large, such that the limit condition σ 2 τ < 2 is only met tightly, the splitting and log-ODE schemes perform well compared to the Euler-Maruyama, Milstein and piecewise linear method in terms of preserving the asymptotic variance. In particular, even though the Strang and log-ODE schemes perform slightly worse than the Euler-Maruyama and Milstein schemes in terms of the asymptotic mean, they clearly outperform them in terms of the asymptotic variance. For this reason, when, for example, analysing the asymptotic coefficient of variation, i.e., CV( which is a measure of dispersion that allows to simultaneously study the error impinging on both quantities, the Strang and log-ODE schemes are superior to all other schemes. Moreover, if the limit condition σ 2 τ < 2 is only met tightly, the first Strang scheme outperforms the second one in terms of the CV.

Preservation of the boundary properties
As discussed in Section 2, the boundary 0 of the IGBM may be of entrance, unattainable and attracting or exit type, depending on the parameter µ. Corresponding properties motivated by this classification have been introduced at the end of Section 2. A numerical schemeỸ (t i ) is said to preserve these properties if the following discrete versions are fulfilled: • Discrete unattainable property: If µ ≥ 0, then P(Ỹ (t i ) > 0|Ỹ (t i−1 ) > 0) = 1.
It is well known that the Euler-Maruyama and Milstein schemes may fail in meeting such conditions. For example, the Euler-Maruyama scheme (7) does not fulfil the discrete unattainable property for any choice of ∆, since ξ i−1 assumes all values in R with a positive probability [44]. Moreover, the Milstein scheme (8) may not fulfil this property either, if Y (t i−1 )(σ 2 where y :=Ỹ (t i−1 ) and G ′ (y) denotes the derivative of G with respect to y [44]. Thus, to guarantee positivity, the time step ∆ would need to be updated in every iteration step.
Note that, the discrete absorbing and entrance properties are the only properties which are satisfied by the Euler-Maruyama and Milstein schemes, for any time step ∆. In contrast, the ODE methods and the derived splitting schemes preserve the different boundary properties for any choice of time step ∆ > 0, as shown below. Moreover, their boundary behaviour depends only on the parameter µ, as it is the case for the IGBM. The discrete boundary properties can be verified from (15)-(18), (22) and (23), using the corresponding assumptions on the parameter µ and the positivity of the exponential function. A detailed proof of Proposition 5 is given in Appendix F.

Simulation results
We now illustrate the theoretical results introduced in the previous sections through a series of simulations. First, we represent graphically the mean-square convergence order of the different numerical methods and discuss their required computational effort. Second, we focus on the conditional and asymptotic moments. Third, we compare the ability of the different methods to estimate the density of the limit distribution of the process. Finally, we consider the boundary properties, and provide a further investigation of the behaviour of the numerical solutions at the boundary.

Mean-square convergence order and computational effort
The mean-square convergence order of the different numerical methods can be approximated via the root meansquared error (RMSE) considered as a function of the time step ∆. In particular, we define where Y k (t max ) andỸ k (t max ) denote the kth realisation and approximation (obtained under a numerical method using the time step ∆) of the process, respectively, at a fixed time t max .
In the left panel of Fig. 7, we report the RMSEs of the different numerical schemes as a function of the time step ∆ and in log10 scale. We use the same parameter setting as in Figure 4.2 in [40], i.e., we fix n = 10 5 , t max = 5, Y 0 = 0.06, µ = 0.004, τ = 10 and σ = 0.6. Since the IGBM is not known explicitly, the values Y k (t max ) are obtained under the log-ODE method, using the small time step ∆ = 2 −10 . The approximated valuesỸ k (t max ) are produced under the considered numerical methods and for different values of ∆, specifically ∆ = 2 −l , l = 0, . . . , 8. Note that the Y k (t max ) andỸ k (t max ) have to be computed with respect to the same Brownian paths, see [40] and its supporting code for how to deal with the rescaled space-time Lévy areas of a Brownian increment.
As expected, we observe a mean-square convergence rate of order 3/2 for the log-ODE scheme, a rate of order 1 for the Milstein, piecewise linear and splitting methods, and a rate of order 1/2 for the Euler-Maruyama discretisation. The log-ODE method yields the smallest RMSEs, and the Euler-Maruyama method produces the largest error estimates. Among the order 1 methods we observe differences in their accuracies. The first Strang scheme yields the smallest RMSEs, with error estimates slightly smaller that of the piecewise linear method. Moreover, the RMSEs of the two Lie-Trotter and second Strang methods are almost the same, the second Strang method performing slightly worse than the Lie-Trotter schemes. The Milstein method yields the largest error estimates in the considered class of order 1 methods.
These results should be considered in relation to the computational effort required by the different schemes to generate a path, see, e.g., [62]. We measure this effort by counting the number of operations, function evaluations and random numbers required per iteration, i.e., required to produceỸ (t i ) givenỸ (t i−1 ), ∆, τ , µ and σ . This is summarised in Table 2.
The log-ODE method requires the largest computational effort, and the Euler-Maruyama method the slightest. While the effort required by the two Lie-Trotter schemes is the same, the effort of the second Strang method clearly exceeds that of the first. Moreover, while the second Strang scheme yields almost the same RMSEs as the Lie-Trotter schemes (cf. Fig. 7), it requires a greater effort to produce these errors (cf. Table 2).

Conditional and asymptotic moments
Here, we illustrate that the conditional and asymptotic means and variances obtained via numerical simulations are in agreement with the previously derived theoretical expressions. To do so, we define the sample meanm t i and variancê v t i as follows whereỸ k (t i ) denotes the kth simulated value of Y (t i ) under each considered numerical method, respectively. We denote by RE(m t i ) and RE(v t i ) the relative biases (36) and (37), estimated replacing E[Ỹ (t i )|Y 0 ] and Var(Ỹ (t i )|Y 0 ) with the sample mean (58) and variance (59), respectively. To investigate the asymptotic case, we fix t i = 100 and denote by RE(m 100 ) and RE(v 100 ) the relative biases (53) and (54) Table 3. The quantities obtained through numerical simulations are in agreement with the theoretical ones. Moreover, we verified that there is no noteworthy difference in the standard deviations of the estimated values across the different numerical schemes.

Limit distribution
As a further illustration, we investigate numerically the limit distribution of the IGBM. If µ > 0, the limit distribution of the IGBM is an inverse gamma distribution [2,3,9,13] with mean (5) and variance (6). The probability density function of this distribution, which we denote by f Y∞ , is given by where Γ (·) denotes the gamma function, α = 1 + 2/σ 2 τ and β = 2µ/σ 2 . In Fig. 9, we report the true density f Y∞ (60) (grey solid lines) and the densitiesf Y∞ , estimated from n = 10 7 simulated values of Y (100), for µ = 1, τ = 5, σ = 0.55 and Y 0 = 10, using the different numerical schemes. The densities are calculated with a kernel density estimator, i.e., where the bandwidth h is a smoothing parameter and K is a kernel function (here Gaussian). If ∆ = 0.5 (left panels), To quantify the distance between the true and the estimated densities under the considered numerical schemes for different time steps, we consider their Kullback-Leibler (KL) divergences given by where the integral is approximated using trapezoidal integration. The results shown in Fig. 9 are confirmed by the KL divergences (61) reported in Table 3. In particular, the best performance is achieved by the log-ODE method, which yields a very accurate estimate of the density (60), even for ∆ = 1, and even though for τ = 5 and σ = 0.55 the first Strang  Comparison of theoretical quantities with simulated values (n = 10 7 ). We report relative conditional mean and variance biases (36) and (37), asymptotic mean and variance biases (53) and (54)

Boundary properties
An illustration of the preservation of the boundary properties by the splitting and ODE schemes is provided in Fig. 10, where we report trajectories generated with the first Lie-Trotter, first Strang, piecewise linear and log-ODE schemes when the boundary 0 is of entrance (top panel), unattainable and attracting (middle panel) and exit (bottom panel) type.

Crossing probability
As a further illustration of the boundary behaviour, we investigate the probability that the process Y crosses the boundary 0 in a fixed time interval (0, t max ], with t max > 0 and Y 0 > 0. We define as the first passage (hitting) time of Y through 0, and estimate the probability that T < t max as follows where T k denotes the crossing time (62), which is obtained from the kth simulated path of Y and 1 A denotes the indicator function of the set A. We are interested in situations where the process is in a high noisy regime, i.e., it is perturbed by a large noise intensity σ .
Note that the functions obtained under the splitting and ODE schemes lie close to each other, in spite of the large value of σ . When the boundary 0 is of entrance or unattainable and attracting type, it is known that P(T < t max ) = 0 for all values of t max . However, only the splitting and ODE schemes correctly preserve this property, while the Euler-Maruyama method drastically fails for all considered values of ∆ and the Milstein scheme only preserves it for small values of ∆ (left and middle panels). The latter is in agreement with condition (56). Consider, e.g., y = Y 0 = 1. Then ∆ < 5/122 ≈ 0.0402 and ∆ < 5/127 ≈ 0.0394 is required in the entrance or unattainable and attracting case, respectively. In the exit scenario, the probabilities obtained from the Euler-Maruyama and Milstein schemes lie above those obtained from the splitting and ODE schemes. This suggests that the Euler-Maruyama and Milstein methods yield trajectories that exit from [0, +∞) faster than those generated from the other schemes. Similar results are obtained when studying these probabilities as a function of t max for fixed µ. Moreover, independent of the type of boundary behaviour, the crossing probabilities obtained  Table 3. from the Strang splitting and log-ODE schemes seem not to vary significantly as ∆ increases (a few undetected crossings may occur). This suggests their reliability even for large time steps, while those obtained from the Euler-Maruyama and Milstein schemes change for different choices of ∆. The crossing probabilities derived under the Lie-Trotter and piecewise linear methods deviate slightly as ∆ is increased, the latter one performing a bit better.

Conclusion
Any numerical method, constructed to approximate a process of interest, should preserve its qualitative properties. Here, we focus on the IGBM, a process characterised by a constant inhomogeneous term, commonly applied in mathematical finance, neuroscience and other fields. We compare two Lie-Trotter splitting schemes, two Strang splitting schemes and two schemes based on the ODE approach (the classical piecewise linear method [36] and the recently introduced log-ODE method [40]) with the frequently applied Euler-Maruyama and Milstein methods both analytically and via simulations.
We prove that, in contrast to the frequently applied methods, the splitting and ODE schemes preserve the different boundary properties of the IGBM, independently of the choice of the time discretisation step. We also investigate through simulations the probability that the process crosses the lower boundary. Compared to the splitting and ODE schemes, the Euler-Maruyama and Milstein methods suggest not only a positive crossing probability in the entrance or unattainable and attracting case, but also higher crossing probabilities in the exit scenario.
We also prove that, in contrast to the Euler-Maruyama and Milstein methods, the limit distribution of the splitting and ODE schemes exists for all time steps and model parameters. Moreover, we provide closed-form expressions for the conditional and asymptotic means and variances of all considered numerical solutions, and analyse the resulting biases with respect to the true quantities. The Euler-Maruyama and Milstein schemes are the only methods having an  asymptotically unbiased mean (if an extra condition, unrelated to the features of the model, is fulfilled). However, the splitting and ODE schemes yield better approximations of the variance of the process, and do not require extra conditions for the existence of the asymptotic quantities. We observe that the Strang splitting schemes clearly outperform the Lie-Trotter splitting schemes in terms of preserving the mean and density of the limit distribution of the process, and that the log-ODE method performs better than the piecewise linear method throughout. Both the Strang and log-ODE schemes show a solid performance. The biases introduced by the log-ODE method are even smaller than that of the Strang schemes for many relevant parameter configurations. However, the drawback of the log-ODE method is that its mean bias depends on the noise parameter σ , and, consequently, can deteriorate for large values of σ . In this case, the two Strang methods, which perform comparably good throughout, may be better alternatives.
Moreover, we emphasise that the first Strang scheme requires almost the same computational effort as the Lie-Trotter and standard methods, while the second Strang and log-ODE schemes are more computationally expensive. In particular, they require to generate two random numbers in each iteration and rely on more function evaluations.
Guaranteeing a correct behaviour of the simulated process near or at the boundary is important in a variety of applications such as optimal stopping problems or positive asset pricing models. Moreover, having explicit closedform expressions for the first two conditional and asymptotic moments of the numerical solutions may, for example, play an important role in moment based statistical inference [13,63]. All schemes yield biased moments which will effect the inferential approaches. Knowing them explicitly may help to adjust the inferential procedure accordingly. The explicit closed-form expressions also allow for a direct control of the respective simulation accuracy through the time discretisation step. There is a trade-off between computation time and quality of the simulation. To achieve a reasonable computation time, it may be necessary to avoid very small time steps. This becomes particularly important when the numerical method is embedded, for example, in a simulation-based inference method [20,64].
The considered equation, its properties and their analysis are also meant as a contribution to extend the range of qualitative features that characterise the quality of numerical methods. The presented results on the IGBM may be extended to other numerical methods and to a wider class of SDEs with similar features. For example, one may derive the exact moments of numerical solutions of other Pearson diffusions [13] and analyse their boundary behaviour in a similar fashion. The presented analysis may also be extended to multi-dimensional versions of the IGBM, and to a broader class of equations, e.g., via adapted linearisation and diagonalisation procedures. Finally, the construction of a boundary preserving numerical method for the IGBM which has at least an asymptotically unbiased mean, still remains an open problem.
The first-order mean-square convergence of the splitting methods follows by recalling the fact that the Milstein method is mean-square convergent of order 1 and by applying Milstein's fundamental theorem on the mean-square convergence order [65] (see also Theorem 1.1 in [31]). □

Appendix C. Proof of Proposition 3
Proof. Since the underlying X j , j = 1, . . . , k, are iid, the mean and variance of W k , k ∈ N, are given by Using the independence of W k and H k+1 , the mean of Z i conditioned on Z 0 is given by (33).
To compute the variance of Z i conditioned on Z 0 , using the independence of W k and H k+1 , we have that Further, using the independence assumption of W k and H k+1 , we obtain for k < l
The splitting schemeỸ L1 (t i ) (26) can be rewritten as (32) with W k defined via  Table 1. Since ξ j ∼ N(0, ∆), the random variable − , and thus the X j are iid random variables with log-normal distribution, mean µ x and second moment r given by Similarly, the splitting schemesỸ L2 (t i ) (27) andỸ S1 (t i ) (28) can be rewritten as (32) using X j given by (D.1), as forỸ L1 , and the values reported in Table 1.
Further, since ϕ j and ψ j are iid random variables distributed as N(0, ∆/2), we have that ξ j := ϕ j +ψ j ∼ N(0, ∆). Setting using an index shift and splitting off the 0-th element, the splitting schemeỸ S2 (t i ) (29) can be rewritten as (32) with the values reported in The piecewise linear methodỸ Lin (t i ) (22) can be rewritten as (32) with W k defined via X j , j = 1, . . . , k, as in (D.1), and the values reported in Table 1. In particular, the mean of H k+1 is given by (35) ) .