# Renewal theory with fat tailed distributed sojourn times: typical versus rare

###### Abstract

Renewal processes with heavy-tailed power law distributed sojourn times are commonly encountered in physical modelling and so typical fluctuations of observables of interest have been investigated in detail. To describe rare events the rate function approach from large deviation theory does not hold and new tools must be considered. Here we investigate the large deviations of the number of renewals, the forward and backward recurrence time, the occupation time, and the time interval straddling the observation time. We show how non-normalized densities describe these rare fluctuations, and how moments of certain observables are obtained from these limiting laws. Numerical simulations illustrate our results showing the deviations from arcsine, Dynkin, Darling-Kac, Lévy and Lamperti laws.

###### pacs:

02. 50. -r, 05. 20. -y, 05. 40. -a^{†}

^{†}preprint: APS/123-QED\definechangesauthor

[name=Author, color=red]Author

## I Introduction

Renewal processes Godreche2001Statistics (); Mainardi2004fractional (); Mainardi2007Beyond (); Godreche2015Statistics (); Niemann2016Renewal (); Akimoto2016Langevin () are simple stochastic models for events that occur on the time axis when the time intervals between events are independent and identically distributed (IID) random variables. This idealised approach has many applications, ranging from the analysis of photon arrival times to queuing theory. In some models the sojourn time probability density function (PDF) has fat tails, and this leads to fractal time renewal processes. In the case when the variance of the sojourn time diverges, we have deviations from the normal central limit theorem and/or the law of large numbers. Such fat tailed processes are observed in many systems, ranging from blinking quantum dots Simone2005Fluorescence (), to diffusion of particles in polymer networks Edery2018Surfactant (), or diffusion of particles on the membrane of cell weron2017Ergodicity () to name a few. In these systems the renewal process is triggering jumps in intensity or in space. The continuous time random walk model Metzler2000random (), the annealed trap model, the zero crossing of Brownian motion, the velocity zero crossing of cold atoms diffusing in momentum space Barkai2014From (), are all well known models which use this popular renewal approach (see however Boettcher2018Aging (); Nyber2018Zero ()). Heavy tailed renewal theory is also used in the context of localization in random waveguides Fernandez2014Beyond (). The number of renewals, under certain conditions, is described by Lévy statistics, and the fluctuations in these processes are large. Hence it is important to explore the rare events or the far tails of the distributions of observables of interest. As mentioned in Touchette2009large (); Whitelam2018Large () the large deviation principle, with its characteristic exponential decay of large fluctuations, does not describe this case, and instead the big jump principle Alessandro2018Single () is used to evaluate the rare events in Lévy type of processes.

The main statistical tool describing observables of interest are non-normalised states, which are limiting laws with which we may obtain statistical information on the system, including for example the variance, which in usual circumstances is the way we measure fluctuations. These non-normalised states were previously investigated, in the context of Lévy walks Rebenshtok2014Infinite (), spatial diffusion of cold atoms Erez2017Large (), and very recently for Boltzmann-Gibbs states when the underlying partition function of the system diverges Erez2018From (). These functions describing the statistical behavior of the system are sometimes called infinite densities or infinite covariant densities, and they appear constantly in infinite ergodic theory Aaronson1997introduction ().

Our goal in this paper is to investigate the statistics of rare events in renewal theory. Consider for example a non-biased ordinary random walk on the integers. The spatial jump process is Markovian hence the zero crossing, where the zero is the origin, is a renewal process. Here like Brownian motion, the waiting time PDF between the zero crossings is fat tailed, in such a way that the mean return time diverges. The distribution of the occupation time , namely the time the random walker spends in the positive domain is well investigated Redner2001guide (); Godreche2001Statistics (). Naively one would expect that when the measurement time is long the particle will spend half of its time to the right of the origin. Instead one finds that this is the least likely scenario, and the PDF of the properly scaled occupation time reads

(1) |

Here and all along this manuscript the subscript denotes the observable of interest, e.g. we consider the PDF of which attains values . This arcsine law, which describes also other features of Brownian motion Morters2010Brownian (); Akimoto2016Distributional (); Sadhu2018Generalized (), exhibits divergences on or . Here a particular scaling of is considered. However, in cases studied below we show that other limiting laws are found when a second time scale is considered and these may modify the statistical properties of the occupation time when is either very small or very large. This in turn influences the anticipated blow up of the arcsine law at its extremes. Notice that here the least likely event, at least according to this law is the case , so our theory is not dealing with corrections to the least likely event, but rather corrections to the most likely events. This is because of the fat tailed waiting times, which make the discussion of deviations from familiar limiting laws a case study in its own right. While the theory deals with most likely events, from the sampling point of view these are still rare, as the probability of finding the occupation time in a small interval close to the extremes of the arcsine law is still small.

The organization of the paper is as follows. In section II, we outline the model and give the necessary definitions. The behavior of the probability of observing renewals in the interval , is analyzed in section III. In sections IV, V and VI, the densities of the forward and backward recurrence time, and the time interval straddling , denoted , and respectively, are derived. In order to see the effects of the typical fluctuations and large deviations, the fractional moments, e.g., , are considered and bi-fractal behavior is found. In section VII, the behavior of the occupation time is studied. In the final section, we conclude the paper with some discussions. All along our work we demonstrate our results with numerical experiments and compare between the statistical laws describing typical fluctuations to those found here for the rare events.

## Ii Model

Renewal process, an idealized stochastic model for events that occur randomly in time, has a very rich and interesting mathematical structure and can be used as a foundation for building more realistic models Metzler2004The (); Brokmann2003Statistic (). As mentioned, the basic mathematical assumption is that the time between the events are IID random variables. Moreover, renewal processes are often found embedded in other stochastic processes, most notably Markov chains.

Now, we briefly outline the main ingredients of the renewal process Godreche2001Statistics (). It is defined as follows: events occur at the random epochs of time , , , , , from some time origin . When the time intervals between events, , , , , , are IID random variables with a common PDF , the process thus formed is a renewal process (see the top panel of the Fig. 1). We further consider the alternating renewal process in which the process alternates between and states. A classical example is a Brownian motion in dimension one, where we denote state with and state for . Generically, we imagine that a device, over time, alternates between on and off states, like a blinking dot Simone2005Fluorescence (); Gennady2005Power (). Here we suppose the process starts in state and stays in that state for a period of time , then goes to state and remains for time ; see bottom panel of the Fig. 1. Clearly, it is natural to discuss the total time in state or . and are called the occupation times in the and state, respectively and . For Brownian motion, and the distribution of time in state is the well known arcsine law.

Motivated by previous studies of complex systems, we consider here PDFs with power law tails, i.e., for large

In this case, the first moment of is divergent for . Here the index and is a time scale. As we show below, the full form of is of importance for the study of the large fluctuations. An example is the fat tailed PDF Metzler2000random ()

(2) |

Using the Tauberian theorem Feller1971introduction (), in Laplace space

(3) |

for small , where is conjugate to , , and . In order to simplify the expression, we denote as the Laplace transform of . When , the first moment is finite and the corresponding Laplace form Metzler2000random () is

(4) |

for small . Notice that , which means that the PDF is normalized. We would like to further introduce the one sided Lévy distribution with index , which is used in our simulations to generate the process; see Appendix A. In Laplace space, one sided stable Lévy distribution is Metzler2004The ()

(5) |

and the small expansion is given by with . For specific choices of , the closed form of the is tabulated for example in MATHEMATICA Burov2012Weak (). In particular, a useful special case is

(6) |

It implies that for large , so the first moment of the sojourn time diverges.

In the following we will draw on the research literature given by Godrèche and Luck Godreche2001Statistics (), which is recommended for an introduction. The number of renewal events in the time interval between and time is

(7) |

Then we have the following relation . Now we introduce the forward recurrence time , the time between and the next event

see Fig. 1. While the corresponding backward recurrence time, the length between the last event before and the observation time , is defined by

Utilizing the above two equations, we get the time interval straddling time , i.e., , which is

For simplification, we drop the subscript, denoting the time dependence of the random quantities, from here on.

## Iii Number of renewals between and

We recap some of the basic results on the statistics of the number of renewal events. The probability of the number of events up to time is

(8) |

and is the probability to have an event at time , defined by

(9) |

Here is the survival probability

(10) |

that is, the probability that the waiting time exceeds the observation time . For power law time statistics and large ,

Using Eqs. (8, 9) and convolution theorem leads to Montroll1965Random ()

(11) |

with .

### iii.1 Number of renewals between and with

Rewriting Eq. (11), using the convolution theorem of Laplace transform, and performing the inverse Laplace transform with respect to , we get a formal solution

(12) |

where is a discrete random variable and means the inverse Laplace transform, from the Laplace space to real space .

Summing the infinite series (summation over ), the normalization condition is discovered as expected. We notice that Eq. (12) can be further simplified when is one sided Lévy distribution Eq. (5). Then the inverse Laplace transform of Eq. (12) gives

(13) |

As usual the large time limit is investigated with the small behavior of . Utilizing Eq. (5), the behavior of Eq. (11) in the large limit and small is,

Here note that . This means that with this approximation is treated as a continuous variables, which is fine since in fact we consider a long time limit, and the limiting PDF of is approaching a smooth function. Hence, we have to denote the continuous approximation. First, using the property of Laplace transform, i.e., , secondly, performing the inverse Laplace transform on the above equation, we find the well known result Aaronson1997introduction (); Godreche2001Statistics (); Schulz2014Aging ()

(14) |

Eq. (14) is customarily called the inverse Lévy PDF. Furthermore, using , we find that can be expressed as Mittag-Leffler probability density

where is the Laplace transform of with respect to and a two-parameter function of the Mittag-Leffler type is defined by the series expansion Podlubny1999Fractional ()

with and . Eq. (14) describes statistics of functionals of certain Markovian processes, according to the Darling-Kac theorem. It was also investigated in the context of infinite ergodic theory Aaronson1997introduction () and continuous time random walks.

The well known limit theorem Eq. (14) is valid when and are large and the ratio is kept fixed. Now we consider rare events when is kept fixed and finite, say , , , and is large. Using Eq. (12) we find

(15) |

Note that is a probability, while is a PDF. To make a comparison between Eq. (14) and Eq. (15) we plot in Fig. 2, the probability that is in the interval versus and compare these theoretical predictions to numerical simulations. Integrating Eq. (14) between and , gives what we call the typical fluctuations. While the result Eq. (15) exhibits a staircase since according to this approximation

(16) |

where gives the greatest integer less than or equal to . From Fig. 2 we see that, besides the obvious discreteness of the probability, deviations between the two results can be considered marginal and non-interesting. Luckily this will change in all the examples considered below, as the statistical description of rare events deviates considerably from the known limit theorems of the field.

### iii.2 Number of renewals between and with

Based on Eq. (11), we obtain a useful expression

(17) |

Here we consider the random variable, , and explore its PDF denoted . Applying Fourier-Laplace transform, and , the PDF of in Fourier-Laplace space is

(18) |

First, we consider the limit of small and small , and the ratio is fixed. As we discuss below this leads to the description of what we call bulk or typical fluctuations, and these are described by standard central limit theorem. Substituting into the above equation and taking inverse Laplace transform

(19) |

Fourier inversion of the above equation yields the PDF , written in a scaling form Godreche2001Statistics ()

(20) |

with and ; see Fig. 3. We see that for fixed observation time the parameter measures the PDF’s width. Furthermore, the function is defined by

where is the asymmetric Lévy PDF; see Appendix C. Compared with the one sided Lévy distribution, holds two sides with the right hand side decaying rapidly. Moreover, the second moment of diverges for .

As well known the central limit theorem (here of the Lévy form) describes the central part of the distribution, but for finite though large it does not describe the rare events, i.e., the far tail of the distribution. So far we investigated the typical or bulk statistics and as we showed they are found for . Technically this was obtained using the exact Laplace-Fourier transform, and then searching for a limit where and are small their ratio finite, as mentioned. However, it turns out that this limit is not unique. As we now show we can use the exact solution, assume both and are small, but their ratio finite and obtain a second meaningful solution. This in turn, leads to the description of rare events, i.e., the far tail of the distribution of the random variable . Roughly speaking, in this problem (and similarly all along the paper) we have two scales, one was just obtained and it grows like , the second (with this example) is , as we now show. This means that we have two ways to scale data, one emphasizing the bulk fluctuations (explained already) and the second the rare events.

We now consider a second limiting law capturing the rare events valid when is of the order of . From Eq. (18), we have

Keep in mind that and are small and they are the same order. After performing inverse Fourier-Laplace transform, the asymptotic behavior of is

(21) |

This is the main result of this section. Here the above equation is only valid for negative . We see that decays like for small negative . Moreover, the scaling behavior of yields

(22) |

with ; see Fig. 4. It means that decays like for , thus is not normalized. Furthermore, for , the dominating term matches the left tail of Eq. (20); see Eq. (92) in Appendix C.

For fixed observation time , the central part of the PDF is well illustrated by the typical fluctuations Eq. (20). While, its tail is described by Eq. (22), exhibiting the rare fluctuations. In order to discuss the effect of typical fluctuations and large deviations, we further consider the absolute moment of Schulz2015Fluctuations (), defined by

(23) |

(24) |

Here we use the fact that is a finite constant for . Note that to derive Eq. (24) we use the non-normalized solution Eq. (21) for , indicating that Eq. (21) while not being a probability density, does describe the high order moments. In the particular case (high order moment), we have . While this result is known Godreche2001Statistics (), our work shows that the second moment, in fact any moment of order , stems from the non-normalized density describing the rare fluctuations Eq. (21). Other examples of such infinite densities will follow.

###### Remark 1

From simulations of the number of renewals , Fig. 4, we see deviations from typical results when . As mentioned, our theory covers the case , so there is a need to extend the theory further. Note that for , the typical fluctuations decay rapidly, while power law decay, for intermediate values of , is found on the left (see Fig. 4). This intermediate power law behaviors can not continue forever, since , and hence when a new law emerges, Eq. (22). Possibly the large deviation principle can be used to investigate the case .

## Iv The forward recurrence time

Several authors investigated the distribution of both for , meaning is of the order of , for and also for ; see Refs. Dynkin1961Selected (); Feller1971introduction (); Schulz2014Aging (). These works considered the typical fluctuations of , while we focus on the events of large deviations. This means that we consider for and for . The forward recurrence time is an important topic of many stochastic processes, such as aging continuous time random walk processes (ACTRW) Schulz2014Aging (); Kutner2017continuous (), sign renewals of Kardar-Parisi-Zhang Fluctuations Takeuchi2016Characteristic () and so on. The forward recurrence time, also called the excess time (see schematic Fig. 1), is the time interval between next renewal event and . In ACTRW, we are interested in the time interval that the particle has to wait before next jump if the observation is made at time t. The PDF of the forward recurrence time is related to according to

(25) |

see Eq. (85) in Appendix B. In double Laplace space, the PDF of Godreche2001Statistics () is

(26) |

Based on the above equation, we will consider its analytic forms and asymptotic ones. In general case, the inversion of Eq. (26) is a function that depends on and . While, for , the above equation can be simplified as , which is independent of the observation time . As expected, for this example we do not have an infinite density, neither multi-scaling of moments, since and more generally thin tailed PDFs, do not have large fluctuations like Lévy statistics.

### iv.1 The forward recurrence time with

First, we are interested in the case of . In Laplace space, this corresponds to . From Eq. (26)

(27) |

We notice that Eq. (27) can be further simplified for a specific , namely Mittag-Leffler PDF Podlubny1999Fractional (); Kozubowski2001Fractional (). In order to do so, we consider

(28) |

with . In Laplace space, has the specific form

(29) |

This distribution can be considered as the positive counterpart of Pakes’s generalized Linnik distribution Jose2010Generalized () with the PDF having the form , , . Plugging Eq. (29) into Eq. (27) leads to

Taking the the double inverse Laplace transform yields

(30) |

see Fig. 6. Notice that , so for , the PDF of for gives .

More generally, using Eq. (3), we have

Performing inverse double Laplace transform leads to the main result of this section, and the density describing the large deviations is

(31) |

which exhibits interesting aging effects Schulz2014Aging (). Here for large is equal to , namely is increasing with measurement time , and for reasons that become clear later we may call it the effective average waiting time (recall the is a constant only if ). The large deviations shows that for large the forward recurrence time decays as . Furthermore, the integration of Eq. (31) over diverges since is not integrable for large . Hence Eq. (31) is not a normalised density. For that reason, we may call in Eq. (31) an infinite density Rebenshtok2014Infinite (), the term infinite means non-normalizable, hence this is certainly not a probability density. Even though Eq. (31) is not normalized, it is used to obtain certain observables, such as averages of observables integrable with respect to this non-normalized state. Besides, infinite densities play an important role in infinite ergodic theory Thaler2006Distributional (); Akimoto2012Distributional () and intermittent maps Korabel2009Pasin ().

Using Eq. (26), we find a formal solution to the problem

(32) |

where ‘’ is the Laplace convolution operator with respect to and the double Laplace transform of the function is

We further discuss a special choice of , i.e., . After some simple calculations, Eq. (32) gives

(33) |

For Mittag-Leffler waiting time Eq. (28), we obtain

(34) |

from which we get the PDF of plotted in Fig. 6.

We now focus on the typical fluctuations, namely the case and both are large. This means that and are small but of the same order. Plugging Eq. (3) into Eq. (26), then taking double inverse Laplace transform, leading to the normalized solution Godreche2001Statistics (); Dynkin1961Selected ()

(35) |

which is plotted by the dashed (black) lines in Figs. 6 and 6. The well known solution Eq. (35) describes the typical fluctuations when .

To summarize, the forward recurrence time shows three distinct behaviors: for , the infinite density Eq. (31) rules, and only in this range, the PDF of depends on the behavior of ; for , both Eqs. (31) and (35) are valid and predict ; for , we use Eq. (35) and then . Note that for certain observables, for example and , when , their PDFs are also governed by the shape of ; see below.

### iv.2 The forward recurrence time with

For , according to Eq. (26)

where as mentioned is finite. This can be finally inverted, yielding the typical fluctuations Feller1971introduction (); Tunaley1974Asymptotic (); Godreche2001Statistics ()

(36) |

Since , Eq. (36) is a normalized PDF and independent of the observation time , which is different from Eq. (31), but they have similar forms. This is the reason why in the previous section we called the effective average waiting time.

Next we discuss the uniform approximation, which is valid for varieties of and large , namely within uniform approximation, we have the only condition that is large but the ratio of and arbitrary. It can be noticed that Eq. (26) can be arranged into the following formula

For , we may neglect the second term, then using and inverting we get

(37) |

which captures both the infinite density and the bulk fluctuations; see Fig. 7. Here, Eq. (37) is true for large without considering the relation between and . If , Eq. (37) can be approximated by Eq. (36).

For the rare fluctuations, i.e., both and are small and comparable, inserting Eq. (4) into Eq. (26), yields

For , taking the double inverse Laplace transform, we find

(38) |

which is consistent with Eq. (36) for large and . Besides, for , decays as independent of the observation time . On the other hand, if , grows linearly with , namely . In addition, using the asymptotic behavior of , for large the uniform approximation Eq. (37) reduces to Eq. (38). Still as for other examples in this manuscript, we may use Eq. (38) to calculate a class of high order moments (for example ), i.e., those moments which are integrable with respect to this infinite density.

###### Remark 2

For simulations presented in Fig. 7, we use particles on a standard workstation, taking about day. We see that in this case we do not sample the rare events. In Ref. Rebenshtok2014Infinite (), simulations of the Lévy walk process with particles is performed, in order to explore graphically the far tails of the propagator of the Lévy walk. When increasing the number of particles, we will observe rare events, however clearly in our case realizations are simply not sufficient for meaningful sampling.

With the help of the above equations, now we turn our attention to the fractional moments, defined by

(39) |

Utilizing Eq. (39) and integration by parts, yields

(40) |

This is to say, for , is a constant, namely, it does not depend on the observation time . Moments of order are determined by the known result Eq. (36), which describes typical fluctuations when is of the order of . The rare fluctuations, described by Eq. (38), give information to events with , and this non-normalized density Eq. (38) yields the moments of ; see Eq. (94) in Appendix D. Especially, if , so in this case the mean is determined by the infinite density. When , is divergent. This is expected since the moment of order of diverges.

## V The backward recurrence time

Compared with the forward recurrence time, one of the important difference is that can not be larger than . In some cases, is called the age at time . Because in the lightbulb lifetime example, it represents the age of the light bulb you find burning at time t. Similar to the derivation of the forward recurrence time

In Laplace space, let and . Using the convolution theorem of Laplace transform and Eq. (84), this gives

(41) |

which was derived in Ref. Godreche2001Statistics () using a different method.

### v.1 The backward recurrence time with

First of all, we study the behaviors of large deviations. For , i.e.,

(42) |

In the long time limit, i.e., , . Performing the double inverse Laplace transform with respect to and , respectively, yields

(43) |

where is defined below Eq. (31). It implies that , which is confirmed in Fig. 9. Furthermore, note that , which means that Eq. (43) is non-normalized.

Now we construct a uniform approximation, which is valid for a wider range of , though is large. We rewrite Eq. (41) as

For simplification, let be the one sided Lévy stable distribution Eq. (5). Expanding the above equation, i.e., , and then using the convolution theorem of the Laplace transform

(44) |

where represents the Heaviside theta function Polyanin2007Handbook (), which is equal to for and for . The in Eq. (44) yields as expected. In addition, for , Eq. (44) reduces to Eq. (43). Note that, for , and comparable and , Eq. (44) is consistent with the arcsine law, while, let go to either or (the extreme cases), the arcsine law does not work anymore; see Fig. 9.

Now we turn our attention to the case of , using the random variable . In Laplace space, the PDF of is

(45) |

According to Eq. (41)

(46) |

For , performing the double inverse Laplace transform and using

(47) |

Let us consider a situation in which is the Mittag-Leffler distribution Eq. (28) with . Next, plugging Eq. (29) into Eq. (47) yields

(48) |

It demonstrates that decays like . Thus, if tends to the observation time , we discover an interesting phenomenon that , verified in Fig. 10. In general case, Eq. (47) is not easy to calculate in real time exactly, though we use the numerical inversion of Laplace transform by MATLAB. Expanding the above equation, we find

Consider a specific , namely one sided Lévy stable distribution

(49) |

which can be used for plotting. To summarize, large deviations are observed for and , Eq. (43) and Eqs. (48, 49) respectively (see Fig. 9), and these are non-symmetric for one sided Lévy distribution. Only when , we find a non-normalized density, Eq. (43).

Next we discuss the typical fluctuations when . Combining Eqs. (3) and (41), yields Godreche2001Statistics ()

(50) |

In a particular case , Eq. (50) reduces to the arcsine law , which is plotted by the dashed (black) line in Figs. 9, 9 and 10. Let , we get a well known formula Godreche2001Statistics ()