On the Estimation of the Persistence Exponent for a Fractionally Integrated Brownian Motion by Numerical Simulations

Abstract

For a fractionally integrated Brownian motion (FIBM) of order $\alpha \in (0,1],$$X_{\alpha }\left(t\right),$ we investigate the decaying rate of $P({\tau }_S^{\alpha }>t)$ as $t\to +\infty ,$ where ${\tau }_S^{\alpha }=inf\{t>0:X_{\alpha }\left(t\right)\ge S\}$ is the first-passage time (FPT) of $X_{\alpha }\left(t\right)$ through the barrier $S>0.$ Precisely, we study the so-called persistent exponent $\theta =\theta \left(\alpha \right)$ of the FPT tail, such that $P({\tau }_S^{\alpha }>t)=t^{-\theta +o\left(1\right)},$ as $t\to +\infty ,$ and by means of numerical simulation of long enough trajectories of the process $X_{\alpha }\left(t\right),$ we are able to estimate $\theta \left(\alpha \right)$ and to show that it is a non-increasing function of $\alpha \in (0,1],$ with $1/4\le \theta \left(\alpha \right)\le 1/2.$ In particular, we are able to validate numerically a new conjecture about the analytical expression of the function $\theta =\theta \left(\alpha \right),$ for $\alpha \in (0,1].$ Such a numerical validation is carried out in two ways: in the first one, we estimate $\theta \left(\alpha \right),$ by using the simulated FPT density, obtained for any $\alpha \in (0,1];$ in the second one, we estimate the persistent exponent by directly calculating $P\left({max}_{0\le s\le t}{X}_{\alpha }\left(s\right)<1\right).$ Both ways confirm our conclusions within the limit of numerical approximation. Finally, we investigate the self-similarity property of $X_{\alpha }\left(t\right)$ and we find the upper bound of its covariance function.

1. Introduction

The study of integrals with respect to the time-parameter of assigned stochastic processes constitutes one of the main chapters of stochastic calculus and one of the main tools for designing phenomenological models (see, for instance, [1,2,3] and references therein). Fractionally integrated stochastic processes are the natural extensions of the above processes in the context of the fractional calculus applied to the stochastic one (see, for instance, [4,5]). They are a rather new topic which appears to be of interest both from a theoretical point of view and for application (see [1] and references therein). Here, in particular, we focus on the (Liouville) fractionally integrated Brownian motion (FIBM) of order $\alpha \in (0,1],$ denoted by $X_{\alpha }\left(t\right)$ rigorously defined below. Our aim is to study the distribution of ${\tau }_S^{\alpha }=inf\{t>0:X_{\alpha }\left(t\right)\ge S\},$ also denoted only by $\tau$ if there is no ambiguity about the specified boundary S. This is the first-passage time (FPT) of $X_{\alpha }\left(t\right)$ through a boundary (otherwise called a barrier or threshold) $S>0.$ We specifically address the problem of studying the decaying rate of $P({\tau }_S^{\alpha }>t)$ as $t\to +\infty .$

The study of the distribution of the FPT $\tau$ of a stochastic process through a boundary is a classic problem in probability theory; generally, it is difficult to obtain an explicit expression of this law. However, it has been observed that in many interesting cases, the survival function has a polynomial decay which does not depend on the boundary:

$$P(\tau >t)=t^{-\theta +o\left(1\right)},\mathrm{as}t\to +\infty ,$$

(1)

or, equivalently:

$$\underset{t\to +\infty }{lim}\frac{logP(\tau >t)}{logt}=-\theta ,$$

(2)

where $\theta$ is a positive constant called persistence exponent and characterizes the profile of the tail of the distribution of $\tau$ for large t values. The computation of this exponent turns out to have connections with various problems in probability and mathematical physics. In general, for self-similar processes, the persistent exponent $\theta$ does not depend on the boundary, so this fact is natural in our case, since $X_{\alpha }\left(t\right)$ is actually self-similar (see Section 4.1).

A well-known result concerns Brownian motion: in this case, the persistence exponent turns out to be equal to $1/2;$ another important result is that of Goldman-Sinai, regarding the case of integrated Brownian motion (see [6,7]), for which the persistence exponent is $\theta =1/4.$ A generalization of this result regards the study of the persistence exponent for twice integrated, or more generally n-th time integrated, Brownian motion (see [3,8] and the references therein). Another example is the study of the persistence exponent for integrated fractional Brownian motion with Hurst parameter H (it was conjectured in [9] that $\theta$ should be $H(1-H))$. Moreover, the persistence exponent for the integrated stable Lévy process was studied in [8].

Furthermore, the persistence exponent was studied for an $\alpha$—fractionally integrated centered Lévy process; in [10], it was proved that the corresponding persistence exponent is a non-increasing function of the fractional order $\alpha ;$ the class of processes considered includes FIBM.

1.1. The Motivation

To our knowledge, none of the known results in the literature regard the theoretical computation of the persistence exponent for FIBM $X_{\alpha }\left(t\right)$, nor have numerical estimates of $\theta \left(\alpha \right)$ been previously obtained. Thus, the aim of the present article is to numerically estimate $\theta \left(\alpha \right)$ as a function of $\alpha \in (0,1)$ by using simulated trajectories of $X_{\alpha }\left(t\right).$ In particular, we are able to validate the following conjecture numerically:

$$\mathrm{for}\alpha \in [0,1]\mathrm{the}\mathrm{persistence}\mathrm{exponent}\mathrm{of}\mathrm{FIBM}{X}_{\alpha }\left(t\right)\mathrm{is}\mathrm{given}\mathrm{by}\theta \left(\alpha \right)=\frac{1}{2(\alpha +1)}.$$

(3)

Note that this formula agrees with the known results in the cases $\alpha =0$ and $\alpha =1,$ and it is also a non-increasing function of $\alpha ,$ according to results in [10]. The idea and some motivations of the conjecture are given in Section 2; the numerical validation of the conjecture will be illustrated in detail in Section 3. Note that our validation is strongly based on numerical simulations of long enough trajectories of the process $X_{\alpha }\left(t\right),$ which require a lot of computation time; so, our analysis could be improved by using a more powerful computer dedicated to the purpose. Unfortunately, we cannot compare our results with those of other authors, since to our knowledge no numerical result of this kind is actually present in the literature.

We emphasize that our study about the persistent exponent of $X_{\alpha }\left(t\right)$ with regard only to the values of $\alpha \in [0,1];$ actually, (3) is only a local conjecture, holding for $\alpha$ belonging to the unit interval. We have not considered extensions of the process $X_{\alpha }\left(t\right)$ to negative values of $\alpha$ or to $\alpha >1,$ nor have we studied the persistent exponent for these values of $\alpha .$ In fact, the conjecture (3) cannot hold for negative $\alpha$ (see [3]) and for $\alpha >1$ (see [11]). The reason we have limited ourselves to study the process $X_{\alpha }\left(t\right)$ and its persistent exponent for $\alpha \in [0,1]$ is due to the fact that we are mainly interested in stochastic processes, such as the FIBM, that model neuronal activities, for which the appropriate range of the fractionally integration parameter $\alpha$ is the unit interval.

The mathematical interest of a such study relies essentially on the need to further investigate the probability laws of $X_{\alpha }\left(t\right)$ and its FPT in order to refine and complete the mathematical setting of the FIBM. Furthermore, this study presents the possibility of shedding light upon a wider class of fractionally integrated stochastic processes and their applications. Indeed, the FIBM has interesting applications in the description of the time evolution of stochastic systems: it appears, e.g., in the framework of certain diffusion models for neuronal activity (see [1], but also [12] for similar models with different processes), where one expects that the inter-spike instants will have a heavy tail distribution, i.e., a power-law decaying rate. The specific choice of fractionally integrated stochastic processes (or, specifically, diffusions) in neuronal modeling allows us to devise models that are more adherent to phenomenological evidence, such those affecting the neuronal spike activity “with memory”for which, after a sequence of short inter-spike times, sequences of long inter-spike times are detected, due to a sort of “adaptation ”([12]).

We can essentially remark that this article is motivated by the aim to study the behavior of the persistent exponent for stochastic processes “with memory” such as the FIBM, by varying the order $\alpha \in [0,1]$ of fractional integration.

1.2. The Results

Our study of the decaying rate of the tail distribution of ${\tau }_S^{\alpha }$ (namely its persistence exponent) is essentially based on numerical simulations. Indeed, by using an ad hoc algorithm implemented in an R-script, we perform simulations of long enough trajectories of the process $X_{\alpha }\left(t\right)$, and the results confirm that for $\alpha =0$ (i.e., case of BM), one has $\theta \left(\alpha \right)=\theta =1/2$, while for $\alpha =1$ (case of integrated BM ), one has $\theta =\theta \left(\alpha \right)=1/4$ (see [6,7]).

For $0<\alpha <1$, our numerical investigation shows that the exponent $\theta \left(\alpha \right)$ decreases as $\alpha$ increases (cf., for instance [3], in which the persistence exponent is revealed to be as a non-increasing function of $\alpha$). We provide numerical estimates and comparisons by means of graphs and tables (see Section 3 for details and discussion of results).

1.3. In Summary

The paper is organized as follows. In the next section, we give the essential ingredients of our study and the main known results from which it starts, and we explain our conjecture, successively supported by simulation results. In Section 3, we describe the specialized algorithm for the simulation paths of the process $X_{\alpha }\left(t\right)$ and of its FPT. We provide graphic results in some figures in order to show and compare the profiles of the FPT density approximations for different values of fractional order $\alpha$. We explain our method to obtain the estimation of the persistence exponent for the simulated cases. Our results are in agreement with the well-known result for the case of $\alpha =1$, and they provide quantitative approximations for cases of $\alpha =0.75,0.5,0.25$, suitably justified under our conjecture. Additionally, we also provide numerical estimates of the probability in (1) in order to show the agreement between the study of the density and of the distribution of the FPT under the conjecture assumption. In Section 4, we highlight some properties of the fractionally integrated processes such as self-similarity that can be useful for implementation purposes of numerical simulations, and long-range dependence that, together with the magnitude of persistence exponent, makes such processes suitable tools for modeling biological dynamics with “memory”. In Section 5, we discuss the possible strategy to be adopted for special Gaussian processes, including Gauss–Markov (GM) processes such as Ornstein–Uhlenbeck (OU). Indeed, thanks to the fact that a GM can be transformed into an OU process (which in turn can be written in terms of BM), one finds that the FPT of a GM process is finite with probability one, and information about the tail behavior of the FPT may be analytically obtained.

2. The Persistence Exponent for Integrated BM: Known Results and a Conjecture

Now, we recall the definition of the fractional integral of order $\alpha \in (0,1)$ of Brownian motion (FIBM):

$$X_{\alpha }\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}B\left(s\right)ds,$$

(4)

where $B_t=B\left(t\right)$ is a standard Brownian motion (BM), and $\mathrm{\Gamma }$ denotes the Gamma Euler function, i.e., $\mathrm{\Gamma }\left(z\right)={\int }_0^{+\infty }{t}^{z-1}{e}^{-t}dt,z>0.$

Taking the limit for $\alpha \to 0^{+},$ one finds that $X_0\left(t\right)$ is BM itself, while for $\alpha =1$, one obtains the ordinary integral of BM. The process $X_{\alpha }\left(t\right)$ starts from zero at $t=0$ with probability 1 (w.p.1) and turns out to be Gaussian with mean zero (cf. [4]); its covariance function as well its variance were studied at length in [1]. Actually, in [1] we have also performed numerical simulations of trajectories of $X_{\alpha }\left(t\right),$ and the probability distribution of ${\tau }_S^{\alpha }$ was numerically studied.

As the case of BM is concerned, it is well-known that its FPT through the barrier $S,$ say ${\tau }_S^B,$ is finite with probability one, though the expectation $E\left({\tau }_S^B\right)=+\infty$ and the exact formula holds (see, e.g., [13]):

$$P({\tau }_S^B\le t)=2(1-\mathrm{\Phi }(S/\sqrt{t}))$$

(5)

where $\mathrm{\Phi }$ denotes the cumulative distribution function of a standard Gaussian variable. Then (see also [14]):

$$P({\tau }_S^B>t)\sim \sqrt{\frac{2}{\pi }}\frac{S}{\sqrt{t}},\mathrm{as}t\to \infty ,$$

(6)

i.e., the persistence exponent is $\theta =1/2$ (see also Example 2.2.2. in [12]). Instead, for non-Markov Gaussian processes such as $X_{\alpha }\left(t\right)$ is, very few results are known about the FPT through a barrier $S.$

Actually, there is an objective difficulty in numerically estimating the FPT distribution using simulated trajectories of the process, since detecting the instant of the first passage through the barrier S can be an arduous task, because the trajectory considered could hit the barrier, but only after a number of simulation time steps which possibly exceed the maximum allowed by the computer algorithm. Therefore, this kind of trajectory is disregarded in the computation of those crossing the boundary within that maximum time of simulation.

Now, we recall the behavior of the tail of the FPT ${\tau }^{\alpha }$ of $X_{\alpha }\left(t\right)$ through the boundary S for well-known cases.

• CASE $\alpha =0$

This is the case when $X\left(t\right):=X_0\left(t\right)=B\left(t\right).$ From Formula (2), it follows that:

$$P({\tau }^B>t)=2\mathrm{\Phi }\left(\frac{S}{\sqrt{t}}\right)-1;$$

then, as easily seen by using the Hospital rule, one gets:

$$\underset{t\to +\infty }{lim}P({\tau }^B>t)\sqrt{t}=\underset{t\to +\infty }{lim}\sqrt{t}\left(2\mathrm{\Phi }\left(\frac{S}{\sqrt{t}}\right)-1\right)=S\sqrt{\frac{2}{\pi }},$$

that is $P({\tau }^B>t)\sim c/\sqrt{t},$ as $t\to +\infty ,$ and the persistence exponent is $\theta =1/2.$

• CASE $\alpha =1$

Now, we have $X\left(t\right):=X_1\left(t\right)={\int }_0^{t}B\left(s\right)ds$ i.e., the ordinary integral of BM. The exact result (see, specifically, [6,15]) is that $\theta =1/4.$ Actually, the estimation of the FPT of integrated BM through the barrier $S,$ numerically obtained by computer simulation, indicates that its probability density behaves as $c·t^{-5/4},$ as $t\to +\infty ,$ that is, the persistence exponent is $\theta =1/4.$ According to this, we find that the tail of the FPT distribution of integrated BM is heavier than that of ${\tau }^B,$ which behaves as $t^{-1/2}.$

The constant c in front of $t^{-5/4}$ was exactly calculated in [15] (see the last formula at pg. 1292), and it is:

$$c=\frac{S^{1/6}{3}^{25/12}\mathrm{\Gamma }(5/12)}{2^{65/12}\pi \sqrt{\pi }}{}_2{F}_1(5/12;7/4;3/2;3/4).$$

(7)

This is the product of the result of Goldman ([6]) and $1/6.$ We also evaluated it by R functions; for the case of $S=1$, we obtained the value $c=0.1795595.$

Remark 1.

One can observe that $X\left(t\right)\sim \mathcal{N}(0,t^3/3),$; so, for fixed $t,X\left(t\right)$ has the same distribution as $\tilde{X}\left(t\right)=B(t^3/3).$ Note that $X\left(t\right)$ is different by $\tilde{X}\left(t\right),$ which is a time-changed BM. Indeed, if $\tilde{\tau }$ denotes the FPT of $\tilde{X}\left(t\right),$ then ${\tilde{\tau }}^3/3=inf\{r>0:B_r>S\}={\tau }^B$, and so:

$$P(\tilde{\tau }>t)=P({\tilde{\tau }}^3/3>t^3/3)=P({\tau }^B>t^3/3)\sim \frac{c}{t^{3/2}},\mathrm{as}t\to +\infty .$$

(8)

This is expected, since the process $\tilde{X}\left(t\right)$ reaches the barrier S more quickly than in the case of BM (being $t^3/3$ much greater than $t,$ for t large).

• CASE $\alpha \in (0,1)$: A conjecture

For the other values of $\alpha \in (0,1),$ our numerical estimations show that the tail of the FPT ${\tau }^{\alpha }$ of the fractionally integrated BM $X_{\alpha }$ through the barrier S is heavier than that of ${\tau }^B$ (corresponding to the case of $\alpha =0$). Precisely, we find that, as $\alpha$ increases in $(0,1]$, the tail becomes heavier and heavier; that is, the persistence exponent $\theta \left(\alpha \right)$ does not increase. Finally, we confirm that the persistence exponent $\theta \left(\alpha \right)$ is a non-increasing function of the fractional order $\alpha$. Indeed, we are confident that our following conjecture holds:

$$\begin{gathered} \mathit{The} \mathit{persistence} \mathit{exponent} \mathit{for} X_{\alpha }\left(t\right) \mathit{is}: \\ \theta =\theta \left(\alpha \right)=\frac{1}{2(\alpha +1)},\alpha \in (0,1). \end{gathered}$$

(9)

Our conjecture is born from the following reasoning. First, we recall the results for the Brownian motion and its FPT ${\tau }^B$ through, e.g., the boundary $S=1$ (recall that the persistent exponent is independent of the boundary), i.e.,

$$P({\tau }^B>t)=P(\underset{0\le s\le t}{max}B\left(s\right)<1)\sim t^{-1/2},\mathrm{for}\mathrm{large}\mathrm{t}$$

(10)

and for the integrated Brownian motion and its FPT $\tau$, i.e.,

$$P(\tau >t)=P\left(\underset{0\le s\le t}{max}{\int }_0^{s}B\left(z\right)dz<1\right)\sim t^{-1/4},\mathrm{for}\mathrm{large}\mathrm{t}.$$

(11)

By recalling the following distribution equality (see, for instance, [3]):

$$P\left(\underset{0\le s\le t}{max}{\int }_0^{s}B\left(z\right)dz<1\right)=P\left(\underset{0\le s\le t}{max}{\int }_0^s(s-z)dB\left(z\right)<1\right)$$

(12)

we also have that

$$P(\tau >t)=P\left(\underset{0\le s\le t}{max}{\int }_0^s(s-z)dB\left(z\right)<1\right)\sim t^{-1/4}\mathrm{for}\mathrm{large}\mathrm{t}.$$

(13)

Note that the last approximation can be interpreted in the two following ways:

$$P\left(\underset{0\le s\le t}{max}{\int }_0^s(s-z)dB\left(z\right)<1\right)\sim {\left(t^{1/2}\right)}^{-1/2}$$

(14)

or

$$P\left(\underset{0\le s\le t}{max}{\int }_0^s(s-z)dB\left(z\right)<1\right)\sim {\left(t^{-1/2}\right)}^{1/2}.$$

(15)

 From this, we do our conjecture for $X_{\alpha }\left(t\right)$. We consider that

$$P(\tau >t)=P\left(\underset{0\le s\le t}{max}{X}_{\alpha }\left(s\right)<1\right)=P\left(\underset{0\le s\le t}{max}{\int }_0^s{(s-z)}^{\alpha -1}B\left(z\right)dz<1\right)$$

(16)

and, by using the following distribution equality

$$P\left(\underset{0\le s\le t}{max}{\int }_0^s{(s-z)}^{\alpha -1}B\left(z\right)dz<1\right)=P\left(\underset{0\le s\le t}{max}\frac{1}{\alpha }{\int }_0^s{(s-z)}^{\alpha }dB\left(z\right)<1\right)$$

we conjecture that, for large $t,$

$$P(\tau >t)=P\left(\underset{0\le s\le t}{max}\frac{1}{\alpha }{\int }_0^s{(s-z)}^{\alpha }dB\left(z\right)<1\right)\sim {\left(t^{\frac{1}{\alpha +1}}\right)}^{-1/2}=t^{-\frac{1}{2(\alpha +1)}}.$$

(17)

The conjecture (17) is equivalent to:

$$P(\tau >t)=P\left(\underset{0\le s\le t}{max}\frac{1}{\alpha }{\int }_0^s{(s-z)}^{\alpha }dB\left(z\right)<1\right)\sim {\left(t^{-\frac{1}{\alpha +1}}\right)}^{1/2}=t^{-\frac{1}{2(\alpha +1)}}.$$

(18)

 In particular, the conjecture, expressed as in (17), can also be explained by means of (2) and by interpreting the persistence exponent $\theta$ for large t as a function of the time t and $\alpha$; i.e., $\theta (t,\alpha )$, such that:

$$\theta (t,0)=\frac{1}{2}=\frac{1}{2·1}\mathrm{in}\mathrm{case}\mathrm{of}{\left(t^{1/(0+1)}\right)}^{-1/2}\sim P\left(\underset{0\le s\le t}{max}{\int }_0^s{(z-s)}^{0}dB\left(z\right)<1\right)$$

$$\theta (t,1)=\frac{1}{4}=\frac{1}{2·2}\mathrm{in}\mathrm{case}\mathrm{of}{\left(t^{1/(1+1)}\right)}^{-1/2}\sim P\left(\underset{0\le s\le t}{max}{\int }_0^s{(z-s)}^{1}dB\left(z\right)<1\right)$$

and consequently,

$$\theta (t,\alpha )=\frac{1}{2(\alpha +1)}\mathrm{in}\mathrm{case}\mathrm{of}{\left(t^{1/(\alpha +1)}\right)}^{-1/2}\sim P\left(\underset{0\le s\le t}{max}\frac{1}{\alpha }{\int }_0^s{(z-s)}^{\alpha }dB\left(z\right)<1\right).$$

 Note that in the conjecture we include the case of $\alpha =0$ (i.e., that of the BM) with persistence exponent $\theta =\frac{1}{2(\alpha +1)}=\frac{1}{2}$, and the case of $\alpha =1$ (i.e., that of the integrated BM) with persistence exponent $\theta =\frac{1}{2(\alpha +1)}=\frac{1}{2(1+1)}=\frac{1}{4}$; all other cases for $\alpha \in (0,1)$ have a persistence exponent $\theta$, such that $\frac{1}{4}\le \theta \le \frac{1}{2}$ with $\theta$ non increasing function of $\alpha .$

 Unfortunately, we are not able to show an analytical reason for the conjecture (17); our heuristic motivation comes by comparing the above equations, and it is confirmed by our numerical computations.

 About the FPT density:

 Actually, by taking the derivative in the expression $P({\tau }^{\alpha }>t)\sim t^{-\theta \left(\alpha \right)},$ the conjecture (9) implies that the FPT density of ${\tau }^{\alpha }$ behaves as $const·t^{-\left(\theta \right(\alpha )+1)},$ as $t\to \infty ,$ where $\theta \left(\alpha \right)$ is the persistent exponent of FIBM $X^{\alpha }\left(t\right).$ Then, inspired by Goldman (see [6,15]), we will suppose that the density of the FPT ${\tau }^{\alpha }$ through the boundary S behaves as

$$c_{\alpha }·S^{1/6}{t}^{-(3+2\alpha )/(2\alpha +2)},$$

(19)

as $t\to \infty ,$ where the multiplicative constant $c_{\alpha }$ is also estimated as a suitable constant multiplied by $\alpha .$ We will validate our conjecture, that is, (19), by means of long trajectories of the process obtained by computer simulation, samples of their FPT and the approximation of the respective densities.

 About the FPT distribution:

 In addition, we will also work with the purpose

(i)

to obtain a numerical estimate of the following probability for t “large enough”

$$P\left(\underset{0\le s\le t}{max}{\int }_0^s{(s-z)}^{\alpha -1}B\left(z\right)dz<1\right)$$

(ii)

to compare it with the function $t^{-\frac{1}{2(\alpha +1)}}$ and from their ratio to derive an estimate of the multiplicative $c_{\alpha }$ such that the asymptotic (in time) tail behavior of the FPT distribution for $S=1$ can be finally characterized as (see Section 3.2 for details)

$$c_{\alpha }·t^{-\frac{1}{2(\alpha +1)}}.$$

 All numerical validations of such a conjecture and approximation results are described in detail in the next section.

3. The Simulation Algorithm and Numerical Results

In the following, we show and compare graphically the behaviour of the tail of the density of FPT ${\tau }^{\alpha }$ for some samples of simulated paths. The simulation algorithm partially follows the R script adopted in [1]. Here, the script is specialized in order to investigate the decaying profile of the FPT densities as the value of the fractional order $\alpha$ varies. The algorithm requires knowledge of the covariance matrix of $X_{\alpha }\left(t\right)$, whose mathematical expression can be found here in Section 4.2. After some statistical comparisons taking into account the execution time of the simulation algorithm, the accuracy of the results and the theoretical expected number of the sample paths crossing the specified boundary, the time discretization step h is chosen as $h=0.05$ here. Moreover, by means of the same sequence of pseudo-random Gaussian numbers, we perform several cases corresponding to different $\alpha$ values.

Specifically, and referring to the methods for generation of pseudo-random numbers (see also [16]), we give the scheme of our R-algorithm:

STEP1.

IN INPUT: we specify the size of the sample of paths we want to simulate, the boundary, the fractional order, the time step, the seed for the random generation routine, and the maximum time length of each path, i.e., the maximum number N of time steps for each path;

STEP2.

we calculate the covariance matrix $C(t_i,t_j)$ with dimension $N\times N$ in $h-$ equi-spaced times $t_i,t_j,i,j=1,\dots ,N$;

STEP3.

we construct the lower triangular matrix $L(i,j)$ such that $C=L{L}^T$ by applying the Cholesky decomposition algorithm to the covariance matrix C;

STEP4.

by the ad hoc R-function (rnorm), for each path, we generate an N-dimensional array z of standard pseudo-Gaussian numbers;

STEP5.

we construct the simulated path with the specified correlation matrix C as the array $x=Lz;$

STEP6.

at each step of the path, we check if crossing of the boundary occurs, and if this is the case, the corresponding FPT is recorded;

STEP7.

we repeat STEPS 4-6 for each path in order to simulate the specified number of paths;

STEP8.

we provide the array of simulated FPTs of dimension equal to the size of the total number of simulated paths (note that if a simulated path does not cross the boundary in the specified (in input) maximum number of steps, a zero value is recorded as its FPT).

STEP9.

IN OUTPUT: the array of FPT is cleaned from zero values and it is used for histograms and kernel density approximations for visualizations and comparisons in figures.

3.1. Graphical and Approximation Results

On the left of Figure 1, we plotted the empirical FPT densities of $X_{\alpha }\left(t\right)$ (on the horizontal axes time $t)$ through the threshold $S=0.1$ for values of $\alpha$ specified in the figure. The empirical density has been obtained by means of 5000 simulated trajectories with a time step equal to $0.05;$ to obtain it, we have counted only the trajectories which have crossed the threshold within the time bound of 300 time steps $h=0.05$, i.e., until time 15. The simulated paths are obtained by applying the algorithm described at the beginning of this section, and their first passage times have been recorded. The plotted curves in the figure are the results of the $plot$ R-function of the $density$ R-function applied to the array of the FPT values of the simulated paths. We underline that the function $density$ computes (Gaussian) kernel density estimates with values for bandwidth “$adjust$ ”parameter between 4 and 12, in such a way that greater values determine smoother curves. The plots of such approximations are provided in the same figure to show how the fractional order $\alpha$ affects the profiles of the FPT densities, and to allow comparisons.

Finally, on the right of Figure 1, we give an enlarged visualization of the tail of the densities on the left of the same figure. We note that as $\alpha$ increases, the tails become heavier, determining smaller values for the corresponding persistence exponents.

First, we show our simulation results for the well-known case of $\alpha =1$ in order to validate our investigation approach. In Figure 2, (left) we show the simulated FPT density (in red) of $X_1\left(t\right)$ through the threshold $S=0.1$. The black curve is the plot of the function we denote as

$$a\left(t\right)=S^{1/6}·t^{-(3+2\alpha )/(2\alpha +2)}=S^{1/6}{t}^{-5/4}$$

with $S=0.1$. Then, we estimated the value of the constant $c_{\alpha }$ as the ratio between the simulated density and the black curve $a\left(t\right)$ at times $\ge 12$ obtaining $c_{\alpha }=0.177.$ Actually, if we denote the simulated density by $d\left(t\right)$, we estimate:

$$c_{\alpha }=\frac{d\left(t\right)}{a\left(t\right)}$$

in values of $t\ge \overline{t}$, $\overline{t}$ being the starting time at which we observe a quite constant difference between $d\left(t\right)$ and $a\left(t\right)$ values. The choice of times “large enough”depends on the value of $\alpha$. In this case, we apply this strategy in order to validate it, because we already know the value of the constant $c_1$, is, as previously specified, $0.1795595$.

Moreover, the blue curve is obtained as the product of the black curve and the constant $c_{\alpha }$, i.e., by denoting as $b\left(t\right)$ the blue curve, we plot, in blue color, the function

$$b\left(t\right)=c_{\alpha }·a\left(t\right)=c_{\alpha }·S^{1/6}·t^{-(3+2\alpha )/(2\alpha +2)}.$$

that, for $\alpha =1$, $S=0.1$ and $c_{\alpha }=0.17$, specifically is

$$b\left(t\right)=0.17·{\left(0.1\right)}^{1/6}·t^{-5/4}.$$

Ultimately, the matching between the two curves (visible in the enlarged visualization on the right of Figure 2) is very good for large enough $t.$ This means that our simulation results are in according to well-known behavior of the FPT tail distribution, i.e., $const·t^{-1/4},$ corresponding to the persistence exponent $\theta \left(1\right)=1/4.$

In Figure 3, Figure 4 and Figure 5, we show the results of the corresponding investigations, as performed to provide the results shown in Figure 2 for the specified cases of $\alpha =0.75,0.5,0.25$, respectively.

In these figures, the simulated FPT densities $d\left(t\right)$ are plotted in red for processes $X_{\alpha }\left(t\right)$ through the same threshold $S=0.1$. The black curves are the plot of the function

$$a\left(t\right)=S^{1/6}·t^{-(3+2\alpha )/(2\alpha +2)}$$

(20)

with the specific value of $\alpha =0.75,0.5,0.25$, respectively. The estimates of $c_{\alpha }$ are obtained as $c_{\alpha }=\frac{d\left(t\right)}{a\left(t\right)}$ for $t\ge {\overline{t}}_{\alpha }$, whose values are specified in

Table 1. For the specified values of fractional order $\alpha$, the persistence exponent $\theta$, the estimates of $c_{\alpha }$, the power decaying function $a\left(t\right)·S^{-1/6}$, the value ${\overline{t}}_{\alpha }$ and the value of the $adjust$ parameter are reported. Note that here, $S=0.1,$ and the value of $c_{\alpha }$ is evaluated as $d\left({\overline{t}}_{\alpha }\right)/a\left({\overline{t}}_{\alpha }\right),$ with $d\left(t\right)$ the simulated density and $a\left(t\right)$ is as in (20). The values of $adjust$ are used in the evaluation of the $density$ R-function to obtain the smoothed simulated density $d\left(t\right)$.

Numerical Approximations for FPT Density
$\mathit{\alpha }$	$\mathit{\theta }$	${\mathit{c}}_{\mathit{\alpha }}$	$\mathit{a}\left(\mathit{t}\right)·{\mathit{S}}^{-\mathbf{1}/\mathbf{6}}$	${\overline{\mathit{t}}}_{\mathit{\alpha }}$	$\mathit{adjust}$
1.0	1/4	0.17	$t^{-5/4}$	14	4
0.75	2/7	0.15	$t^{-9/7}$	13.5	4
0.5	1/3	0.10	$t^{-4/3}$	13.9	8
0.25	$2/5$	0.075	$t^{-7/5}$	14.5	4

. From the presented study, we hypothesize that $c_{\alpha }$ is approximately equal to the product $\alpha ·c_1\approx \alpha ·0.177$, even if we believe further investigations needed to validate this. Finally, our approximation curves, i.e., the blue curves, are obtained as

$$b\left(t\right)=c_{\alpha }·a\left(t\right)=c_{\alpha }·S^{1/6}·t^{-(3+2\alpha )/(2\alpha +2)}$$

(21)

for the values of $\alpha =0.75,0.5,0.25$, respectively.

Furthermore, the values of the persistence exponents for the considered cases are reported in Table 1.

We also recall that our approximation function (the blue curve in Figure 3, Figure 4 and Figure 5) for the specified values of $\alpha =0.75,0.5,0.25$ is finally obtained as $b\left(t\right)$ in (21) by means of values of Table 1, respectively.

Ultimately, for the considered cases, the obtained numerical results, the visible good agreement between the simulated densities (red) and the proposed approximation curves (blue) present in the figures support the proposed conjecture.

3.2. A Further Numerical Approximation

Now, specifically for the tail of the FPT distribution, we aim to give an additional validation of our conjecture, but also a verification tool for previous results inherent to the FPT density. We have already explained the idea of such numerical verification at the end of Section 2. To this end, we realized a specific algorithm in an ad hoc R-script that is implemented with the following steps:

• For a fixed value of $\widehat{t}$ and for a given integer $n,$ we consider the equispaced times $0<t_1<\cdots <t_k<\cdots <t_n=\widehat{t},$ and we generate the pseudo-random values:

  $$E_{t_k}=\left\{{\int }_0^{t_k}{(t_k-z)}^{\alpha -1}B\left(z\right)dz\right\}\mathrm{for}k=1,\dots ,n$$

  by means of the following approximation procedure (with m is specified in input):

  $${\widehat{E}}_{t_k}=\frac{t_k}{m}\sum _{i=1}^{m-1}{(t_k-{\vartheta }_i)}^{\alpha -1}z\left(\sqrt{{\vartheta }_i}\right)$$

  (22)

  where $0<{\theta }_1<\cdots <{\theta }_m=t_k$ are equispaced times and $z\left(\sqrt{{\vartheta }_i}\right),$ with $i=1,\dots ,m-1,$ are pseudo-random numbers generated with Gaussian distribution, having zero mean and $\sqrt{{\vartheta }_i}$ standard deviation.
• We record the value

  $$M_n(\widehat{t})=\underset{k=1,\dots ,n}{max}{\widehat{E}}_{t_k}.$$

  (23)
• We repeat the previous steps for N times in order to obtain a sample of N pseudo-random values

  $$M_n^{\left(j\right)}(\widehat{t})=\underset{k=1,\dots ,n}{max}{\widehat{E}}_{t_k}^{\left(j\right)},j=1,\dots ,N;$$

  we count how many $M_n^{\left(j\right)}(\widehat{t})$ are such that

  $$M_n^{\left(j\right)}(\widehat{t})<1,j=1,\dots ,N,$$

  and we denote such number as ${\mathbb{M}}_n(\widehat{t})$.
• Hence, we approximate the following probability as

  $$P\left(\underset{0\le s\le \widehat{t}}{max}{\int }_0^s{(s-z)}^{\alpha -1}B\left(z\right)dz<1\right)\approx \frac{{\mathbb{M}}_n(\widehat{t})}{N}$$

  and we denote such estimated probability as $p(\widehat{t})=\frac{{\mathbb{M}}_n(\widehat{t})}{N}.$
• Under the conjecture $\theta =-\frac{1}{2(\alpha +1)},$ we compute the value ${\widehat{c}}_{\alpha }$ as follows

  $${\widehat{c}}_{\alpha }=\frac{p(\widehat{t})}{{\widehat{t}}^{-\frac{1}{2(\alpha +1)}}}=\frac{p(\widehat{t})}{{\widehat{t}}^{-\theta }}.$$

  (24)
  We did many executions of such algorithm with different $\widehat{t}$ in order to find the value of $\widehat{t}$ for which we can obtain the agreement between the estimated probability $p(\widehat{t})$ and the expected value, i.e.,

  $$p(\widehat{t})\approx \alpha c_1·{\widehat{t}}^{-\theta }\mathrm{with}{c}_1=0.179.$$
  Indeed, from the previous investigation by means of simulations, we already knew that the approximation $c_{\alpha }\approx \alpha c_1$ can be adopted for the considered cases.

Note that the value of ${\widehat{c}}_{\alpha }$ in (24) depends on the value of $\widehat{t}$. This numerical strategy also provides the value of the time $\widehat{t}$ at which a quite satisfactory agreement of the above approximations is obtained.

The above algorithm has been implemented many times in order to also find the optimal parameters $n,m,N$ to guarantee a sufficient accuracy of the results (we limit ourselves to the second digital digit) and not too long execution times. The results are reported in

Table 2. For the specified values of fractional order $\alpha$, we report the following: the conjectured persistence exponent $\theta$, the number N of trials, the number n for the evaluation of the maximum as in (23), the number m of nodes of the quadrature as in (22), the product $\alpha c_1,$ the estimate ${\widehat{c}}_{\alpha }$ as in (24), the time value $\widehat{t}$ in which the probability is computed, the approximate probability $p(\widehat{t})$ and the expected (conjectured) value $\alpha c_1{\widehat{t}}^{-\theta }$.

Numerical Approximations for FPT Distribution
$\mathit{\alpha }$	$\mathit{\theta }$	$\mathit{N}$	$\mathit{n}$	$\mathit{m}$	$\mathit{\alpha }{\mathit{c}}_{\mathbf{1}}$	${\widehat{\mathit{c}}}_{\mathit{\alpha }}$	$\widehat{\mathit{t}}$	$\mathit{p}(\widehat{\mathit{t}})$	$\mathit{\alpha }{\mathit{c}}_{\mathbf{1}}{\widehat{\mathit{t}}}^{-\mathit{\theta }}$
1.0	1/4	100	75	30	0.179	0.177	15	0.09	0.0909
0.75	2/7	100	70	20	0.13425	0.1275	14	0.06	0.0631
0.5	1/3	500	60	20	0.0895	0.0732	12	0.032	0.0390
0.25	$2/5$	2000	20	20	0.0425	$0.0425$	8	0.0185	0.01849

for specified values of $\alpha$.

We remark that the above probability approximation $p(\widehat{t})$ is a punctual estimation; indeed, it is computed in the time instant $\widehat{t}$ for which we found some results in agreement. We also note that the values of $\widehat{t}$ in Table 2 are different from those of $\overline{t}$ in Table 1, even if these times are those for which we obtain the validation of our conjecture. This depends on the two different numerical strategies used; however, we can say that both times, $\overline{t}$ and $\widehat{t},$ provide information about the starting time of the power decaying of the FPT probability laws (density and distribution, respectively) whose order has been conjectured here.

4. Some Further Properties of the Process ${\mathbf{X}}_{\mathbf{\alpha }}\left(\mathbf{t}\right)$

4.1. On the Self-Similarity Property of $X_{\alpha }\left(t\right)$

Let c be a positive constant; from the definition of the FIBM $X_{\alpha }\left(t\right),$ it follows:

$$X_{\alpha }\left(ct\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^{ct}{(ct-s)}^{\alpha -1}B\left(s\right)ds$$

$$=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^{ct}{\left(ct-c\frac{s}{c}\right)}^{\alpha -1}B\left(s\right)ds;$$

by the variable change $u=s/c,$ and by using the self-similarity property of Brownian motion ($B\left(ct\right)=\sqrt{c}B\left(t\right))$, one finds the following equality in distribution:

$$X_{\alpha }\left(ct\right)=\frac{c^{\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-u)}^{\alpha -1}\sqrt{c}B\left(u\right)du=\frac{c^{\alpha +1/2}}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-u)}^{\alpha -1}B\left(u\right)du=c^{\alpha +1/2}{X}_{\alpha }\left(t\right).$$

Thus, for $c>0$ we obtain the self-similarity property of the FIBM:

$$X_{\alpha }\left(ct\right)=c^{\alpha +1/2}{X}_{\alpha }\left(t\right),$$

or

$$X_{\alpha }\left(t\right)=c^{-(\alpha +1/2)}{X}_{\alpha }\left(ct\right),$$

(25)

where equality is meant in distribution.

From the self-similarity property (25), we obtain:

$$X_{\alpha }\left({\tau }_S\right)=S\Rightarrow c^{-(\alpha +1/2)}{X}_{\alpha }\left(c{\tau }_S\right)=S,$$

and so $X_{\alpha }\left(c{\tau }_S\right)=c^{\alpha +1/2}S.$

Therefore, if $S^{\prime }=c^{\alpha +1/2}S,$ it holds ${\tau }_{S^{\prime }}=c{\tau }_S,$ namely

$${\tau }_S=\frac{1}{c}{\tau }_{S^{\prime }},$$

(26)

where equality is meant in distribution.

Equation (26) provides a relation between the FPT-distribution of the FIBM $X_{\alpha }\left(t\right)$ through the boundary S and the FPT-distribution through the boundary $S^{\prime }=c^{\alpha +1/2}S;$ note that, if $c>1$, one finds that $S^{\prime }>S.$

Now, setting $a=S^{\prime }=c^{\alpha +1/2}S,$ namely

$$c={\left(\frac{a}{S}\right)}^{\frac{1}{\alpha +1/2}},$$

we obtain:

$${\tau }_a=c{\tau }_S={\left(\frac{a}{S}\right)}^{\frac{1}{\alpha +1/2}}{\tau }_S.$$

(27)

Thus, without loss of generality, one can study, e.g., the distribution of ${\tau }_1,$ that is, the FPT of the FIBM through the barrier $S=1$. In fact, from (27), it follows:

$${\tau }_a=a^{\frac{1}{\alpha +1/2}}{\tau }_1.$$

(28)

This fact can be useful during the numerical procedure to obtain the FPT of $X_{\alpha }\left(t\right)$ through the barrier $a.$

4.2. On the Covariance Function of $X_{\alpha }\left(t\right)$

Let $C_{\alpha }(u,t)=cov(X_{\alpha }\left(u\right),X_{\alpha }\left(t\right))$ be the covariance function of $X_{\alpha }\left(t\right),u,t\ge 0;$ we aim to obtain an upper bound for $C_{\alpha }(u,t)$ for $0\le u\le t.$

We recall the following explicit expressions of $C_{\alpha }(u,t)$ and of $Var\left(X_{\alpha }\left(t\right)\right)$ (see Equations (10) and (13) of [1]):

$$C_{\alpha }(u,t)=\frac{1}{{\mathrm{\Gamma }}^2\left(\alpha \right)}\left[\frac{t^{\alpha +1}{u}^{\alpha }}{{\alpha }^2(\alpha +1)}-\frac{t{H}_{\alpha }(u,t)}{\alpha (\alpha +1)}+\frac{J_{\alpha }(u,t)}{\alpha (\alpha +1)}\right],$$

(29)

$$Var\left(X_{\alpha }\left(t\right)\right)=\frac{t^{2\alpha +1}}{(2\alpha +1){\mathrm{\Gamma }}^2(\alpha +1)},$$

(30)

where

$$J_{\alpha }(u,t)={\int }_0^{u}s{(u-s)}^{\alpha -1}{(t-s)}^{\alpha }ds\mathrm{and}{H}_{\alpha }(u,t)={\int }_0^u{(u-s)}^{\alpha -1}{(t-s)}^{\alpha }ds.$$

(31)

Let us suppose that $0\le u\le t$; then, we obtain:

$$0\le H_{\alpha }(u,t)\le {\int }_0^u{(t-s)}^{2\alpha -1}ds;$$

(32)

moreover:

$$0\le J_{\alpha }(u,t)\le {\int }_0^{u}s{(t-s)}^{2\alpha -1}ds\le u{\int }_0^u{(t-s)}^{2\alpha -1}ds$$

$$=\frac{u}{2\alpha }\left[t^{2\alpha }-{(t-u)}^{2\alpha }\right]\le \frac{u{t}^{2\alpha }}{2\alpha }\le \frac{t^{2\alpha +1}}{2\alpha }.$$

(33)

From (29), by using (33) and the fact that $H_{\alpha }(u,t)\ge 0,$ we get:

$$C_{\alpha }(u,t)\le \frac{1}{{\mathrm{\Gamma }}^2\left(\alpha \right)}\left[\frac{t^{2\alpha +1}}{{\alpha }^2(\alpha +1)}+\frac{t^{2\alpha +1}}{2{\alpha }^2(\alpha +1)}\right],$$

namely:

$$C_{\alpha }(u,t)\le \frac{3t^{2\alpha +1}}{2{\mathrm{\Gamma }}^2\left(\alpha \right){\alpha }^2(\alpha +1)},0\le u\le t.$$

(34)

Thus, for fixed u and $t\ge u,$ the covariance function of $X_{\alpha }\left(t\right)$ increases as t increases, but for $t\to +\infty$ it grows at most as a constant × the power $t^{2\alpha +1},$ as the variance does.

5. Some Details on Gauss–Markov Process Fractionally Integrated

Remark 2.

In principle, one could use the previous arguments to study the tail behavior of the FPT of a fractionally integrated Gauss–Markov process $Y\left(t\right).$ We recall that a continuous GM process (see [1]) is in the form:

$$Y\left(t\right)=m\left(t\right)+h_2\left(t\right)B\left(\rho \left(t\right)\right),$$

(35)

where $m\left(t\right),h_1\left(t\right),h_2\left(t\right)$ are continuous functions of $t\ge 0,$ which are $C^1$ in $(0,+\infty )$, such that $h_2\left(t\right)\ne 0,$ and $\rho \left(t\right)=h_1\left(t\right)/h_2\left(t\right)$ is differentiable non-negative, with $\rho \left(0\right)=0$ and ${\rho }^{\prime }\left(t\right)>0$ for $t>0.$ The process $Y\left(t\right)$ has mean $m\left(t\right)$ and covariance $c(s,t)=cov(Y\left(s\right),Y\left(t\right))=h_1\left(s\right)h_2\left(t\right),$ for $0\le s\le t.$ Besides BM itself, a noteworthy case of the GM process is the Ornstein–Uhlenbeck (OU) process; in fact, any continuous GM process can be represented in terms of an OU process (see, e.g., [17]).

Then, the fractional integral of a GM process $Y\left(t\right)$ is defined by:

$$X_{\alpha }\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}Y\left(s\right)ds.$$

(36)

However, the covariance and variance function of $X_{\alpha }\left(t\right)$ are far more complicated than in the case of FIBM (see [1])).

 We conclude with a remark concerning the FPT of OU process.

Remark 3.

Let $Z\left(t\right)$ be OU process, starting from $Z\left(0\right)=z>0,$ driven by the SDE

$$dZ\left(t\right)=-\mu Z\left(t\right)dt+\sigma dB\left(t\right),$$

for some $\mu ,\sigma >0.$ One has

$$Z\left(t\right)=e^{-\mu t}\left(z+B\left(\rho \right(t\left)\right)\right),$$

where

$$\rho \left(t\right)=\frac{{\sigma }^2}{2\mu }\left(e^{2\mu t}-1\right).$$

Set now ${\tau }_Z\left(z\right)=inf\{t>0:Z\left(t\right)=0\},$ then ${\tau }_Z\left(z\right)=inf\{t>0:B\left(\rho \left(t\right)\right)=z\},$ where equality is meant in distribution, so $\rho \left({\tau }_Z\left(z\right)\right)\equiv {\tau }^B\left(z\right)=inf\{s>0:B\left(s\right)=z\}.$ Then, from (5), one gets $P({\tau }_Z\left(z\right)\le t)=2(1-\mathrm{\Phi }(z/\sqrt{\rho \left(t\right)}),$ and ${lim}_{t\to +\infty }P({\tau }_Z\left(z\right)>t)\sqrt{\rho \left(t\right)}=z\sqrt{\frac{2}{\pi }}.$ From this, it follows that the first passage time of $Z\left(t\right)$ through zero has a much lighter tail than that of BM; in particular,

$$\frac{d}{dt}P(\tau \left(z\right)\le t)=\frac{z}{\sqrt{2\pi }{\rho }^{3/2}\left(t\right)}{e}^{-z^2/2\rho \left(t\right)}{\rho }^{\prime }\left(t\right)\sim const·e^{-\mu t},\mathrm{as}t\to +\infty .$$

Therefore, $E\left(\tau \right(z\left)\right)<+\infty ,$ unlike the case of BM, for which $E\left({\tau }^B\right)=+\infty .$

6. Conclusions and Final Remarks

In this paper, we have considered fractionally integrated Brownian motion (FIBM) of order $\alpha \in (0,1]$, that is, $X_{\alpha }\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}B\left(s\right)ds.$ The FIBM is an interesting process, since it appears, e.g., in the framework of diffusion models for neuronal activity (see [1]), where one expects that the inter-spike instants will have a heavy tail distribution, i.e., a power-law decaying rate.

The goal of this paper was to perform a qualitative study of the decaying rate of the tail distribution of ${\tau }_S^{\alpha },$ where ${\tau }_S^{\alpha }=inf\{t>0:X_{\alpha }\left(t\right)\ge S\}$ is the first-passage time (FPT) of $X_{\alpha }\left(t\right)$ through the barrier $S>0.$ Precisely, we have studied the so-called persistent exponent $\theta =\theta \left(\alpha \right)$ of the FPT tail, such that $P({\tau }_S^{\alpha }>t)=t^{-\theta +o\left(1\right)},$ as $t\to +\infty .$ This study has been carried out by numerical simulation of long enough trajectories of the process $X_{\alpha }\left(t\right).$ In fact, we have estimated $\theta \left(\alpha \right),$ as the order $\alpha$ of fractional integration varies in $(0,1]$, and we have showed that it is a non-increasing function of $\alpha \in (0,1],$ with $1/4\le \theta \left(\alpha \right)\le 1/2$. This means that the tail of the distribution of ${\tau }_S^{\alpha }$ becomes heavier and heavier as $\alpha \in (0,1]$ increases. Note that, to our knowledge, none of the known results in the literature regard the theoretical computation of the persistence exponent for FIBM $X_{\alpha }\left(t\right),$ except for $\alpha =0$ (in the case of BM) and $\alpha =1$ (in the case of integrated BM ). Our numerical estimations confirm that for $\alpha =0$, one has $\theta \left(\alpha \right)=1/2,$ while for $\alpha =1$, one has $\theta \left(\alpha \right)=1/4$ (see [6,7]).

In particular, we have numerically validated a new conjecture about the analytical expression of the function $\theta =\theta \left(\alpha \right),\alpha \in (0,1],$ namely $\theta \left(\alpha \right)=\frac{1}{2(\alpha +1)},\alpha \in (0,1].$ Such a numerical validation has been carried out by simulation of long enough trajectories of the process $X_{\alpha }\left(t\right)$ in two ways. In the first one, we have estimated the persistent exponent $\theta \left(\alpha \right)$ by using the simulated FPT density obtained for any $\alpha$. In the second one, we have estimated the persistent exponent by directly calculating $P\left({max}_{0\le s\le t}{\int }_0^{s}B\left(z\right)dz<1\right),$ which is nothing but $P({\tau }_1^{\alpha }>t).$ Both ways confirm our conclusions within the limits of numerical approximation.

In the final part of the paper, we have investigated the self-similarity characteristics of $X_{\alpha }\left(t\right)$ and we have found an upper bound to its covariance function; moreover, we have given some details on the fractionally integrated Gauss–Markov process.

The arguments of this paper allow us, in principle, to also study the decaying rate of the tail distribution (and therefore of the corresponding persistent exponent) of the FPT of the fractional integral of order $\alpha$ of a Gauss–Markov process.