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
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
, whose infinitesimal generator is given for a smooth function
f on
by
where
is the classical Jacobi operator
with
, where
is a finite, nonnegative Radon measure on
with
.
Throughout this paper, we impose the following assumption that guarantees that
is an entrance boundary:
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
, which is obtained as a nonlocal perturbation of the generator of the classical Jacobi process, is indeed the generator of a Markov process on
with càdlàg trajectories. The hypotheses on the
measure guarantee that
satisfies the positive maximum principle which, together with the Hille–Yosida–Ray theorem for Markov generators, ensures that
is the infinitesimal generator of a Markov semigroup on
(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
at a frequency given by
, 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
. 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
by undergoing the same dynamics of the classical Jacobi process until a random time
, at which the process performs a downward jump. The survival probability up to time
t of
is given by
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
is a Bernstein function if it is infinitely differentiable on
and
for all
and
[
32].
We observe that
is uniquely determined by
and
. In particular, by fixing
, the triplet
constitutes a Lévy triplet of the Bernstein function
defined, for
, by
where
such that for a fixed
, there is a one-to-one correspondence between
and
(see [
27] for more details).
In [
33], the authors wrote
for the solution in the space of positive definite functions for the recurrence equation
with
,
and Re
. For any
, we set for the Bernstein-Gamma function
with the convention
. Note that the gamma function appears as a special case of the Bernstein–gamma function
for
.
Using this function, it is possible to introduce the mapping
with
and
, which generalizes the Gauss hypergeometric function that appears as a special case for
(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
. We are interested in the random time in which the process reaches the threshold
S for the first time (i.e., the random variable):
The direct problem of the first passage time consists mainly of finding the distribution of
. 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
of
is not known even for the classical Jacobi process [
34]. Often, it is convenient to evaluate the Laplace transform of
in order to obtain information on the distribution, the probability of crossing the threshold and the moments of
. For the Jacobi process with jumps, the Laplace transform of
is known to be [
27], for any
and
, the following:
involving the mapping defined in Equation (
6), where
is a Bernstein function and
and
are solutions to the system
In principle, the moments of
of any order can be computed using derivatives of
when they exist. The first moment is known to have the following explicit analytical expression [
27]:
However, the dependence of on the parameters of the process is non-trivial since the contribution of the drift is hidden in the function , 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
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
, and according to the probability in Equation (
3), a jump from
to
may occur. Then, the trajectory moves according to the diffusion from state
if there was a jump; otherwise, it moves according to the diffusion from state
. 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:, S, , , , Ensure: FPT sample while do if j = 1 then ▹ a jump occurs else if j = 0 then ▹ a jump does not occur end if ▹ diffusion end while
|
5. Example: Exponential Distribution
We consider a parametric family of non-local Jacobi operators of the form in Equation (
1), for which
where
is an exponential probability density function with a parameter
. In particular, we will choose
throughout the paper as in [
28].
Moreover, to guarantee from the assumption in Equation (
2) that
is an entrance boundary, we have
The presence in the last inequality of a positive term
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
We observe that, as expected, the downward jumps make the asymptotic mean of Y, , smaller than that of the classical Jacobi process ():
Remark 1. The integro-differential operator from Equation (1) takes the form Moreover, , and simple algebra yields to the explicit expression of the corresponding Bernstein function: 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
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
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
. 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
, 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
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
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
, and fit them with two different models. More specifically, we used an exponential function and a power law function, defined by
The estimation of the parameters 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
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
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
is close to
S.
6. Example: Pareto Distribution
We consider a parametric family of non-local Jacobi operators as in Equation (
1) for which
where
is a Pareto Type I probability density function with a shape parameter
and location parameter
. In order to match the assumption that
, we will choose
throughout this paper. Moreover, to guarantee that
is an entrance boundary, from Equation (
2), we have
In this case, the regime is a suprathreshold if, for
, the following is true:
where the last equality involving the integral representation of the gamma function holds only for Re
, which is a case that cannot be considered here.
An analytical expression of the moments of 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 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
will be sampled by the algorithm if
is large. At the same time, the variance of
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
, 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 , which guarantees finite variance for the Pareto distribution, we consider the two following parameter choices:
- Case 1:
, ;
- Case 2:
, .
These choices guarantee expected values equal to 1 and , 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
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.