Next Article in Journal
On Averaging Principle for Caputo–Hadamard Fractional Stochastic Differential Pantograph Equation
Next Article in Special Issue
On the Estimation of the Persistence Exponent for a Fractionally Integrated Brownian Motion by Numerical Simulations
Previous Article in Journal
Fractional Transformation-Based Intelligent H-Infinity Controller of a Direct Current Servo Motor
Previous Article in Special Issue
On the Absorbing Problems for Wiener, Ornstein–Uhlenbeck, and Feller Diffusion Processes: Similarities and Differences
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Approximating the First Passage Time Density of Diffusion Processes with State-Dependent Jumps

by
Giuseppe D’Onofrio
1,*,† and
Alessandro Lanteri
2,†
1
Dipartimento di Scienze Matematiche, Politecnico di Torino, 10129 Turin, Italy
2
Department of Economics, Management and Quantitative Methods, Università degli Studi di Milano, 20122 Milan, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Fractal Fract. 2023, 7(1), 30; https://doi.org/10.3390/fractalfract7010030
Submission received: 28 November 2022 / Revised: 20 December 2022 / Accepted: 23 December 2022 / Published: 28 December 2022
(This article belongs to the Special Issue Stochastic Modeling in Biological System)

Abstract

:
We study the problem of the first passage time through a constant boundary for a jump diffusion process whose infinitesimal generator is a nonlocal Jacobi operator. Due to the lack of analytical results, we address the problem using a discretization scheme for simulating the trajectories of jump diffusion processes with state-dependent jumps in both frequency and amplitude. We obtain numerical approximations on their first passage time probability density functions and results for the qualitative behavior of other statistics of this random variable. Finally, we provide two examples of application of the method for different choices of the distribution involved in the mechanism of generation of the jumps.

1. Introduction

Recently, there has been growing interest in jump diffusion models in many applied areas, ranging from computational neuroscience [1,2,3] to mathematical biology [4], metrology [5] and queueing theory [6], just to name a few. In particular, they have been popular in financial modeling, starting with the celebrated paper by Merton [7]. Since this, the use of such models has been increasing in real markets and theoretical studies (see, for instance, [8,9,10,11,12]), thanks to their ability to account for some empirically observed effects that otherwise would not be explained by traditional diffusion-based models. A comprehensive discussion on this matter can be found in [13]. Roughly speaking, by choosing the parameters of the jump process appropriately, one can generate a wide variety of dynamics incorporating relevant effects without relying only on a very large amount of noise. In all these application contexts, it is often required to face the problem of the first passage time (FPT) of the process describing the dynamics of the model through a boundary [14,15]. Depending on the context, this crossing is interpreted in a different way, but from a mathematical point of view, its treatment is formally the same. Despite being a classical problem [16], its resolution is non-trivial, and exact analytical results are available only in a few cases even for pure diffusion processes [17,18].
Early attempts to introduce jumps occurring at exponential times can be found in [19,20,21,22,23,24], where to maintain mathematical tractability, the jumps were assumed to be of a constant amplitude or coming from a fixed distribution. More recently, results on generalized mechanisms of the generation of jumps have appeared, assuming that the jump size depends on the value of the process [25,26] or the state dependence is for both size and frequency [27]. Following this direction, we consider a family of processes with state-dependent jumps whose diffusion part evolves according to a Jacobi (or Wright–Fisher) model. In [28], the authors introduced and studied nonlocal Jacobi operators, which generalized the classical (local) Jacobi operators. Apart from some analytical results obtained in [27], the literature on the FPT problem for these processes is scarce. For this reason, the study of approximations and simulations of the involved quantities constitutes a fundamental tool.
In this paper, we use a discretization scheme for simulating sample paths of these jump diffusion processes with state-dependent jump intensities. At each time step of the algorithm, a downward jump can occur with a probability and amplitude that depend on the distribution characterizing the jump component and the actual state of the process. Between each jump epoch, the dynamics of the constructed process are purely diffusive and are simulated using the Milstein’s discretization method [29].
From the simulations of the trajectories of the process, we obtain approximations of the density of the FPT through a constant boundary for different choices of the measure describing the distribution of the jumps. In particular, in the case of the Jacobi process with jumps, we consider exponential and Pareto distributions. From the simulations, we observe that despite the presence of only downward jumps, the decay of the tails of the FPT pdfs is fast, as it happens for the diffusion processes without jumps. Moreover, under specific initial conditions, we might observe a bimodal FPT pdf. This behavior suggests the existence of interactions between the two components of the dynamics, resulting in a mixture of two distributions. From the reiteration of the simulation procedure, we also studied the behavior of some quantities related to the FPT, namely the mean and the variance of the FPT as a function of the parameters characterizing the jump component and the average number of jumps, and  we observed a non-trivial behavior for the average jump size, which had a non-monotone behavior as a result of the state dependence of the jumps for both frequency and amplitude.
Throughout this paper, we use the description of the process involving an infinitesimal generator. This approach is more convenient due to the fact that the usual formulation using a stochastic differential equation is made less intuitive by the presence of possibly complicated random measures. Moreover, the simulation strategy adopted here only needs the knowledge of a distribution Π , which will be defined later, and the values of the drift and diffusion parameters. Even the Bernstein function associated with the generator, which is essential for the calculation of the analytical results, does not need to be known explicitly.
The algorithm that we present is specialized for a non-local Jacobi operator, but since the continuous part of the trajectory is constructed using a classical discretization scheme, the procedure can be applied to other non-local perturbations of classical operators, as long as it is proven that the resulting operator is still the generator of a Markov process.
This paper is structured as follows. In Section 2, we introduce the Jacobi processes with jumps and the qualitative behavior of their trajectories. In Section 3, we present the problem of the FPT through the analytical results available in the literature. Section 4 is devoted to the description of the numerical procedure that simulates the sample paths of the process. In Section 5 and Section 6, we provide two examples of use of the discretization scheme for different choices of the distribution involved in the mechanism of generation of the jumps. Finally, in Section 7, we discuss the obtained results and highlight possible future works.

2. Jacobi Processes with Jumps and Their First Passage Time Problems

Let us denote with Y = ( Y t ) t 0 the generalized Jacobi process with jumps introduced in [28]. For this kind of Levy-type process, the best description is given in terms of infinitesimal generators, which are functional operators whose terms contain the drift part of the process, the diffusion component and the contribution of the jumps using the integral with respect to a random measure (for a complete description of Feller semigroups and generators, see, for instance, [30]). The process Y is a Feller process on [ 0 , 1 ] , whose infinitesimal generator is given for a smooth function f on [ 0 , 1 ] by
J Y f ( y ) = J f ( y ) + 0 ( f ( e r y ) f ( y ) ) Π ( d r ) y jump part
where J is the classical Jacobi operator
J f ( y ) = σ 2 2 y ( 1 y ) f ( y ) diffusion part λ y μ f ( y ) drift part
with σ 2 > 0 , where Π is a finite, nonnegative Radon measure on R + with = 0 r Π ( d r ) < .
Throughout this paper, we impose the following assumption that guarantees that y = 0 is an entrance boundary:
μ > + σ 2 2 .
The latter condition is a standing assumption in [28] and extends the Feller classification of boundaries in the presence of jumps. The results presented here could be obtained with an arbitrary set of parameters satisfying these conditions.
In [28], it is shown that J Y , which is obtained as a nonlocal perturbation of the generator of the classical Jacobi process, is indeed the generator of a Markov process on [ 0 , 1 ] with càdlàg trajectories. The hypotheses on the Π measure guarantee that J Y satisfies the positive maximum principle which, together with the Hille–Yosida–Ray theorem for Markov generators, ensures that J Y is the infinitesimal generator of a Markov semigroup on C 1 [ 0 , 1 ] (for more technical details, see [28]). As a consequence, the jumps are only downward, and both the amplitude and the intensity of the jumps are state-dependent. In fact, the process jumps from state y to state e r y at a frequency given by Π ( d r ) / y , which is inversely proportional to the achieved state. When the process is close to the lower boundary (i.e., zero), the average number of jumps is high, but the corresponding jump size is small, as the support of the distribution of the amplitude of the jumps is [ 0 , y ] . Conversely, for higher values of the state of the process, the probability of jumping becomes smaller, whereas the average jump size depends on Π . See Figure 1 for an example of two possible paths of the process under investigation.
Since this process can perform a finite number of jumps in a finite time, we can derive a path interpretation of this Markov process (see [31]). The stochastic process Y starts from y 0 by undergoing the same dynamics of the classical Jacobi process until a random time T , at which the process performs a downward jump. The survival probability up to time t of T is given by
P ( T > t ) = e Π ( R + ) Y t .
After a jump, the process restarts from the new position, undergoing the diffusion dynamics until the next jump.
Different choices of Π allow different sizes for the jumps. The more mass Π concentrates around zero, the smaller, in principle, the amplitude of the jumps is. If  Π admits large values with high probability, then the state of the process can almost be set to zero after the jump.

Bernstein and Bernstein–Gamma Functions

In this section, we recall a few definitions and results that will be useful in the following.
We recall that a function ϕ : [ 0 , ) [ 0 , ) is a Bernstein function if it is infinitely differentiable on R + and ( 1 ) n + 1 d n d u n ϕ ( u ) 0 for all n = 1 , 2 , and u 0 [32].
We observe that J Y is uniquely determined by σ 2 , Π , μ and λ . In particular, by fixing λ , the triplet ( σ 2 , Π , μ ) constitutes a Lévy triplet of the Bernstein function ϕ defined, for u 0 , by
ϕ ( u ) = u + 2 σ 2 μ 2 σ 2 1 + 2 σ 2 0 ( 1 e u r ) Π ¯ ( r ) d r
where Π ¯ ( r ) = r Π ( d u ) such that for a fixed λ , there is a one-to-one correspondence between ϕ and J Y (see [27] for more details).
In [33], the authors wrote W ϕ for the solution in the space of positive definite functions for the recurrence equation W ϕ ( z + 1 ) = ϕ ( z ) W ϕ ( z ) , with W ϕ ( 1 ) = 1 , z C and Re ( z ) > 0 . For any n N , we set for the Bernstein-Gamma function
W ϕ ( n + 1 ) = k = 1 n ϕ ( k )
with the convention k = 1 0 ϕ ( k ) = 1 . Note that the gamma function appears as a special case of the Bernstein–gamma function W ϕ for ϕ ( n ) = n .
Using this function, it is possible to introduce the mapping
2 F 1 a , b , ϕ ; x = n = 0 ( a ) n ( b ) n n ! x n W ϕ ( n + 1 )
with ( a ) n = Γ ( a + n ) Γ ( a ) , n N { 0 } and a C , which generalizes the Gauss hypergeometric function that appears as a special case for Π 0 (see [27] for more details).
We are now ready to formulate the problem of the first passage time of Y through a constant boundary and resume the existing analytical results.

3. The First Passage Time Problem

Let us consider the evolution of the stochastic process Y in the presence of a constant threshold S ( 0 , 1 ) . We are interested in the random time in which the process reaches the threshold S for the first time (i.e., the random variable):
T S = inf { t > 0 ; Y t S | Y 0 = y 0 < S } .
The direct problem of the first passage time consists mainly of finding the distribution of T S . Although it is a classical and easy-to-state problem, its solution is, for most of the stochastic processes, not available [15,16]. An analytical closed-form expression for the probability density function g ( t ) of T S is not known even for the classical Jacobi process [34]. Often, it is convenient to evaluate the Laplace transform of g ( t ) in order to obtain information on the distribution, the probability of crossing the threshold and the moments of T S . For the Jacobi process with jumps, the Laplace transform of g ( t ) is known to be [27], for any 0 < y 0 < S < 1 and q > 0 , the following:
E y 0 e q T S = 2 F 1 κ ( q ) , θ ( q ) ; ϕ ; y 0 2 F 1 κ ( q ) , θ ( q ) ; ϕ ; S ,
involving the mapping defined in Equation (6), where ϕ is a Bernstein function and κ ( q ) and θ ( q ) are solutions to the system
κ ( q ) θ ( q ) = 2 q σ 2 κ ( q ) + θ ( q ) + 1 = 2 λ σ 2 .
In principle, the moments of T S of any order can be computed using derivatives of E y 0 e q T S when they exist. The first moment is known to have the following explicit analytical expression [27]:
E y 0 [ T S ] = 2 σ 2 n = 0 ( 2 λ / σ 2 ) n n + 1 S n + 1 y 0 n + 1 W ϕ ( n + 2 ) .
However, the dependence of E y 0 [ T S ] on the parameters of the process is non-trivial since the contribution of the drift is hidden in the function W ϕ , which merges the contribution of the deterministic and diffusion components.
Using terminology that comes from the context of computational neuroscience, we distinguish between two possible regimes to characterize the tendency of the process to cross the barrier. If the asymptotic mean value of the process (The process has a stationary distribution, which is a generalized beta distribution.) is larger than S, then the process is in the so-called suprathreshold regime. In the classical case, in this regime, the crossings are regular, and the dynamics are driven mainly by the drift part. If the asymptotic mean is smaller than S, then the process is said to be in the subthreshold regime, and the noise plays a prominent role in the crossing of the threshold. The  Jacobi process with jumps is in the suprathreshold regime if
μ > a λ + 0 e r Π ¯ ( r ) d r .
The derivation of higher-order moments from the Laplace transform is impractical, involving at least second derivatives of the generalized Gaussian hypergeometric function in Equation (6) with respect to the parameters. For this reason, it is fundamental to construct an algorithm to simulate trajectories for the family of one-dimensional jump diffusion processes, with the state-dependent intensity generated by the functional in Equation (1) for different choices of Π .

4. The Discretization Scheme for the First Passage Time

We use a discretization scheme for simulating the sample paths of jump diffusion processes with state-dependent jump intensities. Due to the dependence of the jumps on the current state of the process, in terms of both frequency and amplitude, the times when the jumps occur cannot be drawn in advance in the simulation. For this reason, at each time step of the algorithm, a value for r is sampled from the distribution Π ( d r ) , and according to the probability in Equation (3), a jump from Y t to e r Y t may occur. Then, the trajectory moves according to the diffusion from state e r Y t if there was a jump; otherwise, it moves according to the diffusion from state Y t . Between the jump epochs, the dynamics of the constructed process are purely diffusive and are simulated using Milstein’s discretization method, which is a generalization of the Euler–Maruyama scheme used for stochastic processes with multiplicative noise [29].
When the intensity of the point process driving the jump component is state-dependent, the  error generated in the construction of the continuous component can be amplified in the simulation of the jumps and depends on the size of the time step. However, this algorithm is related to the discretization schemes for which it was proven in [35] that the method converges and the weak convergence order equals the order of the adopted time step. This means that the simulation error is of the same order as that of the discretization schemes for pure diffusions.
To avoid discretization errors, one can think of using an exact algorithm that simulates directly the hitting times without constructing the whole paths, as in [10] (see also [36,37]). For the special case of the Jacobi diffusion see [38].
However, the lack of knowledge for many functions and properties concerning these processes and the difficult implementation of the generalized functions involved in the transition density and the stationary distribution of the process may prevent the use of these exact strategies.
The proposed discretization strategy is light and very simple to implement, and it has an advantage: we can simulate the process just from the distribution Π without knowing explicitly the associated Bernstein function ϕ . Moreover, it relies on the study of the process in terms of its generator. Using generators in this context is more convenient with respect to the usual definition of the process as solution to a stochastic differential equation (SDE). Indeed, for a state-dependent jump process, the theory of SDEs is still incomplete and involves integrals with respect to some random measures that make both the numerical implementation and the interpretation of the dynamics hard.
In this paper, we focus mostly on the simulation of the paths of the process in order to answer the problem of the first passage time. For this reason, we construct the trajectory until it reaches the level S, and we record the time of this crossing. We repeat this procedure n times, and we use the collected FPT times to find an approximation of the FPT density and other statistics for which it is not possible to have analytical results.
In Algorithm 1 we illustrate a scheme of the sampling procedure we use to draw an FPT from the process.
Algorithm 1 Sampling FPT
  • Require: y 0 , S, d t , σ 2 , λ , μ
  • Ensure: FPT sample t S
  •   while  y i < S do
  •        r Π ( d r )
  •        j B e 1 e x p ( r y i 1 d t )
  •       if j = 1 then                               ▹ a jump occurs
  •            y * e r y i 1
  •       else if j = 0 then                           ▹ a jump does not occur
  •            y * y i 1
  •       end if
  •        g ( y * ) σ 2 2 y * ( 1 y * )
  •        y i y * ( λ y * μ ) d t + g ( y * ) Δ W i + 1 2 g ( y * ) ( g ( y * ) ) ( ( Δ W i ) 2 d t )     ▹ diffusion
  • end while
  • t S = i

5. Example: Exponential Distribution

We consider a parametric family of non-local Jacobi operators of the form in Equation (1), for which
Π ¯ ( r ) = r Π ( d u ) = e α r , r > 0
where
Π ( d r ) = α e α r d r
is an exponential probability density function with a parameter α . In particular, we will choose α 1 throughout the paper as in [28].
Moreover, to guarantee from the assumption in Equation (2) that y = 0 is an entrance boundary, we have
μ > σ 2 2 + 1 α .
The presence in the last inequality of a positive term 1 / α suggests that the maximum noise amplitude has to be smaller than in the classical case. A large value of σ can lead the process across the lower boundary, a condition that we want to avoid. Unfortunately, classical approximation schemes cannot preserve the properties of the boundaries independent of the choice of the time discretization step, even if the theoretical assumption in Equation (14) is satisfied [39,40]. In particular, for the Jacobi process, even the splitting methods do not preserve the boundary behavior, and other strategies such as the balanced implicit split step (BISS) method, which is able to preserve the boundary structure, are lacking in accuracy (see [41] for an extensive discussion).
Finally, in this case, the regime is a suprathreshold if
μ > a λ + 1 1 + α ,
We observe that, as expected, the downward jumps make the asymptotic mean of Y, μ λ 1 1 + α , smaller than that of the classical Jacobi process ( μ / λ ):
Remark 1.
The integro-differential operator J Y from Equation (1) takes the form
J f ( y ) = σ 2 2 y ( 1 y ) f ( y ) λ y μ f ( y ) 0 1 ( f ( r ) f ( y ) ) r α y α + 1 d r .
Moreover, = 1 / α , and simple algebra yields to the explicit expression of the corresponding Bernstein function:
ϕ ( u ) = u + 2 σ 2 μ 1 u + α 1 .
However, the application of the proposed discretization scheme does not require the knowledge of the explicit expression of the infinitesimal operator nor of the Bernstein function ϕ. This constitutes a great advantage when working with distributions whose expression prevents easy calculation of the involved quantities.
We want to investigate the behavior of the FPT density and other statistical quantities of this random variable under a change in the parameters characterizing the distribution Π . Precise analysis of these moments is made difficult by the presence of the generalized hypergeometric function in the expression of the Laplace transform (Equation (8)), while a formal analytical study on the FPT pdf is prevented by the lack of explicit results. For these reasons, the information obtained from the simulations are of great importance. In the following, we show some of the qualitative statistical behaviors of T S using the simulation scheme described in the previous section.
Let us consider the jump diffusion process Y generated by the non-local operator in Equation (1), with Π given in Equation (13).
In Figure 1, we show two possible realizations of Y until it reaches the threshold S for the first time. We can see the impact of different choices of α on the trajectories of the considered stochastic process. A lower value of α implies that possibly larger values r are sampled from Π , resulting in a lower frequency of the jumps in principle and a higher size for them. However, the state dependence of the jumps makes the jump generation mechanism more complex. Moreover, these pictures display how the jumps were more dense when the process was close to the lower boundary, as expected from the probability in Equation (3).
In [27], it was shown, using Equation (10), that the mean FPT decreases as α increases. This is due to the fact that for large values of α , a small value r is most likely sampled from Π , affecting the probability (Equation (3)). Using simulations, we observed the same behavior for the variance of T S as a function of α (not shown). This is explained equivalently by the fact that the average number of jumps decreased as α increased (see Figure 2). However, the dynamics was not as simple as it may appear due to the presence of state-dependent jumps. Indeed, it is interesting to observe the behavior of the average length of the jumps as a function of α in Figure 2. Up to some value of α , the average size of the jumps increased with α . This behavior could sound counterintuitive, since lower values of α imply, on average, a higher gap between y and e r y . However, we have to take into account that jumps were more frequent when y was small, and since the support of the distribution of the amplitude of the jumps was [ 0 , y ] , most of the jumps were short. Therefore, roughly speaking, even if there is a large jump that pushes the process to the lower boundary, it will be compensated for by many small jumps close to the zero level. Even for higher values of α , we observed that the average length of the jumps started to decrease. In this case, very large jumps occurred with a small probability, so on average, the size of the jumps decreased.
If a diffusion process admits a stationary distribution, then the corresponding FPT pdf is known to have asymptotically, for large times and large boundaries, an approximately exponential distribution whose mean is related to the average first passage time from the origin to the boundary (see [42] for the Ornstein–Uhlenbeck process, ref. [43] for one-dimensional diffusion processes and [44] for Gauss–Markov processes). Since the presence of downward jumps decreases the asymptotic mean of the process with respect to the case without jumps and increases the mean FTP, it is natural to ask whether the FPT pdf tail decays slower.
In Figure 3, we show approximations of the FPT pdfs obtained from histograms of 5 × 10 4 simulated first passage times of Y through S for 6 different choices of α . As expected, we can appreciate that the threshold S was more likely to be crossed earlier for larger values of α , but the tail decay remained qualitatively the same.
To measure the tail decay, we selected the densities g ˜ of the histograms after the third quartile, which were the heights of the bars on the right-hand side of the histogram that summed to a probability approximately equal to 0.25 , and fit them with two different models. More specifically, we used an exponential function and a power law function, defined by
M e : h = β 0 + β 1 e β 2 t
M p : h = β 0 + β 1 t β 2 .
The estimation of the parameters ( β 0 , β 1 , β 2 ) for the two models was performed using non-linear optimization software.
In the legends of Figure 3, we display the logarithm of the mean squared error (MSE) as a measure of goodness of fit for both curves. We can see that although both curves approximated the tails well, the approximation error of the exponential function was always the smallest. This might indicate that the tails of the FPT pdf had an exponential decay.
Interestingly, the FPT distribution could show a bimodal behavior when the starting point y 0 was close to the threshold S. This behavior followed from the fact that if no jump occurred in the very first moments, then the positive drift quickly pushed the process toward S (first peak). On the other hand, if a large jump took place before the process was absorbed right away in S, then the process would take a longer time to reach the threshold with a distribution showing a longer tail (second peak). In other words, one could see this distribution as a mixture of two distributions of the first passage time: one related to a Jacobi process without jumps and the other related to a Jacobi process with jumps and a starting point y 0 far from S, where the weights of such a mixture depend on the probability of the process to perform a jump before the drift pushes it across the threshold S. An example of this behavior can be appreciated in Figure 4, where to better stress the bimodality of the distribution, we display an histogram of the logarithm of the FPT in a setting where y 0 is close to S.

6. Example: Pareto Distribution

We consider a parametric family of non-local Jacobi operators as in Equation (1) for which
Π ¯ ( r ) = r Π ( d u ) = η r θ if r η , 1 if r < η
where
Π ( d r ) = θ η θ r θ + 1 d r if r η 0 if r < η
is a Pareto Type I probability density function with a shape parameter θ > 0 and location parameter η > 0 . In order to match the assumption that = 0 r Π ( d r ) < , we will choose θ > 1 throughout this paper. Moreover, to guarantee that y = 0 is an entrance boundary, from Equation (2), we have
μ > θ η θ 1 + σ 2 2 .
In this case, the regime is a suprathreshold if, for r η , the following is true:
μ > S λ + 0 e r Π ¯ ( d r ) = S λ + 0 e r η r θ d r = S λ + η θ Γ ( 1 θ ) ,
where the last equality involving the integral representation of the gamma function holds only for Re ( θ ) > 1 , which is a case that cannot be considered here.
An analytical expression of the moments of T S in the case of the Pareto distribution is not known and neither is the expression of the Bernstein function associated with the non-local operator. Using the discretization scheme, we can find an estimation of the mean and variance of T S as a function of the two parameters characterizing the Pareto distribution.
In Figure 5, we show the behavior of the variance of the FPT when tuning simultaneously the scale and location parameters θ and η , respectively. The variance increases with the location parameter η . In fact, the support of the Pareto distribution is the interval [ η , ] , meaning that a large value r [ η , ] will be sampled by the algorithm if η is large. At the same time, the variance of T S decreases as θ increases due to the shape of the distribution for large values of θ that favour small values of r. In order to match Equation (21) for all the couples ( θ , η ) , the chosen drift is relatively strong ( μ = 10.5 ) , which explains why the resulting variance was small in all the considered cases. Qualitatively, the same behavior can be observed for the mean FPT (not shown here), where the mean FPT increases with η and decreases as θ increases.
Under the assumption of θ > 2 , which guarantees finite variance for the Pareto distribution, we consider the two following parameter choices:
Case 1:
θ = 1 + 2 , η = 2 1 + 2 ;
Case 2:
θ = 1 + 2 , η = 2 2 ( 1 + 2 ) .
These choices guarantee expected values equal to 1 and 1 / 2 , respectively, and variances equal to the square of the means, as in the exponential case, allowing a comparison between the two examples.
In Figure 6, we consider the histograms of the FPT for the two cases. As performed in the previous section, we applied the fitting models of Equations (18) and (19) to the tails of the distributions. Additionally, in this case, both models fit well, but the use of an exponential curve resulted in a lower MSE. The Pareto distribution is a classical example of a heavy-tailed probability distribution, meaning that large values of r in the mechanism of generation of jumps can be chosen by the algorithm. It could be natural to expect that the shape of the FPT pdf could be stretched by the heavy tails of the Pareto distribution. However, for our parameter choice, the tail of the Pareto distribution became fatter than the one of the exponential distribution only for values with negligible probabilities.

7. Discussion

Due to the lack of analytical results regarding the FPT of jump diffusion processes, for which the jumps are state-dependent in terms of both frequency and amplitude, the use of simulations is crucial. Using a discretization scheme, we simulated the trajectories of these processes, and we studied the problem of their passage times through a constant boundary. The method is specialized for a Jacobi process with jumps but can be used for any Markov process whose infinitesimal generator is obtained as a non-local perturbation of a classical operator. We obtained approximations on the FPT pdf for different choices of the measure describing the distribution of the jumps. In particular, in the case of the Jacobi process with jumps, we considered exponential and Pareto distributions. From the simulations, we observed that despite the presence of only downward jumps, the decay of the tails of the FPT pdfs was fast, as happened for the diffusion processes without jumps. This was a consequence of the main assumption in Equation (2) for the drift of the process, which guaranteed the Markov property. Analytical results for the pdf’s tails’ decay could be obtained by Tauberian theorems for Laplace transforms, but the presence of the generalized functions in Equation (8) prevented straightforward calculations.
Another interesting feature that might be observed is the multimodality of the FPT pdf, which can appear even if the drift part of the process is non-periodic. The effect is more visible if the starting point of the process y 0 is close enough to the boundary S. In this case, the first peak of the FPT pdf is determined by the drift part of the process, and a second bump is visible, suggesting the existence of interactions between the components of the dynamics, resulting in a mixture of two distributions. In [2], a similar behavior was observed but only in the presence of both positive and negative jumps.
From iterations of the simulation procedure, we also studied the behavior of some moments of the FPT. In particular, we studied the variance of the FPT as a function of the parameters characterizing the jump component in the case of both one- and two-parameter distributions. Finally, we studied the average number of jumps as a function of the parameter of Π , and we observed a non-trivial behavior of the average jump size, which had a non-monotone behavior as a result of the state dependence of the jumps in terms of both frequency and amplitude.
We stress that we used the infinitesimal generator in Equation (1), but in the application context, more specific operators defined on different intervals can be used. The same result follows if one can identify a homeomorphism between the new semigroups and the one of Jacobi processes with jumps on (0, 1) defined in Equation (1). For example, this was performed in the framework of mathematical neuroscience in [27], taking advantage of the intertwining approach.

Author Contributions

Conceptualization, G.D. and A.L.; Methodology, G.D. and A.L.; Software, G.D. and A.L.; Investigation, G.D. and A.L.; Writing—original draft, G.D. and A.L.; Writing—review & editing, G.D. and A.L.; Funding acquisition, G.D. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA—INdAM).

Data Availability Statement

The codes used during this study are available from the authors on reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
FPTFirst passage time
pdfProbability density function
MSEMean squared error
SDEStochastic differential equation

References

  1. Jahn, P.; Berg, R.W.; Hounsgaard, J.; Ditlevsen, S. Motoneuron membrane potentials follow a time inhomogeneous jump diffusion process. J. Comput. Neurosci. 2011, 31, 563–579. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Sirovich, R.; Sacerdote, L. Noise induced phenomena in jump diffusion models for single neuron spike activity. IEEE Int. Jt. Conf. Neural Netw. 2004, 4, 3025–3028. [Google Scholar]
  3. Sirovich, R.; Sacerdote, L.; Villa, A.E.P. Cooperative behavior in a jump diffusion model for a simple network of spiking neurons. Math. Biosci. Eng. 2013, 11, 385–401. [Google Scholar] [CrossRef] [PubMed]
  4. Giorno, V.; Román-Román, P.; Spina, S.; Torres-Ruiz, F. Estimating a non-homogeneous Gompertz process with jumps as model of tumor dynamics. Comput. Stat. Data Anal. 2017, 107, 18–31. [Google Scholar] [CrossRef]
  5. Zucca, C.; Tavella, P. A mathematical model for the atomic clock error in case of jumps. Metrologia 2015, 52, 514. [Google Scholar] [CrossRef] [Green Version]
  6. Dharmaraja, S.; Di Crescenzo, A.; Giorno, V.; Nobile, A.G. A continuous-time Ehrenfest model with catastrophes and its jump-diffusion approximation. J. Stat. Phys. 2015, 161, 326–345. [Google Scholar] [CrossRef]
  7. Merton, R.C. Option pricing when underlying stock returns are discontinuous. J. Financ. Econ. 1976, 3, 125–144. [Google Scholar] [CrossRef] [Green Version]
  8. Briani, M.; Caramellino, L.; Terenzi, G. Convergence Rate of Markov Chains and Hybrid Numerical Schemes to Jump-Diffusion with Application to the Bates Model. SIAM J. Numer. Anal. 2021, 59, 477–502. [Google Scholar] [CrossRef]
  9. Brignone, R.; Sgarra, C. Asian options pricing in Hawkes-type jump-diffusion models. Ann. Financ. 2020, 16, 101–119. [Google Scholar] [CrossRef]
  10. Casella, B.; Roberts, G.O. Exact simulation of jump-diffusion processes with Monte Carlo applications. Methodol. Comput. Appl. Probab. 2011, 13, 449–473. [Google Scholar] [CrossRef] [Green Version]
  11. Kou, S.G. A jump-diffusion model for option pricing. Manag. Sci. 2002, 48, 1086–1101. [Google Scholar] [CrossRef]
  12. Ksendal, B.; Sulem, A. Applied Stochastic Control of Jump Diffusions; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  13. Tankov, P. Financial Modelling with Jump Processes; Chapman and Hall/CRC: Boca Raton, FL, USA, 2003. [Google Scholar]
  14. Metzler, R.; Redner, S.; Oshanin, G. First-Passage Phenomena and Their Applications; World Scientific: Singapore, 2014. [Google Scholar]
  15. Redner, S. A Guide to First-Passage Processes; Cambridge University Press: Cambridge, UK, 2001. [Google Scholar]
  16. Karlin, S.; Taylor, H.E. A Second Course in Stochastic Processes; Elsevier: Amsterdam, The Netherlands, 1981. [Google Scholar]
  17. Buonocore, A.; Caputo, L.; D’Onofrio, G.; Pirozzi, E. Closed-form solutions for the first-passage-time problem and neuronal modeling. Ric. Mat. 2015, 64, 421–439. [Google Scholar] [CrossRef] [Green Version]
  18. Ricciardi, L.M.; Di Crescenzo, A.; Giorno, V.; Nobile, A.G. An outline of theoretical and algorithmic approaches to first passage time problems with applications to biological modeling. Math. Jpn. 1999, 50, 247–322. [Google Scholar]
  19. Abundo, M. On first-passage times for one-dimensional jump-diffusion processes. Probab. Math. Stat. 2000, 20, 399–423. [Google Scholar]
  20. Abundo, M. On the first hitting time of a one-dimensional diffusion and a compound Poisson process. Methodol. Comput. Appl. Probab. 2010, 12, 473–490. [Google Scholar] [CrossRef]
  21. Atiya, A.F.; Metwally, S.A. Efficient estimation of first passage time density function for jump-diffusion processes. SIAM J. Sci. Comput. 2005, 26, 1760–1775. [Google Scholar] [CrossRef]
  22. Herrmann, S.; Massin, N. Exact simulation of the first passage time through a given level of jump diffusions. Math. Comput. Simul. 2023, 203, 553–576. [Google Scholar] [CrossRef]
  23. Kou, S.G.; Wang, H. First passage times of a jump diffusion process. Adv. Appl. Probab. 2003, 35, 504–531. [Google Scholar] [CrossRef]
  24. Xie, J.; Cui, Z.; Zhang, Z. Some new infinite series expansions for the first passage time densities in a jump diffusion model with phase-type jumps. Appl. Math. Comput. 2022, 429, 127251. [Google Scholar] [CrossRef]
  25. Lefebvre, M. The ruin problem for a Wiener process with state-dependent jumps. J. Appl. Math. Stat. Inform. 2020, 16, 13–23. [Google Scholar] [CrossRef]
  26. Lefebvre, M. First-passage problems for diffusion processes with state-dependent jumps. Commun. Stat.-Theory Methods 2022, 51, 2908–2918. [Google Scholar] [CrossRef]
  27. D’Onofrio, G.; Patie, P.; Sacerdote, L. Jacobi processes with jumps as neuronal models: A first passage time analysis. arXiv 2022, arXiv:2205.08237. [Google Scholar]
  28. Cheridito, P.; Patie, P.; Srapionyan, A.; Vaidyanathan, A. On non-local ergodic Jacobi semigroups: Spectral theory, convergence-to-equilibrium and contractivity. J. l’Ecole Polytech.-Math. 2021, 8, 331–378. [Google Scholar] [CrossRef]
  29. Higham, D.J. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev. 2001, 43, 525–546. [Google Scholar] [CrossRef]
  30. Ethier, S.N.; Kurtz, T.G. Markov Processes: Characterization and Convergence; John Wiley & Sons: Hoboken, NJ, USA, 2009. [Google Scholar]
  31. Bass, R.F. Adding and subtracting jumps from Markov processes. Trans. Am. Math. Soc. 1979, 255, 363–376. [Google Scholar] [CrossRef]
  32. Schilling, R.L.; Song, R.; Vondraček, Z. Bernstein Functions; Volume 37 of de Gruyter Studies in Mathematics; Theory and Applications; Walter de Gruyter & Co.: Berlin, Germany, 2010. [Google Scholar]
  33. Patie, P.; Savov, M. Bernstein-gamma functions and exponential functionals of Lévy processes. Electron. J. Probab. 2018, 23, 1–101. [Google Scholar] [CrossRef]
  34. D’Onofrio, G.; Tamborrino, M.; Lansky, P. The Jacobi diffusion process as a neuronal model. Chaos 2018, 28, 103119. [Google Scholar] [CrossRef] [Green Version]
  35. Glasserman, P.; Merener, N. Convergence of a discretization scheme for jump-diffusion processes with state–dependent intensities. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 2004, 460, 111–127. [Google Scholar] [CrossRef]
  36. Beskos, A.; Roberts, G.O. Exact simulation of diffusions. Ann. Appl. Probab. 2005, 15, 2422–2444. [Google Scholar] [CrossRef] [Green Version]
  37. Herrmann, S.; Zucca, C. Exact simulation of the first-passage time of diffusions. J. Sci. Comput. 2019, 79, 1477–1504. [Google Scholar] [CrossRef] [Green Version]
  38. Jenkins, P.A.; Spano, D. Exact simulation of the Wright–Fisher diffusion. Ann. Appl. Probab. 2017, 27, 1478–1509. [Google Scholar] [CrossRef] [Green Version]
  39. Tubikanec, I.; Tamborrino, M.; Lansky, P.; Buckwar, E. Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion. J. Comput. Appl. Math. 2022, 406, 113951. [Google Scholar] [CrossRef]
  40. Dangerfield, C.E.; Kay, D.; Macnamara, S.; Burrage, K. A boundary preserving numerical algorithm for the Wright–Fisher model with mutation. BIT Numer. Math. 2012, 52, 283–304. [Google Scholar] [CrossRef] [Green Version]
  41. Eder, J. Boundary Behaviour of Pearson Diffusion Processes and Numerical Splitting Methods Preserving Them. Ph.D. Thesis, Universität Linz, Linz, Austria, 2020. [Google Scholar]
  42. Nobile, A.G.; Ricciardi, L.M.; Sacerdote, L. Exponential trends of Ornstein–Uhlenbeck first-passage-time densities. J. Appl. Probab. 1985, 22, 360–369. [Google Scholar] [CrossRef]
  43. Giorno, V.; Nobile, A.G.; Ricciardi, L.M. On the asymptotic behaviour of first-passage-time densities for one-dimensional diffusion processes and varying boundaries. Adv. Appl. Probab. 1990, 22, 883–914. [Google Scholar] [CrossRef]
  44. Nobile, A.G.; Pirozzi, E.; Ricciardi, L.M. Asymptotics and evaluations of FPT densities through varying boundaries for Gauss-Markov processes. Sci. Math. Jpn. 2008, 67, 241–266. [Google Scholar]
Figure 1. Examples of trajectories of Y for α = 1 (left) and α = 2 (right) while fixing the other parameters. At the bottom of each plot, each red vertical segment represents the time point of a jump.
Figure 1. Examples of trajectories of Y for α = 1 (left) and α = 2 (right) while fixing the other parameters. At the bottom of each plot, each red vertical segment represents the time point of a jump.
Fractalfract 07 00030 g001
Figure 2. Average jump lengths (black solid line) and average number of jumps (red dashed line) as a function of α . Confidence bands for the jump lengths (gray solid line) were obtained for the 5th and the 95th percentiles. The averages were taken over 5000 simulated sample paths of Y from Equation (1) for μ = 1.1 , σ = 0.1 , λ = 1.1 , x 0 = 5 × 10 4 , S = 0.9 and time step d t = 0.01 .
Figure 2. Average jump lengths (black solid line) and average number of jumps (red dashed line) as a function of α . Confidence bands for the jump lengths (gray solid line) were obtained for the 5th and the 95th percentiles. The averages were taken over 5000 simulated sample paths of Y from Equation (1) for μ = 1.1 , σ = 0.1 , λ = 1.1 , x 0 = 5 × 10 4 , S = 0.9 and time step d t = 0.01 .
Fractalfract 07 00030 g002
Figure 3. Histograms of FPT data obtained from 5 × 10 4 simulations of trajectories of Y for 6 different choices of α . Other parameters were μ = 1.1 , σ = 0.1 , λ = 1.1 , y 0 = 5 × 10 5 , S = 0.9 and d t = 0.01 . The black vertical line indicates the third quartile, marking the beginning of the tail of the approximated FPT pdf. The legend shows the MSE in logarithmic scale of the estimation of the tails made with exponential and power curves.
Figure 3. Histograms of FPT data obtained from 5 × 10 4 simulations of trajectories of Y for 6 different choices of α . Other parameters were μ = 1.1 , σ = 0.1 , λ = 1.1 , y 0 = 5 × 10 5 , S = 0.9 and d t = 0.01 . The black vertical line indicates the third quartile, marking the beginning of the tail of the approximated FPT pdf. The legend shows the MSE in logarithmic scale of the estimation of the tails made with exponential and power curves.
Fractalfract 07 00030 g003
Figure 4. Histogram of the logarithm of FPT data obtained from 5000 simulations of trajectories of Y, with α = 1 and y 0 = 0.85 . All other parameters were chosen as in Figure 3.
Figure 4. Histogram of the logarithm of FPT data obtained from 5000 simulations of trajectories of Y, with α = 1 and y 0 = 0.85 . All other parameters were chosen as in Figure 3.
Fractalfract 07 00030 g004
Figure 5. Variance of the FPT as a function of the scale and location parameters of the Pareto distribution. The heatmap was obtained by simulating 2000 FPTs for each couple of parameters. The trajectories were obtained from Equations (1) and (20) with μ = 10.5 , σ = 0.1 , λ = 1.1 , x 0 = 5 × 10 5 , S = 0.9 and d t = 0.01 .
Figure 5. Variance of the FPT as a function of the scale and location parameters of the Pareto distribution. The heatmap was obtained by simulating 2000 FPTs for each couple of parameters. The trajectories were obtained from Equations (1) and (20) with μ = 10.5 , σ = 0.1 , λ = 1.1 , x 0 = 5 × 10 5 , S = 0.9 and d t = 0.01 .
Fractalfract 07 00030 g005
Figure 6. Histograms of 2 × 10 4 simulated FPT data. The parameters of the Pareto distribution were those of Case 1 (left) and Case 2 (right). All other parameters were chosen as in Figure 3. The black vertical line indicates the third quartile, marking the beginning of the tail of the approximated FPT pdf. The legend shows the MSE in a logarithmic scale of the estimation of the tails made with exponential and power curves.
Figure 6. Histograms of 2 × 10 4 simulated FPT data. The parameters of the Pareto distribution were those of Case 1 (left) and Case 2 (right). All other parameters were chosen as in Figure 3. The black vertical line indicates the third quartile, marking the beginning of the tail of the approximated FPT pdf. The legend shows the MSE in a logarithmic scale of the estimation of the tails made with exponential and power curves.
Fractalfract 07 00030 g006
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

D’Onofrio, G.; Lanteri, A. Approximating the First Passage Time Density of Diffusion Processes with State-Dependent Jumps. Fractal Fract. 2023, 7, 30. https://doi.org/10.3390/fractalfract7010030

AMA Style

D’Onofrio G, Lanteri A. Approximating the First Passage Time Density of Diffusion Processes with State-Dependent Jumps. Fractal and Fractional. 2023; 7(1):30. https://doi.org/10.3390/fractalfract7010030

Chicago/Turabian Style

D’Onofrio, Giuseppe, and Alessandro Lanteri. 2023. "Approximating the First Passage Time Density of Diffusion Processes with State-Dependent Jumps" Fractal and Fractional 7, no. 1: 30. https://doi.org/10.3390/fractalfract7010030

APA Style

D’Onofrio, G., & Lanteri, A. (2023). Approximating the First Passage Time Density of Diffusion Processes with State-Dependent Jumps. Fractal and Fractional, 7(1), 30. https://doi.org/10.3390/fractalfract7010030

Article Metrics

Back to TopTop