Detrended fluctuation analysis of earthquake data

The detrended fluctuation analysis (DFA) is extensively useful in stochastic processes to unveil the long-term correlation. Here, we apply the DFA to point processes that mimick earthquake data. The point processes are synthesized by a model similar to the Epidemic-Type Aftershock Sequence model, and we apply the DFA to time series $N(t)$ of the point processes, where $N(t)$ is the cumulative number of events up to time $t$. Crossover phenomena are found in the DFA for these time series, and extensive numerical simulations suggest that the crossover phenomena are signatures of non-stationarity in the time series. We also find that the crossover time represents a characteristic time scale of the non-stationary process embedded in the time series. Therefore, the DFA for point processes is especially useful in extracting information of non-stationary processes when time series are superpositions of stationary and non-stationary signals. Furthermore, we apply the DFA to the cumulative number $N(t)$ of real earthquakes in Japan, and we find a crossover phenomenon similar to that found for the synthesized data.


I. INTRODUCTION
Although stationarity is one of the most important properties in stochastic processes, non-stationary phenomena are rather ubiquitous in nature, ranging from disordered systems [1][2][3][4][5][6], seismicity [7][8][9][10] to biological systems [11][12][13]. Particularly in point processes, there are two typical types of non-stationary processes. The first (non-stationarity of the first type) is a process in which the probability density function (PDF) for recurrence times depends explicitly on time. Typical examples are rainfalls that exhibit daily and seasonal alterations. The other (non-stationarity of the second type) is a process where a characteristic time scale of the process, such as the mean of the interval between consecutive points, diverges. Owing to divergence, the process never reaches a steady state and thus exhibits non-stationary behaviors [2,6,14]. In this study, we focus on non-stationary processes of the first type.
Earthquakes are an example of the non-stationary point process of the first type. One of the most wellknown statistical laws of seismicity is the Rutenberg-Richter law [15], which states that the magnitude distribution of earthquakes follows an exponential distribution. This statistical law is universal in the sense that the exponential distributions are observed in any region on earth, any periods and any types of earthquakes such as mainshocks and aftershocks. However, the decay constant, the so-called b-value, depends on time [16]; thus, it is a non-stationary law. Additionally, the Omori law describes a non-stationary property for aftershocks [7], which states that the occurrence rate of aftershocks decays with the time elapsed from the mainshock. More precisely, the occurrence rate decays as a power law: λ(t) ∝ t −p for large t, where λ(t) is the occurrence rate * takuma@rs.tus.ac.jp at elapsed time t after a mainshock, and p (> 0) is a parameter. Since the rate of aftershocks λ(t) explicitly depends on t; aftershocks are intrinsically non-stationary, and earthquake occurrences are non-stationary processes of the first type.
Several methods are proposed to analyze nonstationary time series. For diffusion in heterogeneous environments, trajectories of a diffusing particle can be tracked, and using the trajectory data, the diffusion coefficient can be obtained from the time-averaged mean square displacement (MSD) calculated from the trajectory. If the process is non-stationary, the diffusion coefficient depends explicitly on the total measurement time [6,11,12,[17][18][19][20]. Thus, plotting the diffusion coefficient as a function of the measurement time provides us information on how the process ages. In another method of non-stationary data analysis of the second type, the inter-occurrence times are utilized frequently [11,21,22]. For example, the inter-occurrence-time PDFs have also been extensively used for earthquake researches [23][24][25][26][27].
Although the time-averaged MSD and interoccurrence-time PDF are useful analysis methods for non-stationary time series of the second type, these methods are not effective in case of the first type. In non-stationary processes of the first type, a longmeasurement-time limit of a time average may have a definite value, as the process does not age monotonically. Thus, the diffusion coefficient of time-averaged MSD does not depend on the total measurement time. Moreover, the inter-occurrence-time PDF analysis is based on the fact that the inter-occurrence times are independent and identically distributed (IID). However, this assumption is not valid for nonstationary processes of the first type. Therefore, it is important to develop a method to extract information of nonstationary features from the time series.
In this study, we utilize the detrended fluctuation analysis (DFA), to overcome the above-mentioned difficulties of non-stationary time-series analysis of earthquakes [28].
The difference between the DFA and the time-averaged MSD is that the local trends are subtracted in the DFA. The DFA of recurrence times of earthquakes in stationary regimes was studied to unveil the long-term correlation [29]. We note that the DFA can be utilized irrespective of the inter-occurrence times being IID. Here, we apply the DFA to both data synthesized by an earthquake model and data of real earthquakes in non-stationary regimes. In particular, we perform the DFA to a cumulative number N (t) of earthquakes occurred up to time t. As an earthquake model, we employ a simplified version of the Epidemic-Type aftershock sequence model [9].
We find crossover phenomena in DFA for both synthesized and real earthquake data. It is shown that the crossover time represents a characteristic time scale of non-stationary process. Moreover, we present analytical predictions of long time behaviors in DFA for point processes. Until now, relations between the DFA and the long-term correlation in time series have been analytically obtained for fractional Brownian motion (fBm) [30], but it is important to obtain an analytical expression for point processes as well, as these two are totally different stochastic processes [31]. The DFA has been widely used in data analysis, and thus the analytical results for point processes are also useful.
This article is organized as follows. In Sec. II, the earthquake model, which is a superposition of stationary and non-stationary point processes, is proposed, and in Sec. III, the DFA method is briefly reviewed. Sections IV and V present results of DFA for synthesized data and real earthquake data, respectively. Finally, Sec. VI presents a summary and discussion.

II. EARTHQUAKE MODEL
Here, we propose a point process describing occurrences of earthquakes over a vast area, such as the entire extent of Japan. In our model, three types of earthquakes are considered: mainshocks, aftershocks, and stationary earthquakes independent of mainshocks and aftershocks (background earthquakes).
First, we assume that mainshocks occur independently, and thus they are described by a Poisson process. This assumption is quite reasonable because the superposition of large numbers of mutually independent renewal processes becomes a Poisson process in general [32]. In fact, mainshock occurrences have been considered a Poisson process [33,34]. We have partially confirmed this assumption for a region around Japan by analyzing the Japan Meteorological Agency (JMA) catalog [35], where we define the mainshocks as earthquakes with magnitudes greater than 7. Figure 1 shows that survival probability P (τ ) of inter-occurrence times τ between successive mainshocks follows a superposition of exponential distributions. In the data, the mean inter-occurrence time is around 2.20 × 10 7 [s] ∼ = 254 [d]. The exponential distribution with mean 2.20 × 10 7 [s] also well describes the . Therefore, mainshock rate is given by 36 10 [s] in Japan.
Second, we assume that aftershocks are triggered by a mainshock. In particular, we assume the Omori law [7], which states that the rate of occurrence of aftershocks after a mainshock follows a power-law decay: where is elapsed time from a mainshock, is a degree of the aftershock rate, is a parameter characterizing a relaxation of the rate, and is a power-law exponent. For simplicity, all the parameters are fixed for all mainshocks. Third, we assume other earthquakes occur independently of mainshocks. This assumption is di erent from the Epidemic-Type Aftershock Sequence model [9], where every earthquakes are triggered by mainshocks or earthquakes triggered by mainshocks. Analyzing earthquake data in Japan from 2014 to 2015, we obtain the average occurrence rate of earthquakes. This rate is considered to be the rate of stationary earthquakes because almost all earthquakes are independent of mainshocks. Thus, we set the rate as = 1 10 [s]. Figure 2 shows a schematic view of our model.
Third, we assume other earthquakes occur independently of mainshocks. This assumption is di erent from the Epidemic-Type Aftershock Sequence model [9], where every earthquakes are triggered by mainshocks or earthquakes triggered by mainshocks. Analyzing earthquake data in Japan from 2014 to 2015, we obtain the average occurrence rate of earthquakes. This rate is considered to be the rate of stationary earthquakes because almost all earthquakes are independent of mainshocks. Thus, we set the rate as = 1 10 [s]. Figure 2 shows a schematic view of our model.

III. DETRENDED FLUCTUATION ANALYSIS
DFA were invented to analyze data which have local trends and are thus non-stationary [18]. A basic idea is to quantify fluctuations around local trends. In particular, the standard deviation is used as follows: where is a time series that we represents a local trend in [jn+1 we assume that the local trend linear function, which are obta method in [jn + 1 + 1) ]. F with when there increment of , i.e., hand, it increases as increment has a strong correlatio in the correlation. In particula there is anti-correlation in incre hand, 2 implies a positive We are interested in the DF ated by a point process. In p tonically increasing sequence, e.g number of earthquakes until tim in the original DFA, changes o ous times. Thus, in what follow that definition of function ) r DFA for point processes. In a p for sequence was also studi where inter-occurrence times ar tically distributed (IID) random occurrence times are IID rand ity density function (PDF) of th characterizes the whole process is investigated extensively in m However, earthquake time serie sequences. Thus, the inter-occur not capture non-stationary featu vantage of the DFA is that it can time series. Therefore, we use t stationary features embedded in A. earthquake d Here, we apply DFA of , i. quakes until time , using the J is enclosed within 25-50 N latit gitude with magnitude 2. inter-occurrence distribution. Therefore, the mainshock rate λ m for the entire extent of Japan is approximately given by Second, we assume that aftershocks are triggered by a mainshock. In particular, we assume that the Omori law [7], which states that the rate of occurrence of aftershocks after a mainshock follows a power-law decay: where t is the elapsed time after a mainshock, and K is the degree of the aftershock activity, c is a parameter characterizing the relaxation time of the activity, and p is the power-law exponent. IIn particular, it has been shown that the parameter p clearly depends on the magnitude of the mainshock [36]. Third, we assume that background earthquakes occur independently of mainshocks and aftershocks. This assumption differs from the Epidemic-Type Aftershock Sequence model [9], where every earthquake is triggered by mainshocks or aftershocks, which are usually triggered by a mainshock. Furthermore, we assume that the rate of the background earthquakes depends on the magnitudes of the mainshocks; that is, when the magnitude of the mainshock is large, the rate of the background earthquakes is also large. It is difficult to determine whether an earthquake is an aftershock or a background earthquake. As we consider earthquakes over a vast region, we assume that almost all earthquakes are independent of mainshocks, and thus are background earthquakes. Under this assumption, we can determine some of the model parameters. Figure 2 shows a schematic view of our model.
A Poisson process with rate λ can be generated by creating inter-event times following the exponential distribution with mean 1/λ in numerical simulations. In particular, the mainshocks and background earthquakes are generated by Poisson processes with rates λ s and λ m , respectively. As aftershocks are described by the Omori's law, we generate them by using a non-stationary Poisson process (see Appendix A for the details). Unlike a non-Markov model of aftershocks [37], it is a Markov model with the exception that it is non-stationary.

III. DETRENDED FLUCTUATION ANALYSIS
The DFA was invented to analyze data that have local trends and non-stationary features [28]. This method has been used to unravel long-range correlations in stationary as well as non-stationary time series such as heartbeat rates, weather variations, recurrence times of earthquakes, and conformation fluctuations of proteins [29,[38][39][40][41][42].
The basic idea is to quantify fluctuations around local trends as where y i is a time series that are considered, andỹ j i represents a local trend in the time interval [jn+ 1, (j + 1)n]. Thus, n is the length of these time intervals. The local trendỹ j i in the interval i ∈ [jn + 1, (j + 1)n] is given by a linear function obtained by the least-square fit to data y i in the same interval.
Even when there are local trends in data, the function F (n) characterizes a long-term correlation in the data. In particular, F (n) increases as F (n) ∝ n 1/2 when there is no correlation in increments ∆y i of y i , i.e., ∆y i ≡ y i − y i−1 . However, it increases as F (n) ∝ n α with α = 1/2 when the increment has a strong correlation, implying a power-law decay of the correlation function. In particular, α < 1/2 implies that there is an anti-correlation in increments ∆y i and α > 1/2 implies a positive correlation of increments ∆y i Here, we apply the DFA to time series y i generated by a point process. More precisely, y i is a monotonically increasing sequence defined by is the cumulative number of earthquakes up to time i. The variable i is an integer in the original DFA, whereas i in N (i) represents the continuous time; thus, it is a real number. Hence, in what follows, we use t as the argument and use the notation N (t). Note, however, that the definition of the function F (n) in Eq. (2) remains unchanged even for point processes because we only use discrete data points of N (t). In a previous study [43], the DFA for a sequence generated by a point process such as N (t) was studied, where inter-occurrence times are IID random variables. Such a process is called a renewal process. However, the inter-occurrence times may not be IID and the time series are non-stationary of the first type in earthquakes. Here, we investigate the DFA for point processes for such non-stationary time series.

IV. DETRENDED FLUCTUATIONS ANALYSIS ON SYNTHESIZED DATA
In this section, the DFA is applied to three types of data synthesized by the earthquake model: (1) Background earthquakes (Poisson processes), (2) One mainshock and its aftershocks without background earthquakes, and (3) Poissonian mainshocks with aftershocks and background earthquakes. Numerical simulations are carried out for these models and compared with theoretical predictions for small and large n. Derivations of these predictions are presented in Appendices C and D.

A. Background earthquakes (Poisson process)
First, we apply the DFA to Poisson processes, for which inter-occurrence times follow an exponential distribution with rate λ. Figure 3 shows that F (n) increases as F (n) = An 1/2 for any n > 0 and a constant A depends on rate λ of the Poisson process. A theory of the DFA for a Poisson process implies (a proof is given in Appendix C), which is confirmed using numerical simulations (inset of Fig. 3). As a Poisson process is a memory-less process, the scaling of F (n) ∝ n 1/2 is quite reasonable. However, we numerically find that scaling n 1/2 is no longer valid for renewal processes where the PDF of inter-occurrence times follows a power-law distribution with a divergent mean. Therefore, scaling n 1/2 represents the signature of a stationary Poisson process. In a biased continuous-time random walk, the variance of the displacement, which is a quantity similar to the DFA, shows an anomalous scaling [44,45].   However, definition of function F (n) remains unchanged. In a previous was also studied in renewal processes.
A. earthquake data catalog where is a time series that we are interested in and˜represents a local trend in [jn + 1) ]. For simplicity, we assume that the local trend can be described by a linear function, which are obtained by the least-square method in [jn + 1 + 1) ]. Function ) increases with when there is no correlation in the increment of , i.e., . On the other hand, it increases as with = 1 2 when the increment has a strong correlation like a power-law decay in the correlation. In particular, 2 implies that there is anti-correlation in increment . On the other hand, implies a persistent correlation of increment Here, we are interested in time series generated by a point process. In particular, is a monotonically increasing sequence, e.g., , where the number of earthquakes until time . While is discrete in the original DFA, changes of occurs at continuous times.
However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog er scaling for large . More precisely, ) increases as for n > n 0 and 10 s (= 4 5 hours). This crossover is observed for all ime series. Moreover, small-behaviors of ) are overlapped in each data ssover times are slightly changed.
hesized data a physical meaning of the crossover phenomena, we analyze data synthesized rocesses and our earthquake model. Figure 2 shows that ) increases with in Poisson processes and constant depends on rate λ of a Poisson process.

B. One mainshock and its aftershocks without background earthquakes
Second, synthesized data N (t) are generated using the earthquake model. To obtain a deeper understanding of the features of the DFA, we consider a simple situation where a mainshock occurs only once at t = 0, and there are no background earthquakes, i.e., the time series comprises one mainshock and its aftershocks.
As shown in Fig. 4, all the results of the DFA show a crossover from n 1/2 to n α scaling. For small-n behavior, F (n) shows F (n) = Bn 1/2 . By an adiabatic approximation, we approximately obtain B: where and T is the total length of the time series. For large-n behavior, F (n) also shows when p = 1 (see Appendix. D). Equation (6) is a special case of a general result for p < 3/2 which is valid for n → ∞ (see Appendix. D). Thus, the power-law exponent in the DFA is determined by p. In other words, the parameter p can be obtained from the asymptotic behavior of the DFA for N (t). This is one of the most important analytical results of our study. Figure 4 summarizes the results of the DFAs of N (t) for different parameters. Figure 4(a) shows that crossover time n c in F (n) increases with increasing parameter c and short-n behaviors are almost the same. For large-n behavior, F (n) converges to a n 1/2 scaling, which does not depend on c. As shown in Figs. 4(b) and (c), crossover time n c also depends on K and p, but the dependencies are relatively weak compared with the c dependency. Importantly, crossover time n c is thus almost proportional to c. Therefore, the parameter c can be estimated from the crossover time n c . This is significantly important when time series are a superposition of non-stationary and stationary signals, because information of the nonstationary part can be obtained without distinguishing the stationary and non-stationary time series. In Fig. 4, it is clearly shown that the asymptotic behaviors of F (n) exhibit different power-law scaling with exponent 3/2−p. Therefore, the parameter p can be obtained from the asymptotic behavior of the DFA for N (t) if background earthquakes are removed from the time series.

C. Poissonian mainshocks with aftershocks and background earthquakes
For synthesized data N (t) generated by the earthquake model with several mainshocks, we find a crossover phenomenon such that F (n) exhibits a n 1/2 to n α scaling. The superposition of two Poisson processes with rates λ 1 and λ 2 is equivalent to a Poisson process with rate λ 1 + λ 2 . Therefore, the DFA for the superposition of the two Poisson processes with rates λ 1 and λ 2 becomes F (n) = (λ 1 + λ 2 )n/ √ 15. The synthesized time series are composed of background earthquakes and aftershocks triggered by a mainshock. Since background earthquakes and mainshocks are described by Poisson processes with rates λ b and λ m , the above estimation can be utilized. For small-n behavior, F (n) shows that F (n) ≃ Cn 1/2 and C can be approximately obtained as where The small-n behavior of the DFA for N (t) is determined by λ b , λ m , and N a .

V. DETRENDED FLUCTUATION ANALYSIS ON EARTHQUAKE DATA CATALOG
Here, we apply the DFA to the cumulative number N (t) of real earthquakes included in the JMA catalog; this catalog contains data of earthquakes with magnitude M ≥ 2 and the ones that occurred in the area . On the other hand, it increases as with = 1 2 when the crement has a strong correlation like a power-law decay in the correlation. In particular, 2 implies that there is anti-correlation in increment . On the other hand, plies a persistent correlation of increment Here, we are interested in time series generated by a point process. In particular, is a onotonically increasing sequence, e.g., , where the number of earthquakes until ime . While is discrete in the original DFA, changes of occurs at continuous times.
owever, definition of function (n) remains unchanged. In a previous study, DFA for as also studied in renewal processes.
A. earthquake data catalog Figure shows DFA for earthquake time series for several periods in Japan. We find a rossover phenomenon in ); i.e., ) increases as for small and  . On the other hand, it increases as with = 1 2 when the increment has a strong correlation like a power-law decay in the correlation. In particular, 2 implies that there is anti-correlation in increment . On the other hand, implies a persistent correlation of increment Here, we are interested in time series generated by a point process. In particular, is a monotonically increasing sequence, e.g., , where the number of earthquakes until time . While is discrete in the original DFA, changes of occurs at continuous times.
However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog Figure shows DFA for earthquake time series for several periods in Japan. We find a crossover phenomenon in ); i.e., ) increases as for small and   However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function F (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function F (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function F (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog However, definition of function F (n) remains unchanged. In a previous study, DFA for was also studied in renewal processes.
A. earthquake data catalog  Fig. 5(b)]. In this synthesized data, instead of generating mainshocks according to Poisson statistics, mainshocks are assumed to occur at the same times t i (i = 1, 2, . . . ) at which real earthquakes with magnitudes greater than 7 occurred. Furthermore, we applied the DFA to N (t) for specific earthquakes such as the Tohoku and Kumamoto earthquakes. In the Tohoku earthquakes, a mainshock occurred on March 11, 2011, with a magnitude of M = 9.0 and we analyzed earthquakes after the mainshock whose area overlaps the area around the epicenter [see Fig. 5(e)]. In the Kumamoto earthquakes, a mainshock occurred on April 16, 2016, with a magnitude of M = 7.3, and we analyzes earthquakes after the mainshock whose area overlaps the area around the epicenter [see Fig. 5(e)]. In the DFAs for these real earthquakes [ Fig. 5(c)(d)], we find crossover phenomena similar to that found for the synthesized data. Moreover, we successfully generate time series N (t) with our earthquake model that reproduce the DFAs of the two earthquake time series, i.e., the Tohoku and Kumamoto earthquakes [circles in Fig. 5(c)(d)]. In the DFA of the Tohoku earthquakes at large n, however, there is a slight difference between the results of the catalog data and the earthquake model. While we assume that a mainshock occurs only at t = 0, there are a few mainshocks (earthquakes with magnitudes greater than 7) occurring after t = 0 in the catalog data. We believe that such large aftershocks significantly affect subsequent aftershocks. Therefore, our model cannot completely reproduce the DFA of Tohoku earthquakes.

VI. CONCLUSION
We found that crossover phenomena in the DFA of N (t), i.e., the number of earthquakes up to time t, are universally observed in earthquake data. Extensive numerical simulations of the earthquake model show that the crossover phenomena originate from the nonstationarity of the aftershock sequences. In particular, crossover time n c in the DFA characterizes parameter c, which represents the relaxation time of aftershocks in the Omori's law. Although we do not determine if an earthquake is an aftershock or not, we can successfully obtain information regarding the aftershocks. Therefore, our analysis is significantly important when the time series is a superposition of the two types of time series that cannot be distinguished and when one of the two types is non-stationary and the other is stationary. Moreover, we present theories of the DFA for stationary and nonstationary point processes, which are necessary for performing a thorough analysis.
In addition, if the measurement time T is transformed as T → T /c ≡T , the remaining parameters areK,T and p.
Appendix C: Theory of DFA for point process In this study, our objective is the extraction of nonstationary information from point processes using DFA. Here, we provide a theoretical argument regarding the DFA for stationary point processes with constant rate λ.
If ψ(τ ) is given by the exponential distribution the point process is referred to as the Poisson process. For the Poisson process Eqs. (C1) and (C2) is given by Note that equalities hold for the Poisson process. From the definition of σ 2 , it follows that σ 2 = λ holds for the Poisson process.

Theory of DFA for point process
Based on the stationarity of the point process, Eq. (2) can be represented as where y i = N (i) is obtained using the point process described above. We rewrite the DFA as whereỹ is defined asỹ i = y i − λi. Moreover, a and b are the coefficients of the linear fitting ofỹ i by the leastsquare method, i.e.,ỹ 0 i − λi = ai + b. By expanding the summand, we obtain The parameters a and b are given by where S xy and S 2 x are a covariance and a variance given by In Poisson processes, we have N (i)N (i + i ′ ) = N 2 (i) + N (i)[N (i ′ )−N (i)] = N 2 (i) + N (i) N (i ′ )− N (i) , because a Poisson process is a memory-less process, i.e., N (i) and N (i ′ ) − N (i) are independent. It follows that ỹ i = 0, ỹ 2 i = σ 2 i, and ỹ iỹj = σ 2 min(i, j). For general point processes, N (i) and N (i ′ ) − N (i) are not independent. It becomes