Modelling bursty time series

Studying human activity patterns is of central interest due to the wide practical usage. Understanding the dynamics underlying the timing of various human activities—such as starting a phone call, sending an e-mail or visiting a website—are crucial to modelling the spreading of information or viruses [1]. Modelling human dynamics is also important in resource allocation. It has been shown that for many human activities the interevent time distribution follows a power law with exponents usually varying between 1 and 2. Processes with power-law decaying interevent time distribution look very different from the Poisson process, which has been used to model the timing of human activities before [2, 3]. While time series from the latter look rather homogeneous, the former processes produce bursts of rapidly occurring events, which are separated by long inactive periods.

Some examples where power-law decaying interevent time distribution has been observed are email-communication (with exponent β ≈ 1, [4]), surface mail communication (β ≈ 1.5, [5]), web-browsing (β ≈ 1.2, [6]) and library loans (β ≈ 1, [7]). In some other cases a monotonous relation has been reported between the user activity and the interevent time exponent, for example in short-message communication (β∈(1.5,2.1), [8]) or in watching movies (β∈(1.5,2.7), [9]). In a recent paper there can be found a distribution of exponents of various channels of communications [10]. These observations make it necessary to find a model in which the interevent time exponent is tunable.

Similar bursty behaviour has been observed in other natural phenomena, for example in neuron-firing sequences [11] or in the timings of earthquakes [12]. The interevent time distribution does not give us any information about the dependence among the consecutive actions. Correlation between events is usually characterized by the autocorrelation function of the timings of detected events. Bursty behaviour is often accompanied by power-law decaying autocorrelation function [13], which is usually thought to indicate long-range dependence, see e.g. [14]. However, time series with independent power-law distributed interevent times show long-range correlations [15, 16].

This paper is organized as follows. We start with introducing a task-queuing model, which has an advantage compared to the Barabási-model [17], namely that the observable is not the waiting time of an action (from adding the task to the list until executing it) but the interevent time between similar activities. We determine the asymptotic decay of the interevent time distribution in a simple limit of the model. We give a simple proof of the scaling law between the exponents of the interevent time distribution and the autocorrelation function based on Tauberian theorems. Finally, we demonstrate that the scaling law can be violated if the interevent times are long-range dependent.

We assume that people have a task list of size N and they choose their forthcoming activity from this list. The list is ordered in the sense that the probability w_i of choosing the ith activity from the list is decreasing as a function of position i. The task at the chosen position jumps to the front (triggering the corresponding activity) and pushes the tasks that preceded it to the right (figure 1). This mechanism is responsible for producing the bursty behaviour because once a person starts to do an activity, that is going to have high priority for a while.

Zoom In Zoom Out Reset image size

The model is capable of covering many types of activities but now we only concentrate on one type. The tasks of this type are marked with A in figure 1 (e.g. watching movies). For the sake of simplicity we assume that the list contains only one single item of type A, all the others are activities of type B. It turns out that this restriction is important for the interevent time distribution but irrelevant for the autocorrelation function.

We show in our paper that this model produces power-law decaying interevent time distribution for a wide variety of the w_i distributions. First we analyse the model with power-law decaying priority distribution w_i, second we analyse the effect of exponentially decaying priority distribution and finally we discuss the stretched exponential case. If the list is finite, the power-law regime is followed by an exponential cutoff. The cutoff is the consequence of reaching the end of the list, from which a geometrically distributed waiting time follows: P_wt(t) ∼ (1 − w_N)^t. A marginal result of section 4 shows that the expectation value of the interevent time is equal to the length of the list independently on the priority distribution (as long as ergodicity is maintained). For the sake of simplicity we determine the exponent of power-law decaying region in the case of an infinitely long list where the exponential cutoff does not appear.

The interevent time is equal to the recurrence time of the first element of the list. Let q(n,t) (n ⩾ 1) denote the probability of finding the observed element at position n after t timesteps without any recurrences up to time t. We emphasize here that q(n,t) is not a conditional probability and that $1-\sum _{n=1}^{\infty } q(n,t)$ gives the survival function of the interevent time distribution. The initial condition is set as q(n,1) = (1 − w₁)δ_n,2, that is, at the first step the observed element moves from the front of the list to the second position with probability (1 − w₁) (otherwise a recurrence would occur). The restriction not to recur is important because this makes large jumps to the front of the list forbidden for the observed element. Time evolution is given by

$$\begin{equation} q(1,t) \equiv 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} q(n,t+1)=(1-Q_{n-1})q(n,t)+Q_{n-1}q(n-1,t) , \end{equation} \tag{ 2 }$$

where $Q_{n}=1-\sum _{k=1}^{n}w_k$ is the survival function of the priority distribution. Equation (2) expresses that the observed element can be found at position n either if it was there in the previous time step and we choose a position smaller than n or it was at position n − 1 and we choose an activity at position k > n. Our aim is to determine the asymptotic behaviour of $\sum _{n=1}^{\infty } q(n,t)$. The main trick we use in analysing the recursion above is to consider the n = const levels first instead of the t = const levels which intuitively one would do. At every step the time coordinate gets increased while the position of the observed element might remain unchanged or get increased by one as well (figure 2). The path of the element in the (n,t) plane can be divided into sections that start with a step on the bias and are followed by some steps upwards (which can be zero as well). These sections can be characterized by their height-difference, which can take values from $1,2,\ldots$. These height differences are independent and (optimistic) geometrically distributed with parameter depending on the position. Let ξ_k be independent geometrically distributed random variables with parameter Q_k, i.e. $\mathbb {P}(\xi _i=\tau )=(1-Q_{i})^{\tau -1}Q_{i}$. With n fixed q(n,t) is the probability that we find the element at the nth position at time t (without any recurrences). This corresponds to paths with n − 1 steps to the right (started from position 1 at time 0) and q(n,t) is the distribution mass function of the sum of the random heights:

$$\begin{equation} q(n,t)=\mathbb{P}\left(\sum_{i=1}^{n-1}\xi_i = t\right). \end{equation} \tag{ 3 }$$

We could go so far without specifying the priority distribution, now we have to specify Q_n to continue.

Zoom In Zoom Out Reset image size

3.1. Power-law decaying priorities

Here we analyse the model in which the survival function of the priority distribution is power-law decaying, i.e. $Q_{n}= c n^{-\sigma +1} + \mathcal {O}(n^{-\sigma })$. $\mathcal {O}(n^{\omega })$ means that term is asymptotically at most of the order of n^ω. In this case w_i is also power-law decaying with exponent σ. Though the random variables ξ_k are not identically distributed—by checking the Lyapunov condition—one can show that the central limit theorem holds for this situation (see section A.2 of the appendix). In this approximation q(n,t) is Gaussian in variable t with mean $\sum _{i=1}^{n-1}\mathbb {E}\left (\xi _i\right )$ and variance $\sum _{i=1}^{n-1}\mathbb {D}^2\left (\xi _i\right )$, where $\mathbb {E}\left (\xi _i\right )=Q_i^{-1}$ and $\mathbb {D}^2\left (\xi _i\right )=Q_i^{-2}(1-Q_i)$. Using integral approximation to evaluate the sums and keeping only the highest order terms in n yields

$$\begin{equation} q(n,t) \stackrel{\mathrm{CLT}}{\approx} \frac{c\sqrt{{2\sigma-1}}}{\sqrt{2 \pi} n^{\sigma-1/2}} \exp \left\{ -\frac{c^2(2\sigma-1)}{2} \frac{\left(t-\frac{n^{\sigma}}{c\sigma}\right)^2}{n^{2\sigma-1}} \right\} . \end{equation} \tag{ 4 }$$

This formula shows that the probability of finding an element at position n at time t is centred on the $t=\frac {n^{\sigma }}{c\sigma }$ curve which is in agreement with the numerical results (figure 3).

Zoom In Zoom Out Reset image size

The sum of q(n,t) in variable n is the probability of the observed element has not recurred up to time t. We approximate this sum by integral and we apply the following substitution (with new variable r):

$$\begin{equation} n=(c\sigma t)^{1/\sigma} - r (c\sigma t) ^{1/(2\sigma)} . \end{equation} \tag{ 5 }$$

For γ > 0 (here γ may refer to σ − 1/2,σ,2σ − 1):

$$\begin{equation} n^\gamma = (c\sigma t)^{\frac{\gamma}{\sigma}} - \gamma (c\sigma t)^{\frac{2\gamma-1}{2\sigma}} r + o(t^{\frac{2\gamma-1}{2\sigma}})\, , \end{equation} \tag{ 6 }$$

where o(n^ω) means that term is asymptotically of strictly smaller order than n^ω. From equations (4)–(6) it follows that

$$\begin{equation} \int_0^\infty q(n,t)\, \mathrm{d}n = c^{\frac{1}{\sigma}}(\sigma t)^{\frac{1}{\sigma}-1} + o (t^{\frac{1}{\sigma}-1}). \end{equation} \tag{ 7 }$$

Differentiating this equation gives the first main result of our paper, P_ie(t) ∼ t^−β with $\beta =2-\frac {1}{\sigma }$.

3.2. Exponentially decaying priorities

Now we study the model with Q_n = ce^−λn priority survival function. In contrast to the previous case the central limit theorem cannot be used, because the last term in $\sum _{i=1}^{n-1}\xi _i$ is comparable with the complete sum.

However, we can construct a limit theorem even for this situation. First of all, we approximate the geometrically distributed random variables ξ_k ∼ Geom(1 − Q_k) by exponential ones: $\tilde {\xi }_k \sim {\mathrm { Exp}}(Q_k)$. We apply the scaling property of the exponential distribution to express $\tilde {\xi }_k$ by i.i.d. exponential random variables: $\tilde {\xi }_k \sim \mathrm {e}^{\lambda k} \eta _k$, where η_k ∼ Exp(c). In equation (3) we need the distribution of sums of the first n − 1 ξ:

$$\begin{equation} \sum_{i=1}^{n-1}\tilde{\xi}_i \sim \mathrm{e}^{n \lambda} X_{n-1}, \end{equation} \tag{ 8 }$$

$$\begin{equation} X_{n-1}=\sum_{k=1}^{n} \mathrm{e}^{-k\lambda}\eta_k . \end{equation} \tag{ 9 }$$

X_n has a well defined limit as n → ∞, which we denote by X. The probability density function of this (limit) random variable is denoted by ψ(x) (and is given explicitly in section A.1 of the appendix). The finite sum X_n in equation (8) can be approximated by X, yielding $\mathbb {P}(\sum _{i=1}^{n-1}\tilde {\xi }_i < t)\approx \mathbb {P}(X < \mathrm {e}^{-\lambda n}(t))$ and

$$\begin{equation} \fl \int_{0}^{\infty} q(n,t)\, \mathrm{d}n \approx \int_{0}^{\infty} \mathrm{e}^{-\lambda n} \psi(t \mathrm{e}^{-\lambda n})\,\mathrm{d}n = \frac{1}{\lambda t}\int_0^t\psi(y) \,\mathrm{d}y =\frac{1}{\lambda t} + o(t^{-1}) , \end{equation} \tag{ 10 }$$

since $\int _0^t\psi (y)\, \mathrm {d}y$ tends to 1 as t → ∞. By differentiating we get β = 2 independently on parameter λ.

3.3. Stretched exponentially decaying priorities

The last investigated priority distribution is the stretched exponential, i.e. Q_n = c e^−λnν .

Faster than exponential decay (ν > 1).

In this case $\sum _{i=1}^{n-1}{\xi _i}$ is totally dominated by the last term in the sum, hence we approximate the sum with this term

$$\begin{equation} q(n,t) \approx c\, \mathrm{e}^{-\lambda n^{\nu}} \mathrm{e}^{-t c \,\mathrm{e}^{-\lambda n^{\nu}}}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \int_0^{\infty} q(n,t)\, \mathrm{d}n = \frac{1}{t\lambda^{1/\nu} \nu} \int_{0}^{t c} \mathrm{e}^{-y} (\log(tc)-\log(y))^{-(1-1/\nu)} \,\mathrm{d}y \end{equation} \tag{ 12 }$$

$$\begin{equation}\hspace{6pc} =\frac{1}{\lambda^{1/\nu} \nu} \frac{1}{t (\log(tc))^{1-1/\nu}}+ o\left(\frac{1}{t (\log(tc))^{1-1/\nu}}\right). \end{equation} \tag{ 13 }$$

This means that the β = 2 exponent holds for this case but a logarithmic correction appears.

Slower than exponential decay (ν < 1).

In this case the increment of Q_n is just small enough for the central limit theorem to be still applicable (see the appendix). We approximate the mean and variance of $\sum _{i=1}^{n-1}\xi _i$ by integral and we use the property $\int _{0}^{a} \mathrm {e}^{y} y^{1/\nu -1} \stackrel {a\rightarrow \infty }{\approx } \mathrm {e}^{a} a^{1/\nu -1}$. With these we get

$$\begin{equation} q(n,t) \approx \sqrt{\frac{c^2 \lambda \nu}{\pi}}\, \mathrm{e}^{-\lambda n^{\nu}}n^{-\frac{1-\nu}{2}} \exp\left\{-\frac{\left(t-\mathrm{e}^{\lambda n^{\nu}}n^{1-\nu}/(c \lambda \nu)\right)^2}{2 \frac{1}{2 c^2\lambda \nu} \mathrm{e}^{2\lambda n^{\nu}}n^{1-\nu}}\right\} . \end{equation} \tag{ 14 }$$

The t → ∞ asymptotic behaviour of $\int\! q(n,t)\, \mathrm {d}n$ can be determined by introducing a proper new variable

$$\begin{equation} y= n^{\frac{1-\nu}{2}}-c \lambda \nu t \,\mathrm{e}^{-\lambda n^{\nu}}n^{\frac{\nu-1}{2}} . \end{equation} \tag{ 15 }$$

Keeping the leading order term only, one finds that the survival function of the interevent time distribution is

$$\begin{equation} \int_{0}^{\infty} q(n,t)\, \mathrm{d}n =\frac{1}{\lambda^{1/\nu}\nu} \frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}+o\left(\frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}\right) . \end{equation} \tag{ 16 }$$

For more details of the calculations see section A.3 of the appendix. We have got β = 2 with logarithmic correction again. While for ν > 1 the correction makes the decay of the distribution faster than the pure power law, for ν < 1 it is slightly slowed down.

Another characteristic property of the time series is autocorrelation function. Let X(t) denote the indicator variable of the observed activity: X(t) = 1 if the observed activity is at the first position of the list at time t (i.e. an event happened) and X(t) = 0 otherwise. We define the autocorrelation function as usual,

$$\begin{equation} \mathcal{A}(t)=\frac{\left\langle \mathbb{E}\left[X(s)X(t+s)\right]\right\rangle_s-\left\langle \mathbb{E}\left[X(s)\right] \right\rangle_s^2}{\left\langle \mathbb{E}\left[X(s)\right] \right\rangle_s-\left\langle \mathbb{E}\left[X(s)\right] \right\rangle_s^2}. \end{equation} \tag{ 17 }$$

The model defines a Markov-chain and the state space is the space of permutations of the list. The transition matrix for a list of length N is given by

$$\begin{equation} M=\left(\begin{array}{@{}ccccc@{}} 1-Q_1 & Q_1 & 0 & \cdots & 0\\ Q_1-Q_2 & 1-Q_1 & Q_2 & 0 & \vdots\\ Q_2-Q_3 & 0 & 1-Q_2 & Q_3 & 0\\ \vdots & \vdots & 0 & \ddots & Q_{N-1}\\ Q_{N-1}-Q_{N} & 0 & \cdots & 0 & 1-Q_{N-1} \end{array}\right), \end{equation} \tag{ 18 }$$

where Q_N = 0. This matrix is doubly stochastic, hence $\left \langle \mathbb {E}\left [X(t)\right ]\right \rangle _t=\frac {1}{N}$. According to the recurrence formula of Marc Kac [18] we conclude that the expectation value of the interevent time is equal to the size of the list independently on the priority distribution (as long as ergodicity is maintained). The autocorrelation function can be written in simpler form making use of the knowledge of $\left \langle \mathbb {E}\left [X(t)\right ]\right \rangle _t$,

$$\begin{equation} \mathcal{A}(t)=\frac{\mathbb{P}\left(X(t)=1 | X(0)=1\right)-\frac{1}{N}}{1-\frac{1}{N}} . \end{equation} \tag{ 19 }$$

The probability of finding the activity at the first position can be calculated numerically by successive application of the Markov-chain transition matrix (18). For power-law decaying priority distributions numerical computations show that the autocorrelation function is power-law decaying with an exponential cutoff (figure 4). Given σ, the autocorrelation functions for various list sizes can be rescaled to collapse into a single curve (figure 4, inset). This property can be written in a mathematical way:

$$\begin{equation} \mathcal{A}(t,N,\sigma)=N^{-\delta} f_{\!\sigma} \left( \frac{t}{N^{\gamma}} \right) =t^{-\frac{\delta}{\gamma}} \skew3\tilde{f}_{\sigma} \left( \frac{t}{N^{\gamma}} \right). \end{equation} \tag{ 20 }$$

The exponents used for rescaling the autocorrelation functions are listed in table 1. Indicating that in the N → ∞ limit the autocorrelation function is power-law decaying with exponent $\alpha =\frac {1}{\sigma }$ for σ ⩾ 2. With the exact $\beta =2-\frac {1}{\sigma }$ result this can be combined into a scaling law: α + β = 2.

Zoom In Zoom Out Reset image size

If the priority distribution decays faster than power-law, then β = 2 and α = 0. This means that the autocorrelation function decays slower than any power. This is the case when the priority distribution decays like a (stretched) exponential function.

The essential properties of the model for the scaling relation are that the interevent times are independent and power-law decaying. Let T denote the set of recurrence times and let τ be an interevent time. The autocorrelation function can be written in the following form:

$$\begin{equation} \mathcal{A}(t)=\frac{\mathbb{P} \left( t \in T \right)-\frac{1}{\mathbb{E}\left[\tau\right]}}{1-\frac{1}{\mathbb{E}\left[\tau\right]}} , \end{equation} \tag{ 21 }$$

which simplifies to $\mathcal {A}(t)=\mathbb {P} \left ( t \in T \right )$ if β ⩽ 2. The Laplace transform of the autocorrelation function can be expressed by the Laplace transform of the interevent time distribution

$$\begin{equation} g(\lambda) \equiv \sum_{t=0}^{\infty}{\mathcal{A}(t)\, \mathrm{e}^{-\lambda t}}=\left( 1-\mathbb{E} \left[ \mathrm{e}^{- \lambda \tau} \right] \right)^{-1}, \end{equation} \tag{ 22 }$$

where we used the property that the interevent times are independent and identically distributed. Tauberian and Abelian theorems connect the asymptotics of a function with the asymptotics of its Laplace transform [20]. Applying Abelian theorem to the right side of the last equation results in g(λ) ∼ λ^1−β. Then applying the Tauberian theorem yields $\mathcal {A}(t) \sim t^{\beta -2}$ or

$$\begin{equation} \alpha+\beta=2. \end{equation} \tag{ 23 }$$

To be precise we had to use an extended version of the Abelian theorem, which can be derived from the original theorem using integration by parts. With similar train of thought the scaling law can be extended to the 2 < β < 3 region where β − α = 2 holds. These results are in agreement with [16, 19].

Independence of the interevent times was important in the proof of the scaling law. The violation of the scaling relation indicates dependence among the interevent times but we emphasize that the opposite direction does not necessarily hold. In this section we investigate two models that are characterized by dependent interevent times. One of them satisfies the scaling relation (23) and the other does not.

6.1. A model with weak dependence

Our first example is the two-state model of [13], because that model was introduced to capture deep correlations, nevertheless, it obeys the scaling law. We start with a brief introduction of the model. The system alternates between a normal state A, in which the events are separated by longer time intervals, and an excited state B, which performs events with higher frequency. After each event executed in state A the system may switch to state B with probability π or remain in state A with probability 1 − π. In contrast, the excited state has a memory. That is, if the system had executed more events in state B since the last event in the normal state, the transition to state A would be less probable. The probability of performing one more event in state B is given by $p(n)=\left (n/(n+1)\right )^\nu$, where n is the number of events in the current excited state. This memory function was introduced to model the empirically observed power-law decay in the bursty event number distribution (for the definition see the reference). The dynamics is illustrated in figure 5(a) for a better understanding. In both states the interevent times are generated from a reinforcement process resulting independent random interevent times τ_A and τ_B with power-law decaying distributions. We denote the corresponding exponents by β_A and β_B. The parameters of the model were the following: π = 0.1, ν = 2.0, β_A = 1.3, β_B = 5.0, and it was found that the model satisfies the scaling law with exponents β = 1.3 and α = 0.7.

Zoom In Zoom Out Reset image size

The number of events performed in a single period of state A or B are i.i.d. random variables N_A or N_B. N_A is geometrically distributed with parameter π. The distribution of N_B can be easily obtained from p(n), and it is power-law decaying with exponent 1 + ν. Making use of these new random variables the dynamics simplifies: the system executes N_A events with interevent times distributed as τ_A, then switches to the excited state and executes N_B events with interevent times τ_B, then switches back to the normal state, and so on.

It is clear that the decay of the interevent time distribution of the whole process is determined by the smallest of the exponents β_A and β_B. Since both of $\mathbb {E}\left [\tau _B\right ]$ and $\mathbb {E} [N_B]$ are finite, the periods of state B have a characteristic length $\xi _t=\mathbb {E} [\tau _B]\mathbb {E}[N_B]$. Within the periods of state A there is no characteristic time-scale, hence by looking on the time series at a larger scale the effective length of the B periods will be shorter (figure 5(b))). Therefore the existence of the B periods is irrelevant for the asymptotic decay of the autocorrelation function giving rise to the scaling law.

6.2. A model with long-range dependence

Another simple way to introduce dependence between the interevent times is to construct a Markov-chain, i.e. the next interevent time depends only on the actual interevent time. As a counterexample to the scaling law we constructed a long-range dependent set of interevent times by Metropolis-algorithm. The base of the algorithm is constructing a Markov chain on the integers that has power-law decaying stationary distribution P(x) ∼ x^−β. The algorithm uses a proposal density (transition rate) Q(x',x_n) which generates a proposal sample from the current value of the interevent time. This sample is accepted for the next value with probability

$$\begin{equation} \alpha(x',x_n)=\min\left\{\!\frac{Q(x_n,x')P(x')}{Q(x',x_n)P(x_n)},1\!\right\} , \end{equation} \tag{ 24 }$$

otherwise the previous value is repeated. To generate a power-law distributed sample with long-range dependence the mixing of the Markov chain should be slow, i.e. the gap in the spectrum of the Markov chain should vanish. In this order we allow only small differences between consecutive interevent times: $Q(x',x_n)=\frac {1}{2D}\mathbb {I}(x_n-D,x_n+D)$ or equivalently x' as a random variable is discrete uniformly distributed: x' ∼ DU[x_n − D,x_n + D]. The simulation results (figure 6) show that the autocorrelation function decays slower than it would be assumed from the scaling law (23). This example shows that the scaling relation (23) can be violated if there is long range dependence among the interevent times.

Zoom In Zoom Out Reset image size

A trivial extension of the model could be putting more items of the observed activity into the list. Then this parameter would tune the frequency of doing the activity. In this case the interevent times become dependent and the simulation results on finite lists show that the interevent time distribution decays faster than power-law but decays slower than an exponential function. It is easy to prove that in spite of this the autocorrelation function (17) is independent of the number of the observed activity and remains power-law decaying. This is another example for models in which the scaling law does not hold.

When the list is finite, we have much more freedom in the choice of the priority distribution. However, the method based on expressing q(n,t) as a probability density function of sum of n − 1 independent geometrically distributed random variables is general, and can offer a good base for further calculations, even for finite lists.

The model may be interesting not only for human dynamics but also for the mathematical theory of card shuffling [21]. We can define the time reversed version of the model in which we choose a position from the same distribution as in the original model and we move the first element of the list to the random position. This model has similar properties to the original model, e.g. this model has the same interevent time distribution and autocorrelation function. If we think of the list as a deck of cards, then the time reversed model is a generalization of the top-in-at-random shuffle method. The latter is a primitive model of shuffling cards: the top card is removed and inserted into the deck at a uniformly distributed random position [21].

In a proper model one should take into consideration the dependence among the interevent times besides the interevent time distribution. In human dynamics the latter is slowly decaying and as a consequence of this the autocorrelation function of the interevent times is not well-defined (i.e. the second moment does not exist). However, the autocorrelation function of the time series exists and it might be a good measure for long range dependence among the interevent times. Time series of messages sent by an individual in an online community are reported to be not more correlated than an independent process with the same interevent time distribution [16]. Similarly, the exponents in the neuron-firing sequences approximately satisfy the scaling law (α = 1.0, β = 1.1 [13]). For these systems the model we studied might be applicable. However, this is not always the case. For example, the autocorrelation function of the e-mail sequence decays slower (α = 0.75 [13]) than it should do estimated from the scaling law. This indicates long range dependence in the time series in addition to the power-law decaying interevent time distribution. In this case a dependent model should be applied, for example a model based on Metropolis algorithm that is similar to the one we studied as a counterexample to the scaling law. Effects of circadian patterns and inhomogeneity in an ensemble of individuals should also be considered [22].

In real networks interactions between individuals have to be taken into account to reproduce some social phenomena, e.g. temporal motifs [23] observed in a mobile call dataset. Interactions could be incorporated in this model by allowing the actual activity of an individual to modify the priority list of some of his/her neighbour. If this effect is rare and can be considered as a perturbation, our results on the dynamics of the list could be a starting point to further calculations covering for example information flow in a network.

The model discussed here was introduced by Chaoming Song and Dashun Wang. We are grateful to them, László Barabási and Márton Karsai for discussions. This project was partially supported by FP7 ICTeCollective Project (no. 238597) and by the Hungarian Scientific Research Fund (OTKA, no. K100473).

A.1. Limit theorem for a geometric series of exponentially distributed random variables

In this section we determine the probability density function of the limit random variable X. In the main text we defined X_n and X as follows:

$$\begin{equation} X_{n} \sim \sum_{k=1}^{n} \mathrm{e}^{-k\lambda}\eta_k, \end{equation} \tag{ A.1 }$$

$$\begin{equation} X \sim \sum_{k=1}^{\infty} \mathrm{e}^{-k\lambda}\eta_k , \end{equation} \tag{ A.2 }$$

where η_k ∼ Exp(c) are i.i.d. exponential random variables. The probability density function of sums of independent exponential random variables with different parameters can be expressed [24]. Substituting the proper parameters into that formula one gets

$$\begin{equation} f_{X_n}(x)=\prod_{j=1}^{n}{c\,\mathrm{e}^{j\lambda}}\sum_{j=1}^{n}\frac{\mathrm{e}^{-c \,\mathrm{e}^{j\lambda}x}}{\prod_{k=1,k\neq j}^{n}(\mathrm{e}^{k\lambda}-\mathrm{e}^{j\lambda})} . \end{equation} \tag{ A.3 }$$

With some trivial algebraic manipulations this can be written in the following form:

$$\begin{equation} f_{X_n}(x)=c \sum_{j=1}^{n}(-1)^{j-1}\frac{\mathrm{e}^{-j(j-3)\frac{\lambda}{2}}}{\prod_{k=1}^{j-1}(1-\mathrm{e}^{-k\lambda})\prod_{k=1}^{n-j}(1-\mathrm{e}^{-k\lambda})}\mathrm{e}^{-c \,\mathrm{e}^{j\lambda}x} . \end{equation} \tag{ A.4 }$$

The probability density function in the n → ∞ limit reads

$$\begin{equation} f_{X}(x)=\frac{c}{\prod_{k=1}^{\infty}(1-\mathrm{e}^{-k\lambda})} \sum_{j=1}^{\infty}(-1)^{j-1}\frac{\mathrm{e}^{-j(j-3)\frac{\lambda}{2}}}{\prod_{k=1}^{j-1}(1-\mathrm{e}^{-k\lambda})}\mathrm{e}^{-c \,\mathrm{e}^{j\lambda}x} . \end{equation} \tag{ A.5 }$$

The products $\prod _{k=1}^{s}(1-\mathrm {e}^{-k\lambda })$ are usually cited as q-Pochhammer symbols (e^−λ,e^−λ)_s and are convergent in the s → ∞ limit. Equation (A.5) shows clearly that the asymptotic behaviour of the limit random variable X is determined by the first term in the sum.

A.2. Central limit theorem for the power-law decaying case

The Lyapunov condition for the central limit theorem reads

$$\begin{equation} \frac{\sum_{i=1}^n\mathbb{E} \left|\xi_i-\mathbb{E}\xi_i\right|^{2+\delta}}{(\sum_{i=1}^n\mathbb{D}^2\xi_i)^{1+\delta/2}} \rightarrow 0 \quad {\rm as} \ n \quad \rightarrow \infty. \end{equation} \tag{ A.6 }$$

If this condition is fulfilled for any δ, then the central limit theorem can be used. We test this condition at δ = 1 with the following moments:

$$\begin{equation} \mathbb{E}(\xi_k) = \frac{1}{Q_k}=\frac{1}{c} k^{\sigma-1}, \end{equation} \tag{ A.7 }$$

$$\begin{equation} \mathbb{D}^2(\xi_k) = \frac{1}{Q_k^2}=\frac{1}{c^2} k^{2\sigma-2}, \end{equation} \tag{ A.8 }$$

$$\begin{equation} \mathbb{E}(|\xi_k-\mathbb{E}\xi_k|^3) = \frac{12 - 2 e}{e Q_k^3}=\frac{12 - 2 e}{e c^3} k^{3\sigma-3} . \end{equation} \tag{ A.9 }$$

Substituting these and approximating the sums by integrals yield

$$\begin{equation} \frac{\sum\mathbb{E} |\xi_i-\mathbb{E}\xi_i|^3}{(\sum\mathbb{D}^2\xi_i)^{3/2}} \sim \frac{\int_0^{n}k^{3\sigma-3}\,\mathrm{d}k}{\left(\int_0^{n}k^{2\sigma-2}\,\mathrm{d}k\right)^{3/2}} \sim n^{-1/2}\rightarrow 0 \quad {\rm as}\quad n \rightarrow \infty. \end{equation} \tag{ A.10 }$$

A.3. Central limit theorem for the stretched exponential case (ν < 1)

We test the Lyapunov condition (A.6) at δ = 1 similarly to the previous case. The moments are the following:

$$\begin{equation} \mathbb{E}(\xi_k) = \frac{1}{Q_k}=\frac{1}{c} \mathrm{e}^{\lambda k^{\nu}}, \end{equation} \tag{ A.11 }$$

$$\begin{equation} \mathbb{D}^2(\xi_k) = \frac{1}{Q_k^2}=\frac{1}{c^2} \mathrm{e}^{2\lambda k^{\nu}}, \end{equation} \tag{ A.12 }$$

$$\begin{equation} \mathbb{E}(|\xi_k-\mathbb{E}\xi_k|^3) = \frac{12 - 2 e}{e Q_k^3}=\frac{12 - 2 e}{e c^3} \mathrm{e}^{3\lambda k^{\nu}} . \end{equation} \tag{ A.13 }$$

We make integral approximation of the sums of the moments above:

$$\begin{equation}\fl \sum_{k=1}^{n}{\mathrm{e}^{a \lambda k^{\nu}}}\approx \int_{0}^{n}{\mathrm{e}^{a \lambda k^{\nu}}} \mathrm{d}k = \frac{1}{\nu \lambda^{1/\nu}}\int_{0}^{a \lambda n^{\nu}}\mathrm{e}^{y} y^{1/\nu-1} \approx \frac{1}{a\lambda \nu}\mathrm{e}^{a \lambda n^{\nu}} n^{1-\nu} . \end{equation} \tag{ A.14 }$$

The last approximation comes from the property that $\int _{0}^{a} \mathrm {e}^{y} y^{1/\nu -1} \stackrel {a\rightarrow \infty }{\approx } \mathrm {e}^{a} a^{1/\nu -1}$. With these we have

$$\begin{equation} \frac{\sum\mathbb{E} |\xi_i-\mathbb{E}\xi_i|^3}{(\sum\mathbb{D}^2\xi_i)^{3/2}}\sim \frac{\mathrm{e}^{3\lambda n^{\nu}}n^{1-\nu}}{(\mathrm{e}^{2\lambda n^{\nu}}n^{1-\nu})^{3/2}}=n^{-\frac{1-\nu}{2}} \rightarrow 0 \quad {\rm as}\; n\rightarrow \infty \end{equation} \tag{ A.15 }$$

if ν < 1.

Now we give some more details in calculating the β exponent. The mean and variance appearing in the central limit theorem are

$$\begin{equation} \sum_{i=1}^{n-1} \mathbb{E} (\xi_i) \approx \frac{1}{c \lambda \nu}\mathrm{e}^{\lambda n^{\nu}}n^{1-\nu}, \end{equation} \tag{ A.16 }$$

$$\begin{equation} \sum_{i=1}^{n-1} \mathbb{D}^2 (\xi_i) \approx \frac{1}{2 c^2 \lambda \nu} {{\rm e}^{2\lambda n^{\nu}}}n^{1-\nu} \end{equation} \tag{ A.17 }$$

and the approximating probability density function is

$$\begin{equation} q(n,t) \approx \sqrt{\frac{c^2 \lambda \nu}{\pi}} \mathrm{e}^{-\lambda n^{\nu}}n^{-\frac{1-\nu}{2}} \exp\left\{-\frac{\left(t-\mathrm{e}^{\lambda n^{\nu}}n^{1-\nu}/(c \lambda \nu)\right)^2}{2 \frac{1}{2 c^2\lambda \nu} \mathrm{e}^{2\lambda n^{\nu}}n^{1-\nu}}\right\} . \end{equation} \tag{ A.18 }$$

We introduce a new variable

$$\begin{equation} y= n^{\frac{1-\nu}{2}}-c \lambda \nu t\, \mathrm{e}^{-\lambda n^{\nu}}n^{\frac{\nu-1}{2}}, \end{equation} \tag{ A.19 }$$

$$\begin{equation} \mathrm{d}y=\left[\frac{1-\nu}{2}n^{\frac{-\nu+1}{2}}+(n^{\frac{1-\nu}{2}}-y)\left(\lambda \nu n^{\nu-1}+\frac{1-\nu}{2}n^{-1}\right)\right]\mathrm{d}n , \end{equation} \tag{ A.20 }$$

in the second equation $n\equiv n(y)= (\frac {1}{\lambda }\log (c\lambda \nu t))^{1/\nu }+o(\log (c\lambda \nu t)^{1/\nu })+y \, o(\log (c\lambda \nu t)^{1/\nu })$. The dominant term in the Jacobian when t → ∞ is $\lambda \nu n^{\frac {\nu -1}{2}}$. With these the integral $\int _{0}^{\infty }q(n,t)\,\mathrm { d}n$ reads

$$\begin{eqnarray} &&\fl \sqrt{\frac{c^2 \lambda \nu}{\pi}} \int_{-\infty}^{\infty} \mathrm{e}^{-\frac{1}{\lambda \nu} y^2} \left(\underline{\frac{n^{\frac{1-\nu}{2}}}{c \lambda \nu t}}-\frac{y}{c \lambda \nu t}\right) \left[\frac{1-\nu}{2}n^{\frac{-\nu+1}{2}}+(\underline{n^{\frac{1-\nu}{2}}}-y)\left(\underline{\lambda \nu n^{\nu-1}}+\frac{1-\nu}{2}n^{-1}\right)\right]^{-1} \, \mathrm{d}y\nonumber\\ \end{eqnarray} \tag{ A.21 }$$

$$\begin{equation} \fl = \sqrt{\frac{\lambda \nu}{\pi}}\frac{1}{\lambda^{1+1/\nu} \nu^2} \int_{-\infty}^{\infty} \mathrm{e}^{-\frac{1}{\lambda \nu} y^2}\left(\frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}+o\left(\frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}\right)\right) \, \mathrm{d}y \end{equation} \tag{ A.22 }$$

$$\begin{equation} \fl =\frac{1}{\lambda^{1/\nu}\nu} \frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}+o\left(\frac{1}{t \log(c\lambda\nu t)^{1-1/\nu}}\right). \end{equation} \tag{ A.23 }$$

We underlined the terms that determine the decay of the integral for large t.