Simulation of sample paths for Gauss–Markov processes in the presence of a reflecting boundary

Algorithms for the simulation of sample paths of Gauss–Markov processes, restricted from below by particular time-dependent reflecting boundaries, are proposed. These algorithms are used to build the histograms of first passage time density through specified boundaries and for the estimation of related moments. Particular attention is dedicated to restricted Wiener and Ornstein–Uhlenbeck processes due to their central role in the class of Gauss–Markov processes. Subjects: Science; Mathematics& Statistics; Statistics& Probability; Probability; Probability Theory& Applications


Introduction
Diffusion and Gauss-Markov processes in the presence of a reflecting boundary are widely used for modeling continuous time phenomena in many scientific fields, such as neurosciences, mathematical biology, finance, and queueing systems. In some instances, these processes are derived as approximations of discrete-state Markovian models which, although more appropriate in describing the behavior of the real system, are often hard to study from both analytical and computational points of view (cf., for instance, Chen & Whitt, 1993;Harrison, 1985;Kushner, 2001).
Typical situations arise in population dynamics with immigration in which the total number of individuals is bound to take non-negative values, so it is necessary to impose a reflection condition at zero state (cf. Renshaw, 2011;Ricciardi, 1986 and references therein).

ABOUT THE AUTHORS
The research interest of the authors includes theory and simulation of stochastic processes with applications to biomathematical modeling and queueing systems. Specifically, they are involved in the development of analytical and computational methods for the solution of several theoretical and algorithmic problems suggested, or imposed, by the considered models. Their research activity has been particularly intense within the following themes: (i) the study of the first-passage time densities for diffusion processes and Gauss-Markov processes and their asymptotic behavior; (ii) the development of mathematical methods and tools mainly of a probabilistic and computational nature to analyze biological dynamics, among them the input-output behavior of neurons subject to random perturbations.

PUBLIC INTEREST STATEMENT
In a wide variety of different domains of applications, including finance, mathematical biology, physics, neurobiology, queueing systems and engineering, particular attention is dedicated to Gauss-Markov processes restricted from below by particular time-dependent reflecting boundaries. We formulate algorithms to simulate the sample paths of such stochastic processes. Furthermore, we construct histograms of first passage times for the restricted process. Particular attention has been dedicated to restricted Wiener and Ornstein-Uhlenbeck processes for their central role in the class of Gauss-Markov processes.
In heavy traffic conditions, the dynamic of a queue can be approximated by a diffusion process with a zero reflecting boundary (cf. Di Crescenzo, Di Crescenzo & Nobile, 1995;Iglehart & Whitt, 1970;Reiman, 1984;Ward & Glynn, 2005). In the financial applications, diffusion models have been used to capture the stochastic movement of the short-term interest rate in the market, for which negative values are not allowed (cf. Cox, Ingersoll, & Ross, 1985;Goldstein & Keirstead, 1997;Han, Hu, & Lee, 2016;Linetsky, 2005).
In Buonocore et al. (2014), we have considered Gauss-Markov processes, restricted from below by particular time-dependent reflecting boundaries, and we have explicitly determined the transition probability density functions (pdf's), the conditional mean, and the second-order conditional moment. Furthermore, in some special cases, the FPT density through a time-dependent threshold has been computed using both asymptotic methods either numerical or simulation techniques. Instead, in Buonocore et al. (2015), the restricted Gauss-Markov processes have been used to construct inhomogeneous leaky integrate-and-fire stochastic models for single neurons activity in the presence of a lower reflecting boundary and periodic input signal.
The present paper deals with various aspects of the simulation of the restricted Gauss-Markov processes considered in Buonocore et al. (2014). The paper is organized as follows. In Section 2, we start considering a Gauss-Markov process {Y(t), t ∈ T} conditioned on the initial state y at time and we give an exact method for its simulation at desired time instants. In Section 3, we take into consideration the restricted stochastic process X(t), obtained considering Y(t) in the presence of a special time-dependent lower reflecting boundary. The reflecting boundary takes the form where the real constants A and B are chosen in order to ensure that y ≥ ( ). Particular attention is dedicated to restricted Wiener and Ornstein-Uhlenbeck processes, in view of their potential interest in a wide variety of different applications domains (Di Crescenzo & Giorno, 2012;Giorno, Nobile, & di Cesare, 2012;Linetsky, 2005;Ricciardi & Sacerdote, 1987;Wonho, 2009). In Section 4, we formulate two algorithms to simulate the sample paths of the process X(t), distinguishing the case in which A = 0 from the case A ≠ 0. When A = 0, it is easy to realize the sample paths of X(t), using a principle similar to that of reflection of Brownian motion. Instead, when A ≠ 0, a different algorithm has been implemented that constitutes an extension to the Gauss-Markov processes of the one given in Asmussen et al. (1995) and Kroese et al. (2011) to simulate the sample paths of a reflected Wiener process with negative drift. Moreover, in Section 5, we use the algorithms proposed in Section 4 to construct the histogram of FPT density for the process X(t). Finally, the R codes to generate sample paths and to simulate the first passage times for Wiener and Ornstein-Uhlenbeck processes are listed in Appendix 1.
Throughout this paper, the symbol d = denotes equality in distribution.

Simulation of Gauss-Markov processes conditioned on the initial state
In this section, we first generate the sample paths of the conditioned Gauss-Markov process Y(t), according to the stochastic equation (1), using an exact simulation method (cf., for instance, Kroese et al., 2011). Let s and t be two time instants, such that t 0 < s < t; by virtue of (1), we obtain the stochastic equation Hence, we can write with s, t ∼  (0, 1). Substituting (9) in (8), for t 0 < s < t, one has: By applying (10) at times s = t k−1 and t = t k for k = 1, 2, … , n, we get where 1 , 2 , … , n is a sequence of independent and identically distributed (iid) standard normal random variables due to the independence of the increments of a Brownian motion for non-overlapping intervals.
In particular, by setting t k = + k h, where h is the time discretization step, for the Wiener process, Equation (11) becomes whereas for the Ornstein-Uhlenbeck process Equation (11) is: The above exact method for simulate Y(t) at the desired times t 1 , t 2 , … , t n can be formulated as follows.
In Buonocore et al. (2014), we have used the Algorithm 2.1 to generate sample paths of the membrane potential in a neuronal model based on a time inhomogeneous Ornstein-Uhlenbeck process.

Restricted Gauss-Markov processes conditioned on the initial state
In this section, we recall some results obtained in Buonocore et al. (2014); they are the starting point for the development of the simulation techniques considered in the next sections.
We denote by {X(t), t ∈ T} the stochastic process with state space [ (t), +∞), obtained by con-

sidering Y(t) in presence of a reflecting boundary
The choice of the real constants A and B depends on the functions m(t), h 1 (t), h 2 (t) and on the initial state y of the process Y(t) at time ; since (t) is a lower reflecting boundary, the real constants A and B are chosen in order to ensure that the starting point y of the process X(t) is greater or equal to the reflecting boundary at the initial time , i.e. y ≥ ( ). As proved in Buonocore et al. (2014), the pdf of the conditioned process {X(t), t ∈ T} is: We recall that the normal pdf f Y (x, t), with mean and variance (2), satisfies the Fokker-Planck equation, with drift and infinitesimal variance respectively, and the delta initial condition (cf. Di Nardo et al., 2001). Hence, for x > (t) and y > ( ), the pdf f X (x, t), given in (15), satisfies the Fokker-Planck equation and the associated initial condition (see, Buonocore et al., 2014) with the additional requirement that a reflecting condition is imposed on the boundary (t), i.e.
where the drift A 1 (x, t) and the infinitesimal variance A 2 (t) are given in (16).
Since (t) is a reflecting boundary, the total probability mass of We remark that the boundary (14) includes several types of reflecting boundaries. In particular, for the Wiener process, from (4) and (14) one has i.e. the reflecting boundary is linear. When A = 0, the angular coefficient of the reflection line coincides with the drift of the Wiener process. Instead, for the Ornstein-Uhlenbeck process, from (6) and (14) one obtains i.e. the reflecting boundary is an hyperbolic function.

Simulation of the restricted Gauss-Markov process
In this section, we formulate two different algorithms to simulate the sample paths of the process X(t) originated from y at time , in the presence of the reflecting boundary (t), given in (14), with y ≥ ( ).
If A = 0, the problem is simple because the principle of reflection can be used (see Borodin & Salminen, 2002 for Brownian motion).

Proposition 4.1 Let
For t ≥ and y ≥ ( ), the stochastic process {X s (t), t ∈ T}, with is characterized by pdf: is the Gauss-Markov process conditioned to start to y at time and restricted to [ (t), +∞), with the reflecting boundary (t) as in (19).
Proof From (1) and (20) we first note that so that, recalling that y ≥ ( ), it follows X s ( ) = y. For x ≥ (t) and y ≥ ( ), one has We now assume that h 2 (t) > 0, so that recalling (1) and (19), from (22) one has: We note that lim x→+∞ P{X s (t) < x} = 1. Furthermore, recalling that for a < b, from (23) and (24), one has: where M(t|y, ) and V(t| ) are given in (2). Hence, (21) follows immediately from (25). Proceeding in a similar way for h 2 (t) < 0, relation (25) still holds. Since the right-hand side of (21) is equal to that of (15) with A = 0, then X s (t) d = X(t) for all t ∈ T. ✷ Proposition 4.1 shows that, for all t ≥ , X s (t) is identically distributed as X(t) originated at X( ) = y ≥ ( ) in the presence of the reflecting boundary (19). Hence, under the assumptions of Proposition 4.1, a method to simulate X(t) at the desired times t 1 , t 2 , … , t n can be expressed as follows.
In Buonocore et al. (2014), by means of the Algorithm 4.1, we have generated sample paths of the membrane potential in a neuronal model based on a special time inhomogeneous Ornstein-Uhlenbeck process restricted by a lower reflecting boundary.
In Figure 1 The Wiener process Y(t) and the related restricted process X(t) do not admit a steady-state density, as highlighted by the behavior of the trajectory in Figure 1. Instead, as it is also evident in Figure 2(a), when < 0, the Ornstein-Uhlenbeck process Y(t) and the related restricted process X(t) admit the following steady-state densities: Indeed, when < 0, by virtue of (6), from (2) ) is a normal density with mean ∕| | and variance 2 ∕(2| |). Furthermore, recalling that A = 0 and < 0, Note that for the Ornstein-Uhlenbeck process, restricted or not, with > 0, no steadystate densities exist (see Figure 2(b)).
If A ≠ 0, the problem is more complex since one cannot use, as in Proposition 4.1, the principle of reflection.
Proposition 4.2 We consider the stochastic process {X s (t), t ∈ T}, such that where In (27) and (28), we assume that = 1, Hence, for all t ∈ T Proof We assume that h 2 (t) > 0 for all t ∈ T, so that = 1,Ã = −A,B = B. Then, making use of (27) and (28), for x ≥ (t) and y ≥ ( ) one has:
where Z(t) denotes a Wiener process with drift −A and unitary infinitesimal variance, conditioned on z = {y − m[r −1 ( )] − Ah 1 [r −1 ( )]}∕h 2 [r −1 ( )] at time . We now consider the process conditioned on Ẑ ( ) = z ≥ B (cf., for instance, Asmussen et al., 1995). For all t ≥ , the distribution function of Ẑ (t) is: We note that where (a, b, t) is the joint density of Z(t) and inf ≤s≤t [Z(s)] (cf., for instance, Borodin & Salminen, 2002): From (31) and (32), one can derive the pdf of Ẑ (t): (33) The right-hand side of (33) identifies the transition density of a Wiener process Ŵ (t) with drift −A, unitary infinitesimal variance, restricted to the interval [B, +∞) with a reflecting boundary in B, conditioned on z ≥ B at time (cf., for instance, Cox & Miller, 1970). The distribution function of Ŵ (t) is Now, recalling (34), Equation (29) becomes: for x ≥ (t) and y ≥ ( ), with M(t|y, ) and V(t| ) given in (2). The right-hand side of (35) identifies with the transition distribution function of the conditioned Gauss-Markov process X(t), restricted to [ (t), +∞) in the presence of the reflecting boundary (t) as in (14). Hence, for h 2 (t) > 0, X s (t) d = X(t) for all t ∈ T. This completes the proof in the case h 2 (t) > 0. The proof when h 2 (t) < 0 for all t ∈ T follows in a similar way. ✷ Proposition 4.2 allows to construct an exact algorithm for generating X s (t) at the specified time instants. To this end, we assume that h 2 (t) > 0 for all t ∈ T. We set t = t k = + k h (k = 0, 1, … , n; h > 0), where h denotes the time discretization step, and let Since r(t) is a non-negative and monotonically increasing function, we note that the following relation holds: By setting t = t k = + k h in (27), one has: so that, making use of (37) in (38), one obtains: where for k = 1, 2, … , n it results: By virtue of (38), M k−1 can be expressed in terms of X s (t k−1 ) as: [−Z(s)] (k = 0, 1, 2, … , n).
for k = 1, 2, … , n. Substituting (28) with t = t k−1 , = 1 and Ã = −A in (41), the first term of the maximum in (39) becomes: Now we want to simulate the random variable S k that appears in (40). Let Z 1 (t) = −Z(t) for all t ∈ T, with Z(t) given in (28). Hence, We now use the well-known result for the supremum of the Wiener process (cf. Borodin & Salminen, 2002) with drift , infinitesimal variance 2 , starting from z 0 at time 0 : i.e. an exponential distribution that does not depend on . Making use of (44) with = A, 2 = 1, for x ≥ (t k ) and b ≥ z 0 . Let S k be the random variable obtained by conditioning S k with the requirements that Z 1 [r(t k−1 )] = z 0 and Z 1 [r(t k )] = z. Hence, applying the inverse transform algorithm, we set with U k ∼  (0, 1), and solve for S k in terms of random variable U k . Hence, Since we have imposed that Z 1 [r(t k−1 )] = z 0 and Z 1 [r(t k )] = z, by virtue of (28), the value −(z − z 0 ) in (46) can be obtained by simulation of the following random variable where the last identity follows from (9), with k ∼  (0, 1). Hence, making use of (46) and (47), the random variable S k in (40) Since h 2 (t) > 0 for all t ∈ T, it ensures that S k ≥ (t k ) is necessary to take the positive square root in (48). The case h 2 (t) < 0 for all t ∈ T follows in a similar way. Then, the following algorithm can be implemented: where with U 1 , U 2 , … , U n iid =  (0, 1) and 1 , 2 , … , n iid =  (0, 1). In (49) and (50), we assume that = 1,Ã = −A when h 2 (t) > 0 for all t ∈ T, whereas = −1,Ã = A when h 2 (t) < 0 for all t ∈ T.
Note that, by virtue of (9), the first term of the maximum in (49) follows from (42). Furthermore, the second term of the maximum in (49) follows from (48).
We remark that the Algorithm 4.2 is a generalization to restricted Gauss-Markov processes of the simulation algorithm given in Asmussen et al. (1995) and Kroese et al. (2011) for a Brownian motion with a negative drift and unit infinitesimal variance in the presence of a reflecting boundary in zero state.
In Figure 3, we plot a simulated sample path of the Wiener process Y(t) obtained with Algorithm 2.1, by choosing = −1.0, 2 = 1.0, y = 3, = 0 in (5), and a simulated sample path of the related restricted process X(t) in [ (t), +∞), with (t) given in in (17). In Figure 3(a), we set A = 0.5 and B = 1.0, so that from (17) it follows that (t) = −0.5 t + 1.0, whereas in Figure 3(b) we set A = 1.5 and B = 1.0, so that from (17) one has (t) = 0.5 t + 1.0. The sample path of X(t) is obtained with the Algorithm 4.2. The time discretization step is h = 10 −3 . The R code of Figure 3(a) is given in Appendix 1.
Furthermore, in Figure 4, we plot a simulated sample path of the Ornstein-Uhlenbeck process Y(t) achieved with Algorithm by choosing = −0.5, = 1.0, 2 = 1.0, y = 3, = 0 in (7), and a simulated sample path of the related restricted process X(t) in [ (t), +∞), with the reflecting boundary (t) given in (18). In Figure 4(a), we set A = 0.35 and B = 1.35, so that from (18) it follows that (t) = 2 − e −t∕2 + 0.35 e t∕2 ; instead, in Figure 4(b) we set A = −0.01 and B = 2.19, so that from (18) one has (t) = 2 + 0.2 e −t∕2 − 0.01 e t∕2 . The sample path of X(t) is obtained via the Algorithm 4.2. The time discretization step is h = 10 −3 . The R code of Figure 4(a) is given in Appendix 1. In Figures 5 and 6, we show that, starting from the same realization of the Wiener and Ornstein-Uhlenbeck processes Y(t) (green), respectively, Algorithms 4.1 and 4.2 provide for A = 0 and B = 0 different trajectories for the processes X(t) (red). However, as proved in Propositions 4.1 and 4.2, the simulated restricted processes are characterized by the same pdf. We finally note that for A = 0, the Algorithm 4.2 has an higher computational cost than that of the Algorithm 4.1. Specifically, for the comparison of different numerical or simulation methods, the analysis of their accuracy as function   of the discretization step, as well as the study of their numerical complexity, we refer to the books (Asmussen & Glynn, 2007;Kloeden & Platen, 1999;Kroese et al., 2011).

First passage time simulation of restricted Gauss-Markov process
In this section, we consider the first passage time problem for the Gauss-Markov process and for the related restricted process.

FPT simulation of the conditioned Gauss-Markov processes
For the process Y(t), the FPT problem through a boundary S(t) ∈ C 1 (T) can be considered. Let be the random variable that denotes the FPT of the Gauss-Markov process Y(t) from Y( ) = y to the continuous boundary S(t). Let denotes the FPT pdf. In order to simulate g Y [S(t), t], we generate sample paths of the Gauss-Markov process Y(t) using the Algorithm 2.1. Then, the method for the simulation of FPT for the process Y(t) through S(t) can be made as follows: By implementing the previous procedure for N times, one obtains a collection of N simulated first passage times of Y(t) through S(t). Hence, an estimation of the FPT pdf can be achieved by the histogram of such first passage times.

FPT simulation of the restricted Gauss-Markov process
For the process X(t), restricted to [ (t), +∞) with the reflecting boundary (t) given in (14), the FPT problem through a boundary S(t) > (t) can be considered. To this purpose, let S(t) be a C 1 (T)-class function, such that S(t) > (t) for all t ∈ T. For y ≥ ( ), we denote by the random variable FPT of X(t) from X( ) = y to the boundary S(t). Further, let be the FPT pdf.
As proved in Buonocore et al. (2014) , and S(t) are not depend on time t, then the FPT moments through a constant boundary for the process X(t) can be obtained making use of the Siegert formula (cf. Ricciardi, Di Crescenzo, Giorno, & Nobile, 1999;Siegert, 1951). Indeed, with t 0 (S|y) = 1, and where Explicit expressions of moments for the Ornstein-Uhlenbeck process with a constant reflecting boundary are given in Inoue, Sato, and Ricciardi (1997).
The sample paths of the restricted process X(t) are constructed taking into account the presence of the reflecting boundary (t) given in (14), using one of the algorithms described in Section 4. Then, the method for the simulation of FPT for the restricted process X(t) through S(t) can be carried out as follows: By implementing the previous procedure for N times, one obtains a collection of N simulated first passage times X(t) through S(t). Hence, an estimation of the FPT pdf can be achieved by the histogram of such first passage times.
For the restricted Wiener process X(t), with = 0, 2 = 1.0, y = 3.0, = 0, in Figure 7 we plot the histograms of a collection of N = 5000 simulated FPT through the boundary S = 4.5, obtained making use of Algorithm 4.1 on the left and of Algorithm 4.2 on the right; in this case, we set A = 0 and B = 0 in (17). In Figure 7, the histograms of the FPT densities are compared with the exact FPT density (red curve), obtained via (55). It appears from Figure 7 that the FPT histograms fit the exact FPT density reasonably well. Furthermore, in Table 1, we compare the FPT mean t 1 (S|y), the FPT variance V(S|y) = t 2 (S|y) − [t 1 (S|y)] 2 , and the FPT coefficient of variation C(S�y) = √ V(S�y)∕t 1 (S�y), derived from (56) For the restricted Ornstein-Uhlenbeck process X(t), with = −0.1, = 0.0, 2 = 1.0, y = 3.0, = 0, in Figure 8 we plot the histograms of the FPT pdf through the boundary S = 4.5 obtained making use of Algorithm 4.1 on the left and of Algorithm 4.2 on the right; furthermore, note that A = B = 0 in (18). We have considered a collection of N = 5000 simulated FPT. In this case, there is no closed-form expression for the FPT density and the simulation tool is necessary. Furthermore, the FPT mean t 1 (S|y), the FPT variance V(S|y), and the FPT coefficient of variation C(S|y) can be obtained from (56). Hence, in Table 2 are reported t 1 (S|y), V(S|y), and C(S|y) and the related sample values, obtained from the collection of N = 5000 simulated FPT of Figure 8.

Concluding remarks
In conclusion, we propose algorithms for the simulation of the Gauss-Markov processes conditioned to start to y at time and restricted to [ (t), +∞] with the time-dependent lower reflecting boundary (t) = m(t) + A h 1 (t) + B h 2 (t), where the real constants A and B are chosen in order to ensure that y ≥ ( ). Specifically, we formulate two algorithms to simulate the sample paths of the process X(t) for both cases A = 0 and A ≠ 0. When A = 0, the sample paths of the restricted Gauss-Markov processes are obtained using a reflection principle type (see Algorithm 4.1), whereas if A ≠ 0, the sample paths are constructed via Algorithm 4.2. Finally, we use the algorithms proposed to construct the histogram of FPT density for the restricted process. Particular attention has been dedicated to restricted Wiener and Ornstein-Uhlenbeck processes for their central role in the class of Gauss-Markov processes.