On the Simulation of a Special Class of Time-Inhomogeneous Diffusion Processes

: General methods to simulate probability density functions and ﬁrst passage time densities are provided for time-inhomogeneous stochastic diffusion processes obtained via a composition of two Gauss–Markov processes conditioned on the same initial state. Many diffusion processes with time-dependent inﬁnitesimal drift and inﬁnitesimal variance are included in the considered class. For these processes, the transition probability density function is explicitly determined. Moreover, simulation procedures are applied to the diffusion processes obtained starting from Wiener and Ornstein–Uhlenbeck processes. Speciﬁc examples in which the inﬁnitesimal moments include periodic functions are discussed.


Introduction
Time-inhomogeneous stochastic diffusion processes are often used to model real dynamic systems in various scientific areas as neurosciences, population dynamics, queueing systems, finance, and economics (cf., for instance, Buonocore et al. [1,2], Albano and Giorno [3], Giorno and Spina [4], Ricciardi et al. [5], Renshaw [6]), Di Crescenzo et al. [7], Linetsky [8], Glasserman [9]). Specifically, in Buonocore et al. [1,2], restricted Gauss-Markov processes are used to construct inhomogeneous leaky integrate-and-fire stochastic models for single neurons activity in the presence of a reversal hyperpolarization potential and time-varying input signals. In Albano and Giorno [3], a time-inhomogeneous Ornstein-Uhlenbeck is considered as a model for the membrane potential activity of a single neuron and a statistical procedure to fit the constant parameters and the time-dependent functions is proposed. Moreover, in Giorno and Spina [4], a time-inhomogeneous Ornstein-Uhlenbeck diffusion process with jumps is also analyzed. In Ricciardi et al. [5], special emphasis is put on neuronal firing problems and on the description of population dynamics, for which the first-passage time distribution and its statistics carry a fundamental relevance. In Renshaw [6], many aspects of population dynamics are covered, including Wiener and Ornstein-Uhlenbeck diffusion models and simulation techniques. In Di Crescenzo et al. [7], the authors derive a heavy-traffic approximation that allows for approximating the state of the systems by a time-non-homogeneous Wiener process subject to jumps due to catastrophes and random returns to the zero state due to repairs. In Linetsky [8], explicit analytical expressions for transition densities of Brownian motion with drift, the Ornstein-Uhlenbeck process, and affine (square-root) diffusion with one or two reflecting barriers are obtained in the context of economics and finance. In Glasserman [9], the Monte Carlo simulation plays an essential tool in the pricing of derivative securities and in risk management.
To include several aspects of the real phenomena, one is led to consider increasing complex processes (cf., for instance, Di Crescenzo et al. [10], Giorno et al. [11,12], Abundo [13], Veestraeten [14], Molini et al. [15], Lim and Muniandy [16], Jeon et al. [17]). In these contexts, the focus is on specific probabilistic characteristics that define the behavior of the stochastic processes as the transition distributions and the related moments as well as the first passage time through specific time-dependent boundaries. In particular, in Di Crescenzo et al. [10], and in Giorno et al. [11,12], new procedures for constructing transition probability density functions and first passage time densities through constant boundaries are proposed for generally time-inhomogeneous diffusion processes. In Abundo [13], some asymptotic results for diffusions and Gauss-Markov process and their fractional integrals are obtained. In Veestraeten [14], the transition and first hitting time densities and moments for the Ornstein-Uhlenbeck process between exponential thresholds are derived. Moreover, in Molini et al. [15], first passage time statistics, such as the survival probabilities and first passage time densities, are obtained analytically for the Brownian motion driven by timedependent drift and diffusion coefficients. In Lim and Muniandy [16] and Jeon et al. [17], some Gaussian models for anomalous diffusion are considered, and the first passage time problem is discussed. These processes are used to model physical and biological systems.
Moreover, various efforts are direct to computational approaches (Di Nardo et al. [18], Taillefumier and Magnasco [19], D'Onofrio and Pirozzi [20]) that also involve methods of statistical inferences (Albano et al. [21], Albano and Giorno [22], Ramos-Ábalos et al. [23]). Specifically, in Di Nardo et al. [18] and Taillefumier and Magnasco [19] methods to construct first-passage-time probability density functions for Gauss-Markov processes through timedependent boundaries are proposed. In D'Onofrio and Pirozzi [20], the problem of escape times from a region confined by two time-dependent boundaries is considered for a class of Gauss-Markov processes. Moreover, numerical procedures to infer the models based on time-inhomogeneous diffusion processes are proposed in Albano et al. [21], Albano and Giorno [22], and Ramos-Ábalos et al. [23].
However, sometimes these strategies are not actionable or not convenient and an approach based on the simulation of the processes is required (see, for instance, Buonocore et al. [24,25], Tuerlinckx et al. [26], Di Crescenzo et al. [27], Giraudo et al. [28], Herrann and Zucca [29], Headrick and Mugdadi [30], Devroye [31], Iacus [32]). In particular, in Buonocore et al. [24,25], algorithms for the simulation of sample paths and to generate random variates from probability density function of Gauss-Markov processes, restricted by particular time-dependent reflecting boundaries, are proposed. Furthermore, in Tuerlinckx et al. [26], some methods for the simulation of the Wiener process with constant drift and variance are described and compared on two criteria: simulation speed and accuracy of the simulation. In Di Crescenzo et al. [27], Giraudo et al. [28] and Herrann and Zucca [29], simulation procedures to estimate first-passage-time densities are constructed. Moreover, in Headrick and Mugdadi [30], algorithms for extending the class of generalized lambda distributions from univariate to multivariate data generation are presented. Finally, Iacus [32] focuses on the simulation and inference for general stochastic differential equations.
Random variates generation is a fundamental aspect of simulation modeling and analysis (cf., for instance, Devroye [31]). Indeed, in real applications, the distribution of a function of more iid (independent and identically distributed) random variables (statistics, estimators) is often computed via the simulation only, when the resulting distribution can not be deduced analytically. For this purpose, it is necessary to identify suitable methods to simulate the single random variable of interest. Standard procedures, as the inverse transformation method and the acceptance-rejection technique, often cannot be applied successfully so that ad hoc methods are needed to generate a sample of observations and then plot the resulting empirical histogram that can be compared with the related known probability density. The knowledge of effective simulation methods for single random variates generation allows for obtaining the histogram of the distribution of the statistics and the estimators that often have analytically non-computable distributions.
In Giorno and Nobile [33], a special class of time-inhomogeneous diffusion processes is defined and analyzed. These processes are obtained by using the composition of two Gauss-Markov processes conditioned to start from the same initial state. They are useful to model dynamical systems that switch randomly between two different regimes. This class includes diffusion processes with time-dependent infinitesimal moments for which the transition probability density function (pdf) is a mixture of two normal densities with weights depending on the initial condition and on a real parameter ϑ ∈ [0, 1]. General methods to analyze the first-passage time (FPT) are given, and FPT densities are explicitly obtained for suitable time-varying boundaries.
In the present paper, we propose theoretical and computational approaches based on the simulation to investigate on the transition densities, on the related conditional moments and on the FPT densities for the diffusion processes considered in [33].
The paper is organized as follows. In Section 2, starting from two Gauss Markov processes conditioned on the same initial state, we construct a continuous process Y(t) by using the composition method. We formulate an algorithm to obtain the simulated density and a procedure for the simulation of the sample paths of Y(t). In Section 3, we focus on a class of time-inhomogeneous diffusion processes Z(t), whose transition pdf identifies with the density of Y(t). In particular, for ϑ = 0 and ϑ = 1, one has Z(t) = Y(t) for all t, so that an algorithm to simulate the FPT through a general time-dependent boundary is formulated. Furthermore, we show that such algorithm can be generalized to the case 0 < ϑ < 1. In Sections 4 and 5, we apply the theoretical results and the proposed algorithms to diffusion processes obtained by the composition of Wiener processes and Ornstein-Uhlenbeck processes, respectively, by using the statistical environment and the language R. In both these cases, for fixed time instants, we simulate the random variable describing Y(t), and we show that the transition pdf of the related diffusion process Z(t) can be superimposed over the histogram of the process Y(t) obtained via the simulation method. Moreover, for ϑ = 0 and ϑ = 1, making use of the simulation of the sample-paths, we obtain the histogram of first passage times of Z(t) through time-dependent boundaries. Finally, the histogram of first passage times is compared with the closed form FPT density through special boundaries.

Composition Method for Gauss-Markov Processes
Let m(t), h 1 (t), h 2 (t) be C 1 (T)-class functions, where C 1 (T) denotes the set of continuously differentiable functions on T, with a T continuous parameter set, such that h 2 (t) = 0 and r(t) = h 1 (t)/h 2 (t) is a non-negative and monotonically increasing function. Let {X i (t), t ≥ τ} be the Gauss-Markov processes conditioned to start from y at time τ (see, Mehr and McFadden [34]): with mean and variance Starting from the Gauss-Markov processes X 1 (t) and X 2 (t), conditioned to start from y at time τ, we use the composition method (cf., for instance, Ross [35]) to construct a new stochastic process {Y(t), t ≥ τ}.
Let 0 ≤ ϑ ≤ 1 be a real number and let X 1 (t) and X 2 (t) be the Gauss-Markov processes conditioned to start from y at time τ, defined in (1). For any fixed t ∈ T, we assume that U y,τ (t) is a random variable uniform in (0, 1), independent of X 1 (t) and X 2 (t). Then, the stochastic process {Y(t), t ≥ τ}, defined as with where The density f Y (x, t) in (6) is a mixture, or a composition, of the two normal densities f X 1 (x, t) and f X 2 (x, t), with means and variances given in (3). In particular, by virtue of (4) and (6), we note that In the following, we formulate an algorithm to generate random variates from the pdf (6) by using the stochastic equations (1) and the composition method (4). Furthermore, we derive an algorithm for the simulation of the sample paths of Y(t) via the composition method.

Simulated Pdf of the Process Y(t)
Let t ∈ T be a fixed instant and (y, τ) ∈ R × T fixed real constants. For t ≥ τ, we obtain a random sample of N observations of Y(t), and we construct the histogram of the random sample with the density (6) as a function of x.

Algorithm 1
Let (y, τ) ∈ R × T be fixed real constants and t ∈ T be a fixed instant such that t ≥ τ. STEP 1: Generate ξ τ,t ∼ N (0, 1) and U y,τ (t) ∼ U (0, 1), with U y,τ (t) and ξ τ,t independent random numbers; STEP 2: From (1), generate X i (t) as STEP 3: Making use of (7), generate Y(t) via (4); STEP 4: Repeat Steps 1 and 2 for N times, obtaining a random sample of size N from the density (6). The simulated pdf is then depicted by means of a histogram as function of x.
The related sample moments can also be obtained.

Simulation of the Sample Paths of Y(t)
We first generate the sample paths of the Gauss-Markov processes X 1 (t) and X 2 (t), according to the stochastic equations (1) by using an exact simulation method (cf., for instance, Kroese et al. [36]).

Some Time-Inhomogeneous Diffusion Processes
Let {Z(t), t ∈ T} be a time-inhomogeneous diffusion process with infinitesimal drift and infinitesimal variance where 0 ≤ ϑ ≤ 1 is a real number and m(t), h 1 (t), h 2 (t), v(t) are chosen as in Section 2. As proved in Giorno and Nobile [33], if one of the following assumptions is satisfied: The density (6) identifies with the transition pdf f (x, t|y, τ) = ∂P{Z(t) ≤ x|Z(τ) = y}/∂x of the diffusion process Z(t). Then, in the cases (a) and (b), for all fixed t ≥ τ, one

Case (b)
If m(t) = c h 1 (t) for t ∈ T and c ∈ R, the infinitesimal moments (14) of Z(t) become: and, from (6), one obtains the transition pdf: where M 1 (t|y, τ) and V(t|τ) are given in (3), with m(t) = c h 1 (t). Then, the conditional mean and the variance of Z(t) are: Furthermore, for the time-inhomogeneous diffusion process defined in (14), let be the FPT of Z(t) from Z(τ) = y to the boundary S(t) ∈ C 1 (T) and let g Z [S(t), t|y, τ] = ∂P(T Z ≤ t)/∂t be the FPT pdf, with y = S(τ). The FPT pdf g Z [S(t), t|y, τ] is a solution of the first-kind Volterra integral equation: Making use of (18) in (21), one obtains: Similarly, the FPT densities of the processes X 1 (t) and X 2 (t) are also solutions of first-kind Volterra integral equations: In the case (b), the transition densities of the processes X 1 (t) and X 2 (t) satisfy the following symmetry relation: so that, from (23), one has: Hence, for the diffusion process Z(t) with infinitesimal moments (17), from (22) and (24), one obtains where 0 ≤ ϑ ≤ 1, . Equation (25) shows that g Z [S(t), t|y, τ] is a mixture of two FPT densities g X 1 [S(t), t|y, τ] and g X 2 [S(t), t|y, τ]. Therefore, the random variable T Z is identically distributed as where T X i is the FPT of X i (t) through S(t) for i = 1, 2 and U y,τ ∼ U (0, 1) is independent of T X 1 and T X 2 . We note that, for the diffusion process Z(t) with infinitesimal moments (17) if then the FPT density can be expressed in closed form as (see Giorno and Nobile [33]): with f Z (x, t|y, τ) given in (18). Moreover, the first passage of Z(t) through S(t), given in (27), is a certain event if and only if [ Making use of the Algorithm 1, it is possible to obtain a random sample of N observations of Y(t) and to construct the histogram of the random sample. Such histogram can be compared with the transition pdf (15), in the case (a), and with the transition pdf (18) Applying the Algorithm 2, we can produce a sample path for Y(t) at the desired times t 1 , t 2 , ..., t n via the composition method.
For ϑ = 0 or ϑ = 1, recalling that Z(t) = Y(t), we can obtain the following method for the simulation of first passage times for the process Z(t) through S(t) if y < S(τ):

Algorithm 3
Let t 1 < t 2 < . . . < t n be the set of distinct time instants for which the simulation of the process Z(t) is desired.
with ξ 1 , ξ 2 , . . . , ξ n independent normal random numbers. STEP 3: If Z(t k ) ≥ S(t k ), then collect the first passage time t k and stop, else k ← k + 1 and go to Step 2.
The implementation of the previous procedure for N times allows for obtaining a collection of N simulated first passage times of Z(t) through S(t). Then, the histogram of such first passage times can be used to obtain an estimation of the FPT pdf. When y > S(τ), in STEP 3 of Algorithm 3, it is necessary to change Z(t k ) ≤ S(t k ).
In the case (b), when 0 < ϑ < 1, implementing the Algorithm 3 N times for ϑ = 0 and for ϑ = 1, we obtain two collections of N simulated first passage times t 2,1 , . . . , t 2,N for X 2 (t) and t 1,1 , . . . , t 1,N for X 1 (t). Recalling (26), one obtains a collection of N simulated first passage times of Z(t) through S(t) as follows: where u 1 , u 2 , . . . , u N are independent uniform numbers in (0, 1). Hence, an estimation of the FPT pdf of Z(t) through S(t) can be achieved by the histogram of the first passage times t 1 , t 2 , . . . , t N . For the diffusion process Z(t), the estimation of the FPT pdf via Algorithm 3 and its generalization to the case 0 < ϑ < 1 depends on the infinitesimal moments; on the discretization step, on the time-varying boundary and on the choice of the initial state with respect to boundary. In some cases, the FPT densities can present heavy tails as the time increases, so that the simple path of the stochastic process can take a long time to reach and cross the boundary. The simulations run until STEP 3 in Algorithm 3 are satisfied, but, in some cases, can be necessary to set a maximum simulation time in order to avoid very long and time-expensive simulations. Of course, in this case, it is necessary to count the realizations that do not reach the threshold and give an estimation of the success probability in such a way as not to affect the estimation of the FPT density.
In the sequel, particular attention is dedicated to the simulation of the processes generated via the Wiener and the Ornstein-Uhlenbeck processes, with continuous parameter set T = [0, +∞).
We assume that β(t) = 0 for t ≥ 0. By choosing (y, τ) = (0, 0), from (31), one has m(0) = h 1 (0) = 0, so that the assumptions (a) are satisfied; hence, from (15), one obtains: and, making use of (31) in (16), one has: Furthermore, the assumption (b) holds if and only if In this case, from (32), the infinitesimal moments of Z(t) are: and, from (18), it follows: Making use of (31) in (19), one has: Moreover, due to (27) and (28), for the process Z(t), the FPT pdf through the boundary with f Z (x, t|y, τ) given in (37). Finally, the first passage of Z(t) through S(t), given in (39) Example 1. We consider the diffusion process Z(t) in (32) with so that Then, from (32), it follows that the infinitesimal drift and the infinitesimal variance of the process Z(t) are with 0 ≤ ϑ ≤ 1. We suppose that (y, τ) = (0, 0), so that the assumptions (a) are satisfied. In this case, from (7), we obtain: with ξ t ∼ N (0, 1); hence, from (4) and (5), one has and Z(t) d = Y(t) for all fixed t ≥ 0. Algorithm 1 can be used to compare the transition pdf (33) of Z(t), being β(t) and σ 2 (t) given in (41), with the histograms of the random sample of observations of Y(t), obtained via (45), for different instants t. Specifically, in Figure 1, the histogram of a random sample of N = 10 6 observations of Y(t) is obtained for different times t, by choosing β(t) and σ 2 (t) as in (41), with ϑ = 0.1, γ = 1.5, σ 2 = 2, ν = 0.1, Q = 1. refer to the simulation sample mean y t , variance s 2 t and coefficient of variation c t , computed from the same random sample of N = 10 6 observations of Y(t) used in Figure 1. The results of Figure 1 and Table 1 show the good agreement between the probabilistic results and those obtained via simulation. Table 1. The conditional mean, variance, and coefficient of variation of Z(t) are compared with the estimated values y t , s 2 t and c t obtained by means of Algorithm 1 with the same choices of the parameters of Figure 1.  Algorithm 2 permits generating sample paths of the process Y(t). Indeed, for the considered case, from (12), setting t 1 = 0, Y(t 1 ) = 0, we have where ξ 1 , ξ 2 , . . . , ξ n iid = N (0, 1) and with m(t) given in (42). In Figure 2, a sample path of Y(t) is plotted as function of t for ϑ = 0.1 and ϑ = 0.9 by choosing β(t) and σ 2 (t) as in (41), with γ = 1.5, σ 2 = 2, ν = 0.1, Q = 1. When ϑ = 0 or ϑ = 1, Y(t) identifies with a time-inhomogeneous Wiener process Z(t).

Example 2.
We consider the diffusion process Z(t) in (32) with where γ > 2π |ν|/Q, to ensure that B 2 (t) > 0 for t ≥ 0, so that, from (31), one has: Then, from (36), for t ≥ 0 and x ∈ R, one has: The assumption (b) is satisfies with c = 1. In this case, from (7), we obtain: with ξ τ,t ∼ N (0, 1); hence, (4) holds with for all fixed t ≥ 0, we use the Algorithm 1 to compare the transition pdf (37) of Z(t) with the histograms of the random sample of observations of Y(t) for different instants t. Specifically, in Figure 4, the histogram of a random sample of N = 10 6 observations of Y(t) is obtained for different times t, by choosing (y, τ) = (0, 0), β(t) and σ 2 (t) as in (49), with ϑ = 0.4, γ = 1.5, ν = 0.1, Q = 2.  For the same values of the parameters, in Table 2 for t = 1, 5, 10, 15, the conditional mean, variance, and coefficient of variation of Z(t), evaluated via (38), are compared with the simulation sample mean y t , variance s 2 t and coefficient of variation c t , obtained making use of the same random sample of N = 10 6 observations of Y(t) of Figure 4. Table 2. For the same choices of the parameters of Figure 4, the conditional mean, variance, and coefficient of variation of Z(t) are compared with the estimated values y t , s 2 t and c t , obtained by means of Algorithm 1.  Furthermore, due to (39) and (40), for the diffusion process Z(t) having infinitesimal moments (51), the FPT pdf through the boundary is The first passage of Z(t) through (54) is a certain event if and only if [y < S(τ), a 1 ≤ −1] or [y > S(τ), a 1 ≥ 1].
For the same values of the parameters, in Table 3 for t = 1, 5, 10, 15, the conditional mean, variance, and coefficient of variation of Z(t), evaluated via (61), are compared with the simulation sample mean y t , variance s 2 t and coefficient of variation c t , obtained making use of the same random sample of N = 10 6 observations of Y(t) of Figure 6.   Figure 6, the conditional mean, variance, and coefficient of variation of Z(t) are compared with the estimated values y t , s 2 t and c t , obtained by means of Algorithm 1.  Applying the Algorithm 2, we generate sample paths of the process Y(t). For the considered case, from (12), by setting t 1 = 0, Y(t 1 ) = 0, we have: where ξ 1 , ξ 2 , . . . , ξ n iid = N (0, 1) and with m(t) given in (70). In Figure 7, a sample path of Y(t) is plotted as function of t for ϑ = 0.1 and ϑ = 0.9 by using (75) with α = −0.2, β = 1.5 and σ 2 = 2. When ϑ = 0 or ϑ = 1, Y(t) identifies with a time-homogeneous Ornstein-Uhlenbeck process Z(t). In particular, for ϑ = 0, Z(t) has infinitesimal moments B 1 (x) = α x − β and B 2 = σ 2 , whereas, when ϑ = 1, one has B 1 (x) = α x + β and B 2 = σ 2 . We recall that, for the timehomogeneous Ornstein-Uhlenbeck process Z(t) with infinitesimal moments B 1 (x) = α x + β and B 2 = σ 2 , the FTP pdf through the constant boundary S = −β/α is known in closed form: In Figure 8, for ϑ = 1, the histograms of FPT for Z(t), with β = 1.5, σ 2 = 2, α = −0.1 on the left and α = −0.2 on the right, from Z(0) = 0 through the constant boundary S = 15 (on the left) and S = 7.5 (on the right) are plotted as a function of t by using a collection of N = 10 4 simulated first passage times obtained via Algorithm 3. The solid curves in Figure 8 refer to the closed expression of FPT densities for the Ornstein-Uhlenbeck process from Z(0) = 0 through S = −β/α, given in (77).

Conclusions
Starting from two Gauss-Markov processes conditioned on the same initial state, we have constructed a continuous process Y(t) via the composition method. We focused on a class of time-inhomogeneous diffusion processes Z(t), whose transition pdf identifies with the density of Y(t). We have applied the theoretical results and the proposed simulation algorithms to diffusion processes obtained by the composition of Wiener processes and Ornstein-Uhlenbeck processes. For fixed time instants, we have simulated the random variable describing Y(t), and we have shown that the transition pdf of the diffusion process Z(t) can be superimposed over the histogram of the process Y(t), obtained by using the simulation method. Finally, the histogram of first passage times is compared with the closed form FPT density through special boundaries.