Copula-based algorithm for generating bursty time series

TL;DR

This paper introduces a copula-based algorithm to generate bursty time series with correlated interevent times, capturing non-Poissonian patterns and memory effects for more realistic modeling of dynamical processes.

Contribution

A novel copula-based method for generating correlated, heavy-tailed bursty time series, outperforming existing shuffling techniques in certain scenarios.

Findings

Successfully generates event sequences with heavy-tailed IET distributions.

Outperforms shuffling method in specific cases.

Applicable to modeling diverse dynamical processes.

Abstract

Dynamical processes in various natural and social phenomena have been described by a series of events or event sequences showing non-Poissonian, bursty temporal patterns. Temporal correlations in such bursty time series can be understood not only by heterogeneous interevent times (IETs) but also by correlations between IETs. Modeling and simulating various dynamical processes requires us to generate event sequences with a heavy-tailed IET distribution and memory effects between IETs. For this, we propose a Farlie-Gumbel-Morgenstern copula-based algorithm for generating event sequences with correlated IETs when the IET distribution and the memory coefficient between two consecutive IETs are given. We successfully apply our algorithm to the cases with heavy-tailed IET distributions. We also compare our algorithm to the existing shuffling method to find that our algorithm outperforms the shuffling method for some cases. Our copula-based algorithm is expected to be used for more realistic modeling of various dynamical processes.

Equations

$$P(\tau )\sim {\tau }^{-\alpha }.$$

$$M\equiv \frac{1}{n-1}\sum _{i=1}^{n-1}\frac{({\tau }_i-{\mu }_1)({\tau }_{i+1}-{\mu }_2)}{{\sigma }_1{\sigma }_2},$$

$$P({\tau }_i,{\tau }_{i+1})=P({\tau }_i)P({\tau }_{i+1})\left[1+rf({\tau }_i)f({\tau }_{i+1})\right],$$

$$f(\tau )\equiv 2F(\tau )-1,F(\tau )\equiv {\int }_0^{\tau }d{\tau }^{\prime }P({\tau }^{\prime }).$$

$$M\simeq \frac{⟨{\tau }_i{\tau }_{i+1}⟩-{\mu }^2}{{\sigma }^2}=\frac{r}{{\sigma }^2}{\left[{\int }_0^{\infty }d\tau \tau P(\tau )f(\tau )\right]}^2\equiv ar,$$

$$⟨{\tau }_i{\tau }_{i+1}⟩\equiv {\int }_0^{\infty }d{\tau }_i{\int }_0^{\infty }d{\tau }_{i+1}{\tau }_i{\tau }_{i+1}P({\tau }_i,{\tau }_{i+1}),$$

$$P({\tau }_{i+1}|{\tau }_i)=\frac{P({\tau }_i,{\tau }_{i+1})}{P({\tau }_i)}=P({\tau }_{i+1})\left[1+rf({\tau }_i)f({\tau }_{i+1})\right].$$

$$1-x=F({\tau }_{i+1}|{\tau }_i),$$

$$F({\tau }_{i+1}|{\tau }_i)\equiv {\int }_0^{{\tau }_{i+1}}d{\tau }^{\prime }P({\tau }^{\prime }|{\tau }_i).$$

$$c_i\equiv rf({\tau }_i)=r[2F({\tau }_i)-1],$$

$$F({\tau }_{i+1}|{\tau }_i)$$

$$=(1-c_i){\int }_0^{{\tau }_{i+1}}d{\tau }^{\prime }P({\tau }^{\prime })+2c_i{\int }_0^{{\tau }_{i+1}}d{\tau }^{\prime }P({\tau }^{\prime })F({\tau }^{\prime })$$

$$=(1-c_i)F({\tau }_{i+1})+c_{i}F({\tau }_{i+1}{)}^2.$$

$$c_i{y}^2+(1-c_i)y+x-1=0,$$

$${\tau }_{i+1}=F^{-1}\left[\frac{c_i-1+\sqrt{(c_i+1{)}^2-4c_{i}x}}{2c_i}\right],$$

$$P(\tau )$$

$$F(\tau )$$

$$r=4M.$$

$$P({\tau }_{i+1}|{\tau }_i)={\mu }^{-1}{e}^{-{\tau }_{i+1}/\mu }$$

$$\times \left[1+4M(1-2e^{-{\tau }_i/\mu })(1-2e^{-{\tau }_{i+1}/\mu })\right].$$

$${\tau }_{i+1}=\mu \ln \left[\frac{2c_i}{c_i+1-\sqrt{(c_i+1{)}^2-4c_{i}x}}\right],$$

$$c_i=rf({\tau }_i)=4M(1-2e^{-{\tau }_i/\mu }).$$

$$P(\tau )$$

$$F(\tau )$$

$$r=\frac{(2\alpha -3{)}^2}{(\alpha -1)(\alpha -3)}M.$$

$$P({\tau }_{i+1}|{\tau }_i)=(\alpha -1){\tau }_{i+1}^{-\alpha }$$

$$\times \left[1+\frac{(2\alpha -3{)}^2}{(\alpha -1)(\alpha -3)}M(1-2{\tau }_i^{1-\alpha })(1-2{\tau }_{i+1}^{1-\alpha })\right].$$

$${\tau }_{i+1}={\left[\frac{2c_i}{c_i+1-\sqrt{(c_i+1{)}^2-4c_{i}x}}\right]}^{1/(\alpha -1)}.$$

$$c_i=rf({\tau }_i)=\frac{(2\alpha -3{)}^2}{(\alpha -1)(\alpha -3)}M(1-2{\tau }_i^{1-\alpha }).$$

$$P(\tau )$$

$$F(\tau )$$

$$P({\tau }_{i+1}|{\tau }_i)=\frac{{\tau }_c^{\alpha -1}}{\Gamma \left(1-\alpha ,1/{\tau }_c\right)}{\tau }_{i+1}^{-\alpha }{e}^{-{\tau }_{i+1}/{\tau }_c}$$

$$\times \left\{1+\frac{M}{a}\left[1-2\frac{\Gamma \left(1-\alpha ,{\tau }_i/{\tau }_c\right)}{\Gamma \left(1-\alpha ,1/{\tau }_c\right)}\right]\left[1-2\frac{\Gamma \left(1-\alpha ,{\tau }_{i+1}/{\tau }_c\right)}{\Gamma \left(1-\alpha ,1/{\tau }_c\right)}\right]\right\}.$$

$$\frac{\Gamma (1-\alpha ,{\tau }_{i+1}/{\tau }_c)}{\Gamma (1-\alpha ,1/{\tau }_c)}=\frac{c_i+1-\sqrt{(c_i+1{)}^2-4c_{i}x}}{2c_i},$$

Full text

I Introduction

Dynamical processes in various natural and social phenomena have been described by a series of events or event sequences showing non-Poissonian, bursty temporal patterns. Examples include solar flares Wheatland et al. (1998), earthquakes Corral (2004); de Arcangelis et al. (2006), neuronal firings Kemuriyama et al. (2010), and human activities Barabási (2005); Karsai et al. (2018). Temporal correlations in such bursty time series can be understood not only by heterogeneous interevent times (IETs) but also by correlations between IETs Goh and Barabási (2008); Jo (2017). Here the IET, denoted by $\tau$, is defined as a time interval between two consecutive events. On the one hand, heterogeneities of IETs have been characterized by heavy-tailed or power-law IET distributions $P(\tau)$ with a power-law exponent $\alpha$ Karsai et al. (2018):

[TABLE]

On the other hand, the correlations between IETs can be characterized by a memory coefficient Goh and Barabási (2008) among other measures such as local variation Shinomoto et al. (2003) and a bursty train Karsai et al. (2012). The memory coefficient $M$ is defined as the Pearson correlation coefficient between two consecutive IETs, whose value for a sequence of $n$ IETs, i.e., $\{\tau_{i}\}_{i=1,\cdots,n}$, can be estimated by

[TABLE]

where $\mu_{1}$ ($\mu_{2}$) and $\sigma_{1}$ ($\sigma_{2}$) are the average and the standard deviation of the first (last) $n-1$ IETs, respectively. Positive $M$ implies that the small (large) IETs tend to be followed by small (large) IETs. Negative $M$ implies the opposite tendency, while $M=0$ is for the uncorrelated IETs. We mainly focus on the case with positive $M$, based on the empirical observations Goh and Barabási (2008); Wang et al. (2015); Guo et al. (2017); Böttcher et al. (2017).

The dynamical processes, such as spreading and diffusion, taking place in a network of individuals are known to be strongly affected by bursty interaction patterns between individuals Vazquez (2007); Karsai et al. (2011); Miritello et al. (2011); Rocha et al. (2011); Jo et al. (2014); Perotti et al. (2014); Delvenne et al. (2015); Artime et al. (2017); Hiraoka and Jo (2018). This topic has been studied in a framework of temporal networks Holme and Saramäki (2012); Masuda and Lambiotte (2016), where a link connecting two nodes is considered being existent or activated only at the moment of interaction. For example, an infectious disease, information, or a random walker on a node, say $A$, can be transferred to the $A$’s neighboring node, say $B$, only when $A$ interacts with $B$. Therefore, the dynamical behavior, such as spreading speed, would be influenced by how often these nodes interact with each other, more generally, by the temporal interaction pattern between them. This temporal interaction pattern can be described by the event sequence when each event denotes the interaction. Consequently, modeling and simulating realistic temporal networks often requires us to generate event sequences for links (or nodes) showing the empirically observed properties, such as heavy-tailed IET distribution in Eq. (1) and/or memory effects measured by the memory coefficient in Eq. (2).

It is straightforward to generate event sequences characterized only by the IET distribution $P(\tau)$, i.e., when IETs are fully uncorrelated with each other as in the renewal process Mainardi et al. (2007): $n$ IETs are independently drawn from the given $P(\tau)$ to make a sequence of IETs, $\{\tau_{i}\}_{i=1,\cdots,n}$. Provided that the first event occurs at time $t_{0}$, the events are considered to occur at times $t_{i}=t_{0}+\sum_{i^{\prime}=1}^{i}\tau_{i^{\prime}}$ for $i=1,\cdots,n$.

In contrast, the generative methods for event sequences with correlated IETs have not been fully explored, except for a few recent works. We focus on two existing methods: The first method is to shuffle a set of IETs, prepared using a given $P(\tau)$, to implement the desired value of $M$ in Eq. (2) Hiraoka and Jo (2018). By this method one can control $P(\tau)$ and $M$ independently, while the shape of $P(\tau)$ may set bounds on the range of $M$ Guo et al. (2017). Note that for the implementation of this method, the number of IETs should be predetermined, which however is not necessarily the case. For example, let us consider the simulation of a dynamical process using the event sequence; it is often the case that how many IETs are to be needed for reaching a stationary state of the process may not be known a priori. The second method is to generate IETs sequentially using the conditional probability distribution $P(\tau|\tau^{\prime})$ Artime et al. (2017), by which a next IET $\tau$ is drawn given the previous IET $\tau^{\prime}$. One can generate an arbitrary number of IETs without predetermining the number of IETs. However, since the function used in Ref. Artime et al. (2017) is $P(\tau|\tau^{\prime})=\alpha(1+\tau^{\prime})^{\alpha}/(\tau+\tau^{\prime})^{\alpha+1}$, $\alpha$ controls both the heterogeneity of IETs and the degree of memory effects, implying the inevitable dependency between the statistics of IETs and correlations between IETs.

To overcome disadvantages in the previous methods, we propose an alternative generative method using the conditional probability distribution $P(\tau|\tau^{\prime})$ but based on the Farlie-Gumbel-Morgenstern copula Nelsen (2006); Takeuchi (2010), by which one can generate an arbitrary number of IETs without predetermining the number of IETs, while the statistics of IETs and their correlations can be independently controlled. In Sec. II we introduce the copula-based algorithm for generating event sequences with correlated IETs. In Sec. III we apply our algorithm to the cases with exponential and power-law IET distributions as well as power-law IET distribution with exponential cutoff, and we also discuss the performance of our algorithm for different cases in terms of computational times, in comparison to the shuffling method used in Ref. Hiraoka and Jo (2018). Finally, we conclude our work in Sec. IV.

II Copula-based algorithm

We introduce the copula-based algorithm for generating event sequences with correlated interevent times (IETs) for a given IET distribution $P(\tau)$ and memory coefficient $M$ between two consecutive IETs. For this, we model the joint probability distribution $P(\tau_{i},\tau_{i+1})$ by adopting a Farlie-Gumbel-Morgenstern (FGM) copula among others Takeuchi (2010); Nelsen (2006). It is because the FGM copula is simple and analytically tractable, despite the range of correlation being somewhat limited, which will be discussed below. The joint probability distribution based on the FGM copula is written as 111The FGM copula is a function $C$ joining a bivariate cumulative distribution function (CDF) to their one-dimensional marginal CDFs such that $G(x,y)=C[u(x),v(y)]=uv[1+r(1-u)(1-v)]$, where $u$ and $v$ are CDFs of variables $x$ and $y$, respectively Takeuchi (2010); Nelsen (2006). Then the bivariate probability distribution function (PDF) of $x$ and $y$ is obtained by $\frac{\partial^{2}G}{\partial x\partial y}=P_{1}(x)P_{2}(y)[1+r(2u-1)(2v-1)]$, where $P_{1}(x)$ and $P_{2}(y)$ denote PDFs of $x$ and $y$, respectively. The FGM copula for IETs has recently been found to be useful for the analytical approach to autocorrelation functions Jo (2019).

[TABLE]

where

[TABLE]

Here $F(\tau)$ denotes the cumulative distribution function (CDF) of $P(\tau)$, and $P(\tau_{i})$ and $P(\tau_{i+1})$ are assumed to have the same functional form. The parameter $r$, controlling the correlation between two consecutive IETs, is in the range of $-1\leq r\leq 1$ because $P(\tau_{i},\tau_{i+1})\geq 0$ and $0\leq F(\tau)\leq 1$, hence $-1\leq f(\tau)\leq 1$. It is straightforward to relate $r$ with the memory coefficient $M$ in Eq. (2) using the FGM copula in Eq. (3) as follows:

[TABLE]

where

[TABLE]

and $\mu$ and $\sigma$ are the mean and standard deviation of IETs, respectively. The positive constant $a$ is determined only by $P(\tau)$, irrespective of $M$. The upper bound of $a$ is $1/3$ for any $P(\tau)$ as proven in Ref. Schucany et al. (1978), i.e., $|M|\leq a=1/3$. Despite such relatively weak correlations by the FGM copula, one can still study the FGM copula as the empirical values of $|M|$ tend to be relatively small in several cases Goh and Barabási (2008); Wang et al. (2015); Guo et al. (2017).

Using the joint probability distribution in Eq. (3), one can get the conditional probability distribution of $\tau_{i+1}$ for a given $\tau_{i}$:

[TABLE]

In order to draw a value of $\tau_{i+1}$ from $P(\tau_{i+1}|\tau_{i})$, we use the inverse transform sampling or transformation method Clauset et al. (2009): For a random number $x$ drawn from a uniform distribution defined in $0\leq x<1$, we obtain $\tau_{i+1}$ for a given $\tau_{i}$ by solving the equation

[TABLE]

where

[TABLE]

For convenience, we define

[TABLE]

to rewrite Eq. (9) as

[TABLE]

Then from Eq. (8) with $y\equiv F(\tau_{i+1})$, one gets

[TABLE]

leading to

[TABLE]

where $F^{-1}$ denotes the inverse function of $F(\tau)$. The sign of the square root was chosen to be “+” for satisfying $y=1$ for $x=0$ and $y\to 0$ for $x\to 1$.

By our copula-based algorithm, one can generate the event sequence without predetermining the number of IETs as well as with $P(\tau)$ and $M$ being controlled independently. The generating procedure is as depicted in Fig. 1: We randomly draw the first IET from $P(\tau)$, denoted by $\tau_{1}$. Then $\tau_{i}$ is used to generate $\tau_{i+1}$ using $P(\tau_{i+1}|\tau_{i})$ or Eq. (13) with a random number $x$ independently drawn for each $i=1,2,\cdots$. This algorithm can be called Markovian in the sense that $\tau_{i+1}$ depends only on $\tau_{i}$ but not on $\tau_{i^{\prime}}$ with $i^{\prime}<i$.

III Results

We apply the copula-based algorithm to three well-known interevent time (IET) distributions, i.e., exponential and power-law IET distributions as well as power-law IET distribution with exponential cutoff. We also test the performance of our algorithm for these three cases, in comparison to the shuffling method Hiraoka and Jo (2018).

III.1 Exponential IET distribution

As a simple example, we consider the exponential IET distribution with the mean IET $\mu$:

[TABLE]

from which we get

[TABLE]

This implies that $|M|\leq 1/4$ by Eq. (5). The conditional probability distribution in Eq. (7) is written as

[TABLE]

Then from Eq. (13) we obtain $\tau_{i+1}$ for a given $\tau_{i}$ and a random number $x$ as

[TABLE]

where

[TABLE]

For the demonstration, we generate $100$ event sequences of size $n=10^{5}$ for $\mu=100$ and $M=0.1$ to analyze them by measuring the conditional probability distributions and the resultant IET distribution, which are comparable to the corresponding analytical curves, as shown in Fig. 2(a,b). We also measure the memory coefficient $M=0.100(4)$ from the same generated event sequences, which is in good agreement with the input value of $0.1$.

III.2 Power-law IET distribution

Next, we study the case with a power-law IET distribution with power-law exponent $\alpha$:

[TABLE]

where $\theta(\cdot)$ denotes the Heaviside step function and the lower bound of IET has been set to be $1$. We assume that $\alpha>3$ for the finite variance of IETs. The relation between $r$ and $M$ is obtained as

[TABLE]

The conditional probability distribution in Eq. (7) for $\tau_{i}\geq 1$ and $\tau_{i+1}\geq 1$ is written as

[TABLE]

Then from Eq. (13), we obtain $\tau_{i+1}$ for a given $\tau_{i}$ and a random number $x$ as

[TABLE]

where

[TABLE]

Using our copula-based algorithm we generate $100$ event sequences of size $n=10^{5}$ for $\alpha=3.5$ and $M=0.07$. Note that $\alpha=3.5$ leads to $|M|\leq a=5/64\approx 0.078$. From the generated event sequences the conditional probability distributions and the resultant IET distribution are measured and compared to the corresponding analytical curves, as shown in Fig. 2(c,d). We also find the measured memory coefficient $M=0.08(1)$ comparable to the input value of $0.07$.

III.3 Power-law IET distribution with exponential cutoff

As a more realistic case evidenced by empirical results Karsai et al. (2018), we consider a power-law IET distribution with exponential cutoff $\tau_{c}$ as

[TABLE]

where $\Gamma(\cdot,\cdot)$ is the upper incomplete Gamma function and the lower bound of IET has been set to be $1$. Note that the IET distribution in Eq. (26) reduces to that in Eq. (20) in the limit of $\tau_{c}\to\infty$, or to that in Eq. (14) if $\alpha=0$ and the lower bound of IET is set to be [math]. In contrast to the exponential and power-law cases, in the case of power law with exponential cutoff the analytic calculation of the integration in Eq. (5) is not straightforward, hence $a$ in Eq. (5) will be numerically evaluated. With this numerical value of $a$, the conditional probability distribution in Eq. (7) is written as

[TABLE]

Then the value of $\tau_{i+1}$ can be obtained for a given $\tau_{i}\geq 1$ and a random number $x$ by numerically solving Eq. (13), namely,

[TABLE]

where

[TABLE]

We generate $100$ event sequences of size $n=10^{5}$ for $\alpha=2.1$, $\tau_{c}=10^{3}$, and $M=0.015$. Note that for $\alpha=2.1$ and $\tau_{c}=10^{3}$, we have $|M|\leq a\approx 0.020$. The generated event sequences are analyzed to result in the conditional probability distributions and the resultant IET distribution that are comparable to the corresponding analytical curves, as shown in Fig. 2(e,f). The measured memory coefficient $M=0.015(5)$ is found to be close to the input value of $0.015$. All these results indicate that the correlations between two consecutive IETs are successfully implemented.

III.4 Computation times

We discuss the performance of our algorithm in terms of computation times for generating event sequences with correlated IETs. For the exponential and power-law IET distributions, the generation of event sequences is easy and fast, while for the power-law IET distribution with exponential cutoff (“power+cutoff” in short) one needs to numerically solve Eq. (29) for the generation of each IET, implying longer computation times than other cases. To study this issue, we measure the average computation times in seconds for generating event sequences. We use codes written in C on a Linux system with 3.5 GHz Intel Core i5-7600 CPU and 16 GB RAM. We also use GNU Scientific Library for calculating the incomplete Gamma function, e.g., in Eqs. (29) and (30).

We test various combinations of parameter values for the estimation of computation times when $n=10^{5}$ is fixed. For the exponential case, we use $\mu=10,10^{2},10^{3}$ and $M=0.01,0.03,0.1,0.2$. For the power-law case, we use $\alpha=3.3,3.5,4$ and $M=0.001,0.01,0.03,0.07$. For the power+cutoff case, $\alpha=2.1,2.5,2.9$ and $M=0.005,0.01,0.015$ are used for a fixed $\tau_{c}=10^{3}$. Figure 3(a–c) shows that the average computation times in the power+cutoff case are larger than those for other two cases, as expected. We observe that the computation times are linearly increasing with $n$ in Fig. 3(d–f), which is trivial for the copula-based algorithm.

We now compare the computation times of the copula-based algorithm to those of the shuffling method in Ref. Hiraoka and Jo (2018). For implementing the shuffling method, we first draw $n$ random values from $P(\tau)$ to make an IET sequence $T\equiv\{\tau_{i}\}_{i=1,\cdots,n}$. Using the definition of Eq. (2), we measure the memory coefficient from $T$, denoted by $\tilde{M}$. Two IETs are randomly chosen in $T$ and swapped only when this swapping makes $\tilde{M}$ closer to $M$, i.e., the target value. By repeating the swapping, we can obtain the IET sequence whose $\tilde{M}$ is close enough to $M$. Precisely, the swapping stops when $|M-\tilde{M}|<\epsilon$ with a small number $\epsilon$.

The computation times for generating event sequences by the shuffling method are estimated for the same combinations of parameter values as in the copula-based algorithm, together with several values of $\epsilon=0.001,\cdots,0.02$. Note that the computation time for the shuffling method is the sum of the time for generating the initial set of IETs for each given IET distribution and the time for the shuffling procedure. In Fig. 3(a–c) we find that the copula-based algorithm outperforms the shuffling method only for the power+cutoff case by a factor of $10$, while for the exponential and power-law cases the shuffling method is around twice faster than the copula-based algorithm. Therefore, the shuffling procedure itself seems relatively fast. Then much slower generation by the shuffling method in the power+cutoff case could be due to the longer time for generating the initial set of IETs using Eq. (26) 222To generate the IETs from the power-law distribution with exponential cutoff, we have used the method described in Ref. Clauset et al. (2009), which indeed requires one more random number per each IET than the cases with exponential and power-law IET distributions.. We also find that the computation times for the shuffling method are mostly insensitive to the variation of parameter values, except for the effect of $M$ in the exponential and power-law cases. This can be explained by the fact that the initial sequence of IETs drawn from $P(\tau)$ is uncorrelated, i.e., $\tilde{M}\approx 0$, implying that it takes longer time to reach the larger target value of $M$. Such $M$-dependence is not apparent in the power+cutoff case probably due to the dominant effects of the time for generating the initial set of IETs. In addition, the computation times turn out to be linearly increasing with $n$ in Fig. 3(d–f) as in the copula-based algorithm.

Finally the quality of event sequences generated by the shuffling method is tested by measuring the conditional probability distribution $P(\tau_{i+1}|\tau_{i})$ and the IET distribution $P(\tau)$ from the generated event sequences using the same parameter values as in Fig. 2. Here we choose $\epsilon=10^{-3}$ which is of the order of the standard deviations of estimated $M$ in the case with the copula-based algorithm. In Fig. 4 we find almost the same behaviors as in the results from the copula-based algorithm, implying that the copula-based algorithm and the shuffling method can be used interchangeably. It might be due to the fact that both methods are designed to implement only the correlations between two consecutive IETs for a given IET distribution, while randomizing or ignoring all other higher-order correlations between IETs. In this sense, another generative method called the Laplace Gillespie algorithm (LGA) Masuda and Rocha (2018) can also be compared to our algorithm and shuffling method as it can generate event sequences with correlated IETs. However, we leave this comparison as a future work mainly because fine-tuning the value of $M$ by the LGA seems to be difficult.

IV Conclusion

We have proposed the Farlie-Gumbel-Morgenstern (FGM) copula-based algorithm for generating event sequences with correlated interevent times (IETs) for a given IET distribution and a given memory coefficient between two consecutive IETs. This is to overcome the disadvantages in the previous generative methods, i.e., the shuffling method in Ref. Hiraoka and Jo (2018) and the method using conditional probability distribution in Ref. Artime et al. (2017): By adopting the conditional probability distribution $P(\tau_{i+1}|\tau_{i})$ based on the FGM copula Takeuchi (2010); Nelsen (2006), one can generate an arbitrary number of IETs without predetermining the number of IETs, while the statistics of IETs and their correlations can be independently controlled. After deriving the analytical forms of the next IET $\tau_{i+1}$ for a given previous IET $\tau_{i}$, we show that our algorithm successfully generates the event sequences with desired statistical properties.

We also compare the performance of our copula-based algorithm to the shuffling method in terms of computation times: It turns out that the copula-based algorithm outperforms the shuffling method for generating event sequences with power-law IET distributions with exponential cutoff, while for the exponential and power-law IET distributions the shuffling method is around twice faster than the copula-based algorithm. We note that since both methods generate event sequences of the same statistical properties, any of them can be used appropriately. Considering the advantages of our algorithm, we expect our algorithm to be used for modeling and simulating more realistic event sequences, eventually for more realistic temporal networks.

Finally we remark on the limited range of $M$. Apart from the bounds of $M$ set by the shape of $P(\tau)$ Guo et al. (2017), the range of $M$ is also bounded due to the form of the FGM copula, while it is more flexible in the shuffling method. Therefore, other members of the FGM family, e.g., the iterated FGM copula Huang and Kotz (1984), can be investigated to explore a wider range of $M$ for the IET distributions of our interest as a future work. In addition, our copula-based algorithm can be used to generate any other sequence of correlated variables.

Acknowledgements.

H.-H.J. was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2018R1D1A1A09081919). W.-S.J. was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2016R1D1A1B03932590).