Temporal scaling theory for bursty time series with clusters of arbitrarily many events

Complex systems often show complex dynamical behaviors such as long-term temporal correlations, also known as $1/f$ noise [1, 2, 3, 4, 5]. Characterizing those temporal correlations is of utmost importance to understand mechanisms behind such observations. There are a number of characterization and measurement methods in the literature; e.g., one can refer to Refs. [6, 7, 8, 9] and references therein. One of the most commonly used measurements is an autocorrelation function (ACF) [10, 11]. Precisely, the ACF for a time series $x(t)$ is defined with a time lag $t_{\rm d}$ as

where $\langle\cdot\rangle_{t}$ is the time average over the entire period of the time series. The ACF has been extensively used for detecting temporal correlations in various natural and social phenomena [12, 13, 14, 8]. For the time series with long-term correlations, the ACF typically decays in a power-law form with a decay exponent $\gamma$ such that

This power-law behavior is closely related to the $1/f$ noise via Wiener-Khinchin theorem [15] as well as to the Hurst exponent $H$ [16] and its generalizations [17, 18, 19, 20, 21].

A type of time series that has attracted attention is given in a form of a sequence of event timings or an event sequence, which can be regarded as realizations of point processes in time [22]. Temporal correlations in such event sequences have been characterized by the time interval between two consecutive events, namely, an interevent time (IET). In many empirical data sets, IET distributions, denoted by $P(\tau)$, have heavy tails [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]; in particular, $P(\tau)$ often shows a power-law tail as

with a power-law exponent $\alpha$ [8]. Lowen and Teich derived the analytical solution of the power spectral density for a renewal process governed by a power-law IET distribution [37, 38], concluding that the decay exponent $\gamma$ in Eq. (2) is solely determined by the IET exponent $\alpha$ in Eq. (3) as

It is not surprising because the heavy-tailed IET distribution is the only source of temporal correlations in the time series for renewal processes, where there are no correlations between IETs. The same scaling relation in Eq. (4) was derived in other model studies [39, 40, 41].

In general, correlations between IETs, in addition to the IET distribution, should also be relevant to the understanding of asymptotic decay of the ACF. To detect correlations between consecutive IETs, a notion of bursty trains was introduced [42]; for a given time window $\Delta t$, consecutive events are clustered into a bursty train when any two consecutive events in the train are separated by the IET smaller than or equal to $\Delta t$, while the first (last) event in the train is separated from the last (first) event in the previous (next) train by the IET larger than $\Delta t$. The number of events in each bursty train is called a burst size, and it is denoted by $b$. By analyzing various data, Karsai et al. reported that the burst size distribution shows power-law tails as

with a power-law exponent $\beta$ for a wide range of $\Delta t$ [42]. Similar observations have been made in other data [43, 44, 45, 46, 47, 48]. These findings immediately raise an important question: how does the ACF decay power-law exponent $\gamma$ depend on the IET power-law exponent $\alpha$ as well as on the burst-size power-law exponent $\beta$?

It is not straightforward to devise a model or process that can answer this question, because the heavy tail of the burst size distribution typically implies that arbitrarily many consecutive IETs may be correlated with each other. It is worth noting that correlations between two consecutive IETs have been quantified in terms of the memory coefficient [49], local variation [50], and mutual information [51]; they were implemented using, e.g., a copula method [52, 53] and a correlated Laplace Gillespie algorithm in the context of many-body systems [54, 55]. Correlations between an arbitrary number of consecutive IETs have been modeled by means of, e.g., the two-state Markov chain [42], self-exciting point processes [56], the IET permutation method [57], and a model inspired by the burst-tree decomposition method [48]. Although scaling behaviors of the ACF were studied in some of mentioned works [42, 56, 57, 53], the scaling relation $\gamma(\alpha,\beta)$ has not been clearly understood due to the lack of analytical solutions of the ACF.

In this work, we devise a model for generating a time series with arbitrary functional forms of the IET distribution and burst size distribution. Assuming power-law tails for IET and burst size distributions, our model generates correlations between an arbitrary number of consecutive IETs. By theoretically analyzing the model, we derive asymptotically exact solutions of the ACF from the model time series, enabling us to find the scaling relation $\gamma(\alpha,\beta)$ as follows:

We depict the result in Fig. 1. Note that, for the case of $\beta\gg 1$, the scaling relation in Eq. (4) is partly recovered.

The paper is organized as follows. In Sec. II, we introduce the model with arbitrary functional forms of IET and burst size distributions. In Sec. III, we provide an analytical framework for the derivation of the ACF for the model time series. In Sec. IV, by assuming power-law distributions of IETs and burst sizes, we derive analytical solutions of the ACF, hence the decay exponent $\gamma$ as a function of $\alpha$ and $\beta$. We also compare the obtained analytical results with numerical simulations. Finally, we conclude our paper in Sec. V.

Let us introduce the following model for generating event sequences $\{x(1),\ldots,x(T)\}$ in discrete time, where $T$ is the number of discrete times, for given interevent time (IET) distribution and burst size distribution. By definition, $x(t)=1$ if an event occurs at time $t\in\{1,\ldots,T\}$, and $x(t)=0$ otherwise. By construction, the minimum IET is $\tau_{\rm min}=1$. Furthermore, we only consider bursts with $\Delta t=1$ throughout this work. In other words, we regard that events occurring in consecutive times form a burst, including the case in which the burst contains only one event.

The event sequence $x(t)$ is generated using an IET distribution for $\tau\geq 2$, denoted by $\psi(\tau)$, and a burst size distribution $Q(b)$. Note that $\psi(\tau)$ is normalized. To generate an event sequence, we first randomly draw a burst size $b_{1}$ from $Q(b)$ to set $x(t)=1$ for $t=1,\ldots,b_{1}$. Then, we randomly draw an IET $\tau_{1}$ from $\psi(\tau)$ to set $x(t)=0$ for $t=b_{1}+1,\ldots,b_{1}+\tau_{1}-1$. Note that $\tau_{1}\geq 2$. We draw another burst size $b_{2}$ from $Q(b)$ and another IET $\tau_{2}$ from $\psi(\tau)$, respectively, to set $x(t)=1$ for $t=b_{1}+\tau_{1},\ldots,b_{1}+\tau_{1}+b_{2}-1$ and $x(t)=0$ for $t=b_{1}+\tau_{1}+b_{2},\ldots,b_{1}+\tau_{1}+b_{2}+\tau_{2}-1$. We repeat this procedure until $t$ reaches $T$. See Fig. 2(a) for an example.

We analyze the autocorrelation function (ACF) given by Eq. (1) for the time series $\{x(1),\ldots,x(T)\}$ generated by the model described in Sec. II. Because $A(t_{\rm d}=0)=1$ by definition, we consider the positive time lag ($t_{\rm d}>0$) unless we state otherwise. For an event sequence composed of $n$ events, one gets the event rate as

enabling to write $\langle x(t)^{2}\rangle_{t}=\langle x(t)\rangle_{t}=\lambda$ because $x(t)\in\{0,1\}$. The term $\langle x(t)x(t+t_{\rm d})\rangle_{t}$ in Eq. (1) is written as follows:

where $Z(t_{\rm d})$ is the probability that $x(t+t_{\rm d})=1$ conditioned on $x(t)=1$. Note that $Z(0)=1$. Then, the ACF reads

Let us consider the case in which $x(t)x(t+t_{\rm d})$ is non-zero, i.e., $x(t)=1$ and $x(t+t_{\rm d})=1$. As depicted in Fig. 2(a), the time series $x(t)$ in a period of $[t,t+t_{\rm d}]$ is typically composed of several alternating bursts and IETs larger than one. Here the consecutive IETs of $\tau=1$ forms a burst and the sum of such IETs of $\tau=1$ is equal to the burst size minus one. Therefore, the time lag $t_{\rm d}$ is written as a sum of burst sizes (each minus one) of bursts appeared in the period of $[t,t+t_{\rm d}]$ and IETs between those bursts. We note that the first burst contains an event at time $t$ and that the last burst contains an event at time $t+t_{\rm d}$.

We denote by $r$ the number of events from time $t$ to the end of the burst containing the event at time $t$. By definition, we exclude the event at time $t$ when counting $r$, hence $r\geq 0$. For example, if a burst is composed of $5$ events occurring at times $1,\ldots,5$ and we consider $t=2$, then one obtains $r=3$. If we select $t$ such that $x(t)=1$ uniformly at random, $t$ is contained in a burst of size $b$ with a probability proportional to $Q(b)$. Therefore, $r$ is a discrete-time variant of the waiting time or the residual time derived from the IET [8], and the probability distribution of $r$ is given by

Note that $\sum_{r=0}^{\infty}R(r)=1$. See Appendix A for the derivation of the factor $1/\langle b\rangle$.

Next we denote by $r^{\prime}$ the number of events that are in the burst containing the event at time $t+t_{\rm d}$ and occur before or at $t+t_{\rm d}$ [Fig. 2(b)]. It should be noted that $r^{\prime}\geq 1$. The definition of $r^{\prime}$ implies that the first event of the burst containing the event at $t+t_{\rm d}$ occurred at time $t+t_{\rm d}-r^{\prime}+1$. We denote by $q(r^{\prime})$ the probability that $x(t+t_{\rm d})=1$ conditioned that the first event of the burst containing the event at time $t+t_{\rm d}$ is located at $t+t_{\rm d}-r^{\prime}+1$. For this case to occur, the size of the burst starting at time $t+t_{\rm d}-r^{\prime}+1$ has to be larger than or equal to $r^{\prime}$. Therefore, one obtains

To derive $Z(t_{\rm d})$ in Eq. (16), we denote by $p_{k}(t_{\rm d})$ the probability that an event occurs at time $t+t_{\rm d}$ and there are exactly $k$ alterations between the burst and the IET larger than one, conditioned that an event occurs at time $t$. By $k$ alterations, we mean that the event at time $t+t_{\rm d}$ belongs to the $k$th burst after the burst to which the event at time $t$ belongs to. Note that there are then $k$ IETs larger than one between the burst at time $t$ and that at time $t+t_{\rm d}$. Then $Z(t_{\rm d})$ can be written in terms of $p_{k}(t_{\rm d})$ as

For the case with $k=0$, we obtain using Eq. (17)

It should be noted that, when counting $t_{\rm d}$ for $k=0$, we exclude the event at time $t$ and include the event at time $t+t_{\rm d}$ for consistency. For example, consider a burst of $5$ events at times $1,\ldots,5$. If we are considering the two events at times 2 and 4, we set $t=2$ and $t_{\rm d}=2$.

If there are $k$ $(\geq 1)$ IETs larger than one intersecting $[t,t+t_{\rm d}]$, one can write [Fig. 2(b)]:

where we have defined a reduced IET by

for convenience. Note that $\bar{\tau}_{i}\geq 1$ for all $i$ because $\tau_{i}\geq 2$. We assume that all variables on the right-hand side of Eq. (21) are statistically independent of each other. Then, we can write $p_{k}(t_{\rm d})$ for $k\geq 1$ as

It is also remarkable that using the reduced IET, the event rate is obtained as follows:

where

Namely, $\lambda$ is the ratio of the average length of the period of $x(t)=1$ to the sum of those of $x(t)=1$ and of $x(t)=0$ [Fig. 2(b)].

For analytical tractability, we assume that all variables on the right-hand side of Eq. (21) are real numbers. Therefore, $\psi(\bar{\tau})$, $Q(b)$, $R(r)$, $q(r^{\prime})$, and $p_{0}(t_{\rm d})$ are also considered for their respective continuous variables. The continuous versions of Eqs. (17) and (18) are given by

and

respectively. Furthermore, the continuous-time versions of Eqs. (20) and (23) respectively read

and

where $\delta(\cdot)$ is the Dirac delta function.

We take the Laplace transform of Eq. (28) to obtain

where $\tilde{Q}(s)$ denotes the Laplace transform of $Q(b)$. The Laplace transform of Eq. (29) reads for $k\geq 1$

where

and $\tilde{\psi}(s)$ denotes the Laplace transform of $\psi(\bar{\tau})$. Then the Laplace transform of $Z(t_{\rm d})$ in Eq. (19) is obtained as

By taking the inverse Laplace transform of $\tilde{Z}(s)$ and then substituting it into Eq. (16), one can obtain the analytical solution of the ACF for $x(t)$.

To demonstrate the above analytical results, we consider a simple case in which both reduced IETs and burst sizes are real numbers and exponentially distributed. By assuming that

and

we derive an exact solution of the ACF as follows (see Appendix B):

Let us return to our original question on temporal scaling behavior. We assume continuous versions of power-law distributions of reduced IETs and burst sizes as follows:

Here $\alpha,\beta>1$ are power-law exponents, and $\bar{\tau}_{\rm c}$ and $b_{\rm c}$ are cutoffs. Also, $C_{\psi}\equiv\frac{\alpha-1}{1-\bar{\tau}_{\rm c}^{1-\alpha}}$ and $C_{Q}\equiv\frac{\beta-1}{1-b_{\rm c}^{1-\beta}}$ are normalization constants. Then we will derive the analytical result of the decay exponent $\gamma$ of the ACF as a function of $\alpha$ and $\beta$, i.e., $\gamma(\alpha,\beta)$.

We first prove a useful property that $\gamma$ is symmetric with respect to the exchange of $\alpha$ and $\beta$, namely,

To prove this property, let us consider a complementary event sequence $y(t)$ to the original event sequence $x(t)$, which is defined as

The ACF defined for $y(t)$ using the formula in Eq. (1) turns out to be the same as the ACF for $x(t)$:

That is, the decay exponent of the ACF for $x(t)$ must be the same as that for $y(t)$. By the definition of $y(t)$ in Eq. (41), the periods of $x(t)=0$ correspond to those of $y(t)=1$ and vice versa. It means that reduced IETs and burst sizes in $x(t)$ respectively correspond to burst sizes and reduced IETs in $y(t)$, closing the proof.

Under Eq. (38), the average of $\bar{\tau}$ is given by

Similarly, under Eq. (39), the average of $b$ reads

Note that the event rate $\lambda$ in Eq. (24) is determined by the above $\langle\bar{\tau}\rangle$ and $\langle b\rangle$.

We now divide the entire range of $\alpha,\beta>1$ into several cases to derive the analytical solution of the ACF in each case. Considering the symmetric nature of $\gamma$ in Eq. (40), the following cases are sufficient to get the complete picture of the result.