Long-term temporal correlations in time series in a form of an event sequence have been characterized using an autocorrelation function that often shows a power-law decaying behavior. Such scaling behavior has been mainly accounted for by the heavy-tailed distribution of interevent times, i.e., the time interval between two consecutive events. Yet, little is known about how correlations between consecutive interevent times systematically affect the decaying behavior of the autocorrelation function. Empirical distributions of the burst size, which is the number of events in a cluster of events occurring in a short time window, often show heavy tails, implying that arbitrarily many consecutive interevent times may be correlated with each other. In the present study, we propose a model for generating a time series with arbitrary functional forms of interevent time and burst size distributions. Then, we analytically derive the autocorrelation function for the model time series. In particular, by assuming that the interevent time and burst size are power-law distributed, we derive scaling relations between power-law exponents of the autocorrelation function decay, interevent time distribution, and burst size distribution. These analytical results are confirmed by numerical simulations. Our approach helps to rigorously and analytically understand the effects of correlations between arbitrarily many consecutive interevent times on the decaying behavior of the autocorrelation function.

In general, events in complex systems are correlated with each other. The autocorrelation function for the event series measures the overall diminishing impact of one event at a time on another event in the future (or past). The autocorrelation function in complex systems typically shows a long tail, meaning a long-lasting memory. We mathematically solve the autocorrelation function to understand the origin of such long tails in terms of not only the heavy-tailed distribution of the time interval between two successive events, namely, interevent time, but also correlations between arbitrarily many successive interevent times.

## I. INTRODUCTION

^{1–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–9 and references therein. One of the most commonly used measurements is an autocorrelation function.

^{10,11}Precisely, the autocorrelation function for a time series $x(t)$ is defined with a time lag $ t d$ as

^{8,12–14}For the time series with long-term correlations, the autocorrelation function typically decays in a power-law form with a decay exponent $\gamma $ such that

^{15}as well as to the Hurst exponent $H$

^{16}and its generalizations.

^{17–21}

^{22}Temporal correlations in such event sequences have been characterized by the time interval between two consecutive events, namely, an interevent time. In many empirical data sets, interevent time distributions, denoted by $P(\tau )$, have heavy tails;

^{23–36}in particular, $P(\tau )$ often shows a power-law tail as

^{8}Lowen and Teich derived the analytical solution of the power spectral density for a renewal process governed by a power-law interevent time distribution,

^{37,38}concluding that the decay exponent $\gamma $ in Eq. (2) is solely determined by the interevent time exponent $\alpha $ in Eq. (3) as

^{39–41}

^{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 interevent time 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 interevent time 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

^{42}Similar observations have been made in other data.

^{43–48}These findings immediately raise an important question: how does the autocorrelation function decay power-law exponent $\gamma $ depend on the interevent time 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 interevent times may be correlated with each other. It is worth noting that correlations between two consecutive interevent times 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 interevent times have been modeled by means of, e.g., the two-state Markov chain,^{42} self-exciting point processes,^{56} the interevent time permutation method,^{57} and a model inspired by the burst-tree decomposition method.^{48} Although scaling behaviors of the autocorrelation function were studied in some of mentioned works,^{42,53,56,57} the scaling relation $\gamma (\alpha ,\beta )$ has not been clearly understood due to the lack of analytical solutions of the autocorrelation function.

In this work, we devise a model for generating a time series with arbitrary functional forms of the interevent time distribution and burst size distribution. Assuming power-law tails for interevent time and burst size distributions, our model generates correlations between an arbitrary number of consecutive interevent times. By theoretically analyzing the model, we derive asymptotically exact solutions of the autocorrelation function from the model time series, thereby enabling us to find the scaling relation $\gamma (\alpha ,\beta )$. These analytical results will help us to better understand temporal scaling behaviors in empirical time series.

The paper is organized as follows. In Sec. II, we introduce the model with arbitrary functional forms of interevent time and burst size distributions. In Sec. III, we provide an analytical framework for the derivation of the autocorrelation function for the model time series. In Sec. IV, by assuming power-law distributions of interevent times and burst sizes, we derive analytical solutions of the autocorrelation function, 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.

## II. MODEL

Let us introduce the following model for generating event sequences ${x(1),\u2026,x(T)}$ in discrete time, where $T$ is the number of discrete times, for given interevent time distribution and burst size distribution. By definition, $x(t)=1$ if an event occurs at time $t\u2208{1,\u2026,T}$, and $x(t)=0$ otherwise. By construction, the minimum interevent time is $ \tau 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 interevent time distribution for $\tau \u22652$, 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,\u2026, b 1$. Then, we randomly draw an interevent time $ \tau 1$ from $\psi (\tau )$ to set $x(t)=0$ for $t= b 1+1,\u2026, b 1+ \tau 1\u22121$. Note that $ \tau 1\u22652$. We draw another burst size $ b 2$ from $Q(b)$ and another interevent time $ \tau 2$ from $\psi (\tau )$, respectively, to set $x(t)=1$ for $t= b 1+ \tau 1,\u2026, b 1+ \tau 1+ b 2\u22121$ and $x(t)=0$ for $t= b 1+ \tau 1+ b 2,\u2026, b 1+ \tau 1+ b 2+ \tau 2\u22121$. We repeat this procedure until $t$ reaches $T$. See Fig. 1(a) for an example.

^{8}enabling us to treat $\alpha $ and $\beta $ as independent parameters to some extent. Later, we will assume that the interevent time distribution has a power-law tail for $\tau \u22652$, while the burst size distribution is a pure power-law.

## III. ANALYTICAL FRAMEWORK

Let us consider the case in which $x(t)x(t+ t d)$ is non-zero, i.e., $x(t)=1$ and $x(t+ t d)=1$. As depicted in Fig. 1(a), the time series $x(t)$ in a period of $[t,t+ t d]$ is typically composed of several alternating bursts and interevent times larger than one. Here, the consecutive interevent times of $\tau =1$ forms a burst and the sum of such interevent times of $\tau =1$ is equal to the burst size minus one. Therefore, the time lag $ t d$ is written as a sum of burst sizes (each minus one) of bursts appeared in the period of $[t,t+ t d]$ and interevent times 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 d$.

^{8}and the probability distribution of $r$ is given by

## IV. POWER-LAW CASE

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

### A. Case with *α* = *β* = 2

*α*=

*β*= 2

^{60}Using Eq. (15) one obtains

### B. Case with *α* ≠ 2, *β* = 2

*α*≠ 2,

*β*= 2

### C. Case with **1 < ***α*, *β* < 2

*α*,

*β*< 2

### D. Case with *α*, *β* > 2

*α*,

*β*> 2

### E. Case with **1 < ***α* < 2, *β* > 2

*α*< 2,

*β*> 2

### F. Numerical simulation

For the simulations, we use $T=5\xd7 10 7$ and $ \tau \xaf c= b c= 10 6$ to generate $100$ different event sequences $x(t)$ for each combination of $\alpha $ and $\beta $. Then, their autocorrelation functions are calculated using Eq. (89). As shown in Fig. 3, simulation results in terms of $ A ~ ( t d)$ for several combinations of parameters of $\alpha $ and $\beta $ are in good agreement with corresponding analytical solutions of $A( t d)$ in most cases. In some cases, we observe systematic deviations of analytical solutions from the simulation results, which may be due to ignorance of higher-order terms when deriving analytical solutions and/or finite-size effects of simulations.

## V. CONCLUSION

To study the combined effects of the interevent time distribution $P(\tau )$ and the burst size distribution $Q(b)$ on the autocorrelation function $A( t d)$, we have devised a model for generating time series using $P(\tau )$ and $Q(b)$ as inputs. Our model is simple but takes correlations between an arbitrary number of consecutive interevent times into account in terms of bursty trains.^{42} We are primarily interested in temporal scaling behaviors observed in $A( t d)\u221d t d \u2212 \gamma $ when $P(\tau )\u221d \tau \u2212 \alpha $ (except at $\tau =1$) and $Q(b)\u221d b \u2212 \beta $ are assumed. We have derived the analytical solutions of $A( t d)$ for arbitrary values of interevent time power-law exponent $\alpha >1$ and burst-size power-law exponent $\beta >1$ to obtain the autocorrelation function decay power-law exponent $\gamma $ as a function of $\alpha $ and $\beta $ [Eq. (86); see also Fig. 2].

We remark that our model has assumed that interevent times with $\tau \u22652$ and burst sizes in the time series are independent of each other. However, there are observations indicating the presence of correlations between consecutive burst sizes and even higher-order temporal structure.^{48} Thus, it would be interesting to see whether such higher-order structure affects the decaying behavior of the autocorrelation function.

So far, we have focused on the analysis of the single time series observed for a single phenomenon or for a system as a whole. However, there are other complex systems in which each element of the system has its own bursty activity pattern or a pair of elements have their own bursty interaction pattern, such as calling patterns between mobile phone users.^{61,62} In recent years, such systems have been studied in the framework of temporal networks,^{9,63,64} where interaction events are heterogeneously distributed among elements as well as over the time axis. Modeling temporal networks is important to understand collective dynamics, such as spreading or diffusion,^{65,66} taking place in those networks. Some recent efforts for modeling temporal networks are mostly concerned with heavy-tailed interevent time distributions for each element or a pair of elements.^{67,68} Our model can be extended to generate more realistic temporal networks in which activity patterns of elements or interaction patterns between elements are characterized by bursty time series with higher-order temporal structure beyond the interevent time distribution.

## ACKNOWLEDGMENTS

H.-H.J. and T.B. acknowledge financial support by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2022R1A2C1007358). N.M. acknowledges financial support by the Japan Science and Technology Agency (JST) Moonshot R&D (under Grant No. JPMJMS2021), the National Science Foundation (under Grant Nos. 2052720 and 2204936), and JSPS KAKENHI (under Grant Nos. JP 21H04595 and 23H03414).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Hang-Hyun Jo:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Tibebe Birhanu:** Formal analysis (supporting); Validation (equal); Writing – review & editing (equal). **Naoki Masuda:** Formal analysis (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX A: DERIVATION OF THE NORMALIZATION CONSTANT IN EQ. (16)

### APPENDIX B: ANALYSIS FOR THE CASE WITH EXPONENTIAL DISTRIBUTIONS OF REDUCED INTEREVENT TIMES AND BURST SIZES

### APPENDIX C: DERIVATION OF THE LAPLACE TRANSFORMS OF $\psi ( \tau \xaf )$ WITH $\alpha =2$ AND $Q(b)$ WITH $\beta =2$

^{38}Therefore, in the limit of $s\u21920$, we obtain

## REFERENCES

*How Nature Works: The Science of Self-Organized Criticality*

*Cambridge Nonlinear Science Series*, 2nd ed. [Cambridge University Press, Cambridge (GB), 2004].

*Mathematics of Complexity and Dynamical Systems*, edited by R. A. Meyers (Springer New York, New York, NY, 2012), pp. 463–487.

*Bursty Human Dynamics*

*Series on Complexity Science No.*, Vol. 6 (World Scientific, NJ, 2021).

*An Introduction to the Theory of Point Processes*

*Fractal-Based Point Processes*

*Gillespie Algorithms for Stochastic Multiagent Dynamics in Populations and Networks*

*Temporal Network Theory*, edited by P. Holme and J. Saramäki (Springer International Publishing, Cham, 2023), 2nd ed., pp. 165–183.

*Computational Social Sciences*, 2nd ed., edited by P. Holme and J. Saramäki (Springer International Publishing, Cham, 2023).