Brief Communication: Earthquake sequencing: analysis of time series constructed from the Markov chain model

Directed graph representation of a Markov chain model to study global earthquake sequencing leads to a time series of state-to-state transition probabilities that includes the spatio-temporally linked recurrent events in the recordbreaking sense. A state refers to a configuration comprised of zones with either the occurrence or non-occurrence of an earthquake in each zone in a pre-determined time interval. Since the time series is derived from non-linear and nonstationary earthquake sequencing, we use known analysis methods to glean new information. We apply decomposition procedures such as ensemble empirical mode decomposition (EEMD) to study the state-to-state fluctuations in each of the intrinsic mode functions. We subject the intrinsic mode functions, derived from the time series using the EEMD, to a detailed analysis to draw information content of the time series. Also, we investigate the influence of random noise on the data-driven state-to-state transition probabilities. We consider a second aspect of earthquake sequencing that is closely tied to its time-correlative behaviour. Here, we extend the Fano factor and Allan factor analysis to the time series of state-to-state transition frequencies of a Markov chain. Our results support not only the usefulness of the intrinsic mode functions in understanding the time series but also the presence of power-law behaviour exemplified by the Fano factor and the Allan factor.


Introduction
Earthquake sequencing has been the subject of detailed research (Nava et al., 2005;Ünal and Çelebioglu, 2011;Ünal et al., 2014;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Telesca and Lovallo, 2008;Cavers andVasudevan, 2013, 2015;Vasudevan andCavers, 2012, 2013) both in the regional and global sense in recent years.Nava et al. (2005) have introduced the Markov chain model to study the earthquake sequencing in a seismogenically active region where the region is partitioned into zones.The functionality of the method is determined by the characteristics of the state-to-state transitions where each state is described by the earthquake occupancy of the zones.In particular, for a given number of zones, N, a state corresponding to a time interval is expressed as a concatenation of binary digits b N−1 . ..b 1 b 0 , where b L = 1 (or b L = 0) indicates there was (or was not) an earthquake occurrence in zone L during the specified time interval.Thus, states can fall into zones of no occupancy to full occupancy at the extreme and into zones where some are occupied and some are not.The approach of Nava et al. (2005) was immediately extended to other regions (Herrera et al., 2006;Ünal and Çelebioglu, 2011;Ünal et al., 2014).Cavers and Vasudevan (2013) adapted the method of Nava et al. (2005) to a global catalogue which was partitioned into zones on the basis of the tectonic boundaries (DeMets et al., 1990(DeMets et al., , 2010;;Bird, 2003;Kagan et al., 2010).The existing Markov chain model was refined by incorporating the record-breaking recurring events for each event in the catalogue under certain constraints.A directed graph representation of the modified Markov chain model was then subjected to detailed analysis for forecasting purposes (Cavers and Vasudevan, 2015).
One consequence of the approach taken by Cavers and Vasudevan (2015) and Vasudevan and Cavers (2013) is that it results in a time series of state-to-state transition frequencies of the modified Markov chain model, x sstf (t).This time series is for an optimized time interval, t.The fluctuations in state-to-state transitions are t sampled.The time series is a comprehensive representation of earthquake sequencing in which interaction of seismic events within and among zones are considered.Therefore, it can be subjected to a detailed analysis.
Earthquake sequencing may be considered a non-linear and non-stationary process (Kanamori, 2003;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Telesca and Lovallo, 2008;Flores-Marquez and Valverde-Esparza, 2012).In earthquake sequencing, earthquakes are viewed as part of a point process, with earthquake events occurring at some random locations in time.This means that the earthquake sequencing is dictated by the set of event times, and can also be expressed by the set of time intervals between events.The time series of earthquakes for any time interval can be analysed in many ways (Bohnenstiehl et al., 2001;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Telesca and Lovallo, 2008).
We postulate here that the non-linear and non-stationary behaviour in the time series should also be present in the time series of the state-to-state transition frequencies derived from earthquake sequencing.Hence, we consider the approaches of Telesca et al. (2001Telesca et al. ( , 2009Telesca et al. ( , 2011) ) and Telesca and Lovallo (2008) to be appropriate for a study here.
Non-linear and non-stationary time series have been examined in recent years with a method known as empirical mode decomposition (EMD) and the intrinsic mode functions derived from this are useful in this regard (Huang et al., 1998).The present time series of state-to-state transition frequencies is suited for such a study.
In general, the time series has non-zero amplitudes for the state-to-state transition frequencies (Cavers and Vasudevan, 2015).In this particular case, there are instances where there are no earthquakes exceeding the magnitude of 5.6 in all zones for one or more time steps.This introduces "intermittency" in the time series.
However, because of the presence of intermittency in it, an ensemble approach to empirical mode decomposition, EEMD (Wu andHuang, 2004, 2009;Flandrin et al., 2004Flandrin et al., , 2005) ) is applied here.The intermittency problem is handled with the addition of random noise to the time series before carrying out the EEMD (Wu and Huang, 2009).We examine the criteria used for the selection of the added noise and the ensemble number for the EEMD.
Another aspect of the study here is to ask a question if the time series resulting from a directed graph representation of the Markov chain model of earthquake sequences exhibits power-law statistics similar to a description of fractal stochastic point processes (Telesca et al., 2001(Telesca et al., , 2009(Telesca et al., , 2011) ) to model the time-occurrence sequence of seismic events.Quantifying the earthquake sequencing in terms of its fractal properties was done by means of the Fano factor and the Allan factor (Allan, 1966;Barnes and Allan, 1966;Lowen andTeich, 1993, 1995;Thurner et al., 1997;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012;Serinaldi and Kilsby, 2013).Since the fractal properties of the time series studied here has never been investigated, we calculate the Fano factor and the Allan factor for the purpose of quantitative analysis.
The remainder of the paper is divided into three sections.In the next section, we show how the time series of the stateto-state transition frequencies for a modified Markov chain model as described in Cavers and Vasudevan (2015) is generated.In the following section, we describe the EEMD procedure used and the analysis of the results that accrue from this procedure.We extend the approaches of Telesca et al. (2001Telesca et al. ( , 2009Telesca et al. ( , 2011) ) and Telesca and Lovallo (2008) to calculate the Fano factor and the Allan factor with a view to study the fractal properties of the time series.In the last section, we discuss the results of the analysis methods and draw certain inferences about the state-to-state transition frequencies.
To build a Markov chain model we first partition the region, either local or global, into zones.Typically these zones are made up of rectangles that divide the region (Nava et al., 2005;Ünal and Çelebioglu, 2011).Recently, other partitions have been used.In particular, Cavers and Vasudevan (2015) used a simplified five-zone plate boundary template as given by Kagan et al. (2010) to study global seismicity, while Ünal et al. (2014) used a seismic zones map that uses geographic information system analysis to divide Turkey into regions.For this particular study, we used the five-zone model described in Cavers and Vasudevan (2015) and give an overview of its construction here.Kagan et al. (2010) partitioned the shallow (≤ 70 km depth) events with moment magnitude, M w > 5.6 from the Global CMT catalogue (1 January 1982 -31 March 2008) into five-zone sub-catalogues using their grid-assignment schemes (Table 1).The selected catalogue consists of 6752 earthquakes with 4407 from Zone 4 (Trenches), 723 from Zone 3 (Fast-spreading ridges), 487 from Zone 2 (Slowspreading ridges), 898 from Zone 1 (Active continent), and 237 from Zone 0 (Plate interior) respectively.For these five zones, we express a state, corresponding to a time interval t, as a concatenation of binary digits b 4 b 3 b 2 b 1 b 0 , where b L = 1 indicates an earthquake occurrence in zone L during the specified time interval t, and b L = 0 indicates the lack of an earthquake occurrence in zone L during the specified time interval t.We use = [θ ij ] to denote the transition frequency matrix, where θ ij is the number of occurrences of transitions from state i to state j .Letting s(n) represent the A finite-state Markov chain can be depicted using a digraph representation, G, where the set of possible states (binary strings of length 5) are the nodes, and an arc (i, j ) connects two states i and j if and only if p ij > 0 (Jarvis and Shier, 1996).Figure 1 shows an example of a digraph representing a Markov chain with a three zone partition, hence, there are 2 3 = 8 states, {000, 001, 010, 011, 100, 101, 110, 111}, that we write in decimal format {0, 1, 2, 3, 4, 5, 6, 7}, respectively.In this figure, we do not show all of the possible transitions between states and typically an arc (i, j ) is omitted when p ij = 0. We follow the same decimal state labelling format as in Fig. 1 for our 2 5 = 32 states, that is, state "0" (representing 00000 in binary) corresponds to no earthquake occurrence in all five zones in the chosen time interval, t, and state "31" (representing 11111 in binary) points to earthquake occurrences in all five zones.Table 2 shows details for defining all other states, "1" to "29".
For a Markov chain structure given earlier for the five zones, the computation of transition frequencies and hence, transition probabilities, depend on the chosen time interval, t.We use the simple rules outlined by Nava et al. (2005) to choose t: 1. t should be small enough such that the hazard estimations are useful; Table 2. Zone and state definition used in the construction of a directed graph of a Markov chain."0" and "1" refer to the no occurrence or occurrence of an earthquake for a given zone.For five zones, there are 32 states.
State Zone 4 Zone 3 Zone 2 Zone 1 Zone 0 2. t should not be too small that the most frequently occurring transition is from state 0 to state 0; 3. t should not be too large that state 31 to state 31 transitions are dominant.
So, for the threshold magnitudes chosen, t should be large enough to allow interaction among regions and make esti-mates of Markov chain transition probabilities robust.Following the selection rules given elsewhere (Nava et al., 2005;Ünal and Çelebioglu, 2011;Cavers and Vasudevan, 2015), we used a t value of 9 days for the construction of the Markov chain of transition probabilities.The combinatorial structure of a digraph representation of the Markov chain model contains important information for earthquake sequencing (Cavers and Vasudevan, 2015).It is often useful to use a weight, w ij , for each arc (i, j ) of the digraph to get a weighted digraph.The weights have the form w ij = θ ij , w ij = p ij , or can be empirically derived from the Markov chain.To introduce spatial-temporal complexity into the model so that transitions with earthquake occurrences at large distances have less of an impact on our model than transitions with earthquake occurrences at short distances, we follow the approach by Cavers and Vasudevan (2015) to modify the weights w ij in the weighted digraph by considering recurrences.Each earthquake (event) in a zone may have several recurring events in the record-breaking sense (Davidsen et al., 2008).For example, an event j is treated as a record with respect to an earthquake i if no event takes place within the spatial distance, d ij , between i and j around i during the time interval [t i , t j ] with t i < t j .The next recordbreaking event, k, in the catalogue with reference to the original event, i, during the time interval [t i , t k ] with t i < t k will have a spatial distance, d ik , less than d ij .The recurring events for one event in a given zone may fall into other zones or may be in the same zone.This flexibility adds to the possibility of interactions among zones.We first form the network of recurrences as described by Davidsen et al. (2008).The weight applied to each arc in the network of recurrences is derived empirically by using a total count of record breaking events between the corresponding earthquake zones and the distance involved (Cavers and Vasudevan, 2015;Vasudevan and Cavers, 2013).Each recurrence from an earthquake a to an earthquake b in the sequence is given a weight between 0 and 1, with a weight equal to 1 if the distance between a and b is less than 50 km.If the distance is r with r > 50 km and earthquakes a and b occur in Zones j and k respectively, a weight of is given, where L j k (r) defined by Cavers and Vasudevan (2015) is the number of record-breaking events from zone j to zone k at distance at most r in the network of recurrences.
The function in Eq. ( 3) is a decreasing function in r giving a weight close to 0 when the distance r is large.Note that for r = 50 km, an output of 1 is given while for r = 20 000 km, an output of 0 is given.As described by Cavers and Vasudevan (2015), a Markov chain with the inclusion of spatiotemporal complexity of recurring events is derived by summing the weights of the recurrence arcs corresponding to occurrences from state i to state j in consecutive time intervals.Here, we calculated the time series of the resulting stateto-state sequence (Fig. 2a) and the corresponding transition frequency matrix (Fig. 2b).There is one comment in order here.Figures 2a and b provide different representations of the same Markov chain.The first can be considered "dynamic", because it shows the time evolution of the transition from one state to another in consecutive time intervals of 9 days each.The second can be considered "static" because it shows the transition probabilities from one state to another but considering the whole earthquake sequence occurred during the whole observation period.However, they are not equivalent.We can go from the time series data to transition-frequency matrix.We cannot go from the transition-frequency matrix to time series without the additional information such as the catalogue and the record-breaking statistics of recurrences.
Since it is obtained from the non-linear, non-stationary global earthquake sequence, we consider it non-linear and nonstationary as well, and hence, can be subjected to analysis methods.Although it is not shown here, the approach equally applies to earthquake catalogues from localized seismogenic zones.

Analysis methods and results
Each sample in the time series shown in Fig. 2a represents a "zone-configuration" state (Table 2).By definition, a zoneconfiguration has no zone or some zones or all zones highlighted by an earthquake or more in the optimally chosen time interval.Going from one sample to the next does not only represent going from one state to the next but also shows the amplitude fluctuation between them.The adjacent states could represent the same zone-configuration or different zone-configurations.The time series deduced from using the present approach with the five-zones marks the state-tostate fluctuations arising out of the fluctuations of oscillations or earthquake occurrences in the five-zones.We present in the following two analysis methods to glean an insight into the characteristics of the time series.

Ensemble empirical mode decomposition as applied to state-to-state transition frequency sequence
For non-linear and non-stationary time series, the method of empirical mode decomposition (EMD) has been recently proposed as an adaptive time-frequency analysis method (Huang et al., 1998(Huang et al., , 1999) ) to decompose the original data into a basis set of intrinsic mode functions.Since the process that leads to the state-to-state transition frequency sequence or time series is inherently non-linear and non-stationary, it is appropriate to apply the EMD to this data to understand the behaviour of the intrinsic mode functions.The time series (Fig. 2a) reveals the fluctuations in the state-to-state transition frequencies arising out of varying occupancy of the zones from one time interval to the next.A situation would easily arise when two or three successive state-to-state transitions do not have earthquake occurrences in any of the zones studied.This would translate into intermittency in the time series.Recent studies (Flandrin et al., 2004(Flandrin et al., , 2005;;Gledhill, 2003;Wu andHuang, 2004, 2009) support the idea of carrying out noise-added analyses with the EMD.The noiseadded analyses involves multiple realizations of added noises to the time series in question, leading to the ensemble EMD (EEMD), as proposed by Wu andHuang (2004, 2009).
In the EEMD, the signal or the time series in question with the added Gaussian white noise, denoted as one trial, would populate the whole time-frequency space uniformly with the constituting component of different scales.Since the noise added in each trial is different, the ensemble mean of the noise cancels out and, hence, the signal resides in the intrinsic mode functions generated from the EEMD (Wu and Huang, 2009).
The time series of state-to-state transition frequencies of the modified Markov chain model, x sstf (t), is taken as the signal.In each realization of the experiment, white noise, w(t), is added to the signal.One might interpret the added Gaussian white noise as the possible random noise that would be encountered in the measurement process or in certain restrictions applied to the calculation of edge weights in the modified Markov chain.So, for the ith realization, For each realization, we decompose the data with the added Gaussian white noise into intrinsic mode functions (IMFs).
We consider the ensemble means of the IMFs of the decompositions as the final result.Wu and Huang (2009) recommended that the ensemble size should be kept large and the amplitude of the added noise should not be small.We set the ensemble number for the number of realizations in EEMD large such that the noise series cancel each other in the final mean of the corresponding IMFs.For the two parameters, we used an ensemble size of 1000 and added noise with an amplitude of 0.2 times the standard deviation of the original data.We assume that the IMFs resulting from the EEMD represent a substantial improvement over the IMFs of the original EMD in that it utilizes the full advantage of the statistical characteristics of white noise to perturb the signal in its true solution neighbourhood, and to cancel itself after serving its purpose (Wu and Huang, 2009).EEMD results are summarized in Fig. 3a-t with intrinsic mode function followed by their state-to-state relative weight matrix derived from the basis set of the intrinsic mode functions of the time series in a fashion identical to the original time series.By summing the weights of the recurrence arcs corresponding to occurrences from state i to state j in consecutive time intervals, we calculate the weighted matrix for state-to-state transitions for each intrinsic mode function.Since the intrinsic mode functions are the mathematical basis set of the original time series, their static displays or the weighted matrices show negative values.Identical to the sum of the intrinsic mode functions yielding the original time series, the sum of the weighted matrices yields its transition frequency matrix.Similar to what Huang et al. (1998, Fig. 6 in their paper) have observed with the wind data, all of the intrinsic mode functions excluding the trend for the wind data contain both positive and negative values.We observe the same thing with the time series in that the transition probability values for the intrinsic mode functions show both positive and negative values except that the first two intrinsic mode functions have negative values larger than the lowest positive value of the trend.So, it is not surprising that the transition frequency matrices of the intrinsic mode functions  contain the positive and negative numbers.However, viewing each intrinsic mode function with the trend starting with the third IMF will obviate this difficulty in that the high-or low-frequency fluctuations ride on the trend with no negative values and the corresponding transition frequency matrices are positive.A similar observation has been made by Huang et al. (1998, Fig. 7) with their wind data.The observation made with the first two intrinsic mode functions suggests a limit on the proposed method, and it would require further investigation.The decomposition of the original time series into intrinsic mode functions and the trend is dyadic in nature, as shown in Fig. 3.This means that as we go from the first intrinsic mode function to the second and so on, the interval increases by a factor of 2 from t = 9 days to t = 18 days and so on.With an increase in the time interval from one IMF to the next, we observe the relative weights of the state-to-state transitions to vary.We also find that the state-to-state transitions within each IMF occur in packets, and the number of packets progressively decreases.The last packet of state-to-state transitions is persistent over the first eight IMFs corresponding to a time interval of 9-1152 days suggests the importance of the zone 4 earthquakes in understanding the earthquake sequencing.Although zone 4 earthquakes persist in the state-to-state transitions in the first few intrinsic mode functions, the participation of other zones in state-to-state transitions becomes significant in the higher intrinsic mode functions, IMFs 6-9.
In general, the intrinsic mode functions are characterized by (1) a certain number of a pattern of rise and fall of the arc weights and (2) by a systematic decrease in the frequency of the number of such patterns as one goes intrinsic mode function 1 to the intrinsic mode function 9. Since the rise and fall of the arc weights covers the entire catalogue of data, the periodicity that we notice could be intrinsic to earthquake processes.
The Hilbert-Huang amplitude spectrum of the time series, shown in Fig. 4, reveals at least two important features: (1) the temporal fluctuations in amplitudes occur in packets, each packet containing a set of zone to zone interactions.The oscillatory behaviour of packets contains certain periodicity within the earthquake sequence.A periodic trend at low frequencies suggests the role of zone 4 (Trenches) and zone 0 (Intraplate).A higher power at 900 and 950 time interval indicates the importance of zone 4 with earthquakes of larger magnitude prompting a cascade of aftershocks in zone 4 and main shocks in zones that are in close proximity to zone 4.
(2) The frequency-dependence of amplitude packets encapsulates the relative importance of the interaction among multiple zones over different time intervals.We interpret them to mean that certain state-to-state transitions involving zone 4 are important over a range of frequencies.

Evaluation of fractality in a state-to-state transition frequency sequence
Earthquake occurrences have been modelled to be stochastic point processes (Thurner et al., 1997;Telesca, 2005;Telesca et al., 2001Telesca et al., , 2009Telesca et al., and 2011;;Flores-Marquez and Valverde-Esparza, 2012).One representation of the point process is to examine the inter-event time intervals.The resulting interevent interval probability density function says something about the behaviour of the times between events.We do not know anything about the information contained in the relationships among these items.Since successive events do not occur in constant time intervals, another representation of a point process is given by dividing the time axis into equally spaced contiguous counting windows of duration τ , and producing a sequence of counts that fall within each time window.For example, for the kth time window, the expression for the number of counts, N k (τ ), is given by where N k (τ ) is the number of earthquakes in the kth window (Fig. 5a, b, c, d).The correlation in the process {N k (τ )} is the correlation in the underlying point process (Lowen andTeich, 1993, 1995;Thurner et al., 1997;Telesca, 2005;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011) ) have accessed such a representation of the point-processes to underscore the existence or non-existence of fractality in them.They have two calculable measures, Fano factor (FF) and Allan factor (AF), to quantify the fractality of the process (Lowen andTeich, 1993, 1995;Thurner et al., 1997;Telesca, 2005;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012).The Fano factor is a measure of correlation over different timescales (Thurner et al., 1997).It is defined as the ratio of the variance of the number of events in a specified counting time τ to the mean number of events in the counting time, as is given by where denotes the expectation value.Lowen and Teich (1995) point out that the FF of a fractal point process follows a power law with the power-law exponent, α, obeying 0 < α < 1.In other words, the FF is always greater than 1.For Poisson processes, the FF is always near unity for all counting times, and the fractal exponent is approximately equal to zero.
The Allan factor is a relation with the variability of successive counts (Allan, 1996;Barnes and Allan, 1966).It is the ratio of the variance of successive counts for a specified counting time τ divided by twice the mean number of events in the counting time.The expression of AF is given as Similar to the FF, the AF assumes values near unity for Poisson processes.Telesca et al. (2009Telesca et al. ( , 2011; henceforth, referred to as Telesca's approach) and Flores-Marquez and Valverde-Esparza (2012) have shown the power-law exponent for the AF to be 0 < α < 1.
In this paper, we examine both the results of Telesca's approach to the initial catalogue of the data used and of the new representation of the point process with a Markov chain model.For the working model, we compute the state-to-state transition frequencies as described by Nava et al. (2005) and as applied to global seismicity (Vasudevan and Cavers, 2012;Cavers and Vasudevan, 2013).Expressions similar to Eqs. ( 6) and ( 7) can be derived if we know the optimal time interval for the Markov chain model.Since we know the optimal time interval, we introduce a sequence of state-to-state transition frequencies, {N sstf,k (τ )}, with N sstf,k (τ ) referring to the weight of state-to-state transitions over the kth window for the optimal time interval, as is shown in Fig. 5f.For an easy understanding of Fig. 5f, we have included Fig. 5e.There are a few observations to be made.First, N sstf,k (τ ) is not necessarily an integer number for any kth window.Following the definition of a state, in the context of a directed graph of a Markov chain model, a state-to-state transition refers to an edge of a graph.It is the weight associated with the edge of the directed graph that plays an important role.Since we have used a modified Markov chain model which includes the influence of the event recurrences in the record-breaking sense, the above expression includes their weights as well in the computation of N sstf,k (τ ).The sequence of state-to-state transition frequencies, {N sstf,k (τ )}, yields a time series.This time series is the new expression of the point-process where the weighted edges of directed graph of the modified Markov chain represent the significance of the earthquakes between states.This new alternative representation signifies the behaviour of the state-to-state transition frequencies over a large time window.Here, seeking to find the time-correlative behaviour of the time series would be of great importance since this would give us an opportunity to see the interaction of zones considered in a collective sense.
Here, we seek to understand the correlative behaviour by looking at the two statistical measures, FF sstf and AF sstf , as defined below: The behaviour of the two measures, FF sstf and AF sstf , with respect to the optimal time interval should shed some light on the correlative behaviour of the time series but also on the selective clustering of the certain state-to-state transitions.
We consider this knowledge to be useful for forecasting purposes.
In our adaptation of the sum of edge weights for the stateto-state transition frequencies as a new representation of a point-process embedded in the modified Markov chain here, the arguments of Thurner et al. (1997), Telesca (2005) and Telesca et al. (2001Telesca et al. ( , 2009Telesca et al. ( , 2011) ) would apply.This means that the FF of the modified Markov chain sequence would follow a power law with the power-law exponent, α, satisfying 0 < α < 1.
Extending this to FF sstf and AF sstf , as is shown in Fig. 6c  and d, we find that the power law exponent calculated, corresponding to the least-squares fit of the data is greater than zero (0.27 and 0.30 respectively).They suggest not only the fractality of the modified Markov chain sequence for optimal time interval but also the deviation from the Poissonian behaviour of earthquake sequencing considered in this present study.
4 Discussion and conclusions Thurner et al. (1997) pointed out that the sequence of counts, generated by recording the number of events in successive counting time-windows of certain length, contained information about the point process depicted by the set of event times.This idea was further tested in understanding the dynamics of earthquake sequencing (Telesca et al., 2009(Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012), and in particular, the fractal behaviour of the sequence of counts.We know that this idea was initially restricted to the sequence of counts for varying windows of interval times.However, for comparison purposes, we calculated the Fano factor and the Alan factor for the initial catalogue of data using Eq. ( 6) and ( 7).We include their graphs in Fig. 6a and b.Similar to observations made by Telesca et al. (2009) with the earthquake data from the Taiwan region, we find the presence of two distinctly different regions of scaling behaviour.For small time intervals, we also observe the Poisson behaviour.Since very The Fano and Allan factor graphs, respectively, derived from the earthquake catalogue data using the approach of Telesca et al. (2001Telesca et al. ( , 2009Telesca et al. ( , 2011) ) and Telesca and Lovallo (2008) ; (c, d) the Fano and Allan factor graphs, respectively, for the time series of the state-to-state transition frequencies of the modified Markov chain model of the earthquake sequencing.
poor statistics at time-scales larger than 10 8 s would influence the Fano factor and the Allan factor, we have restricted our analysis to 3.17 years or roughly 10 8 s.
In our description of the directed graph of the Markov chain model of any earthquake sequencing, regional or global, we stress the significance of the state-to-state transition probabilities for multiple zones that span the sequence of earthquakes over an optimal time window (Cavers and Vasudevan, 2013;Vasudevan and Cavers, 2013).In other words, the edges of the directed graph carry weights.We conjecture that these weights represent a new definition of the point process.Furthermore, a consideration of the earthquake recurrences within each zone and among zones, following the concept of recurrences in the record-breaking sense (Davidsen et al., 2008), leads to an empirically determined distancedependent weights for the edges.Unlike extending the idea of the sequence of counts where every event occurrence augments the counting value by unity (Thurner et al., 1997;Telesca et al., 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012), we consider the summing of the weights for each edge such that the sum represents a "pulse" for each state-to-state transition.We analyse the resulting time series from the point of view of its Fano factor and Allan factor.
There is evidence for fractality of the multi-state modified Markov chain to represent the earthquake sequencing, as is revealed by the power-law scaling behaviour present in the Fano and Allan factors with their respective exponents of 0.27 and 0.30 (Fig. 6c, d).However, it is important to note that the exponents of the power laws in both cases have a smaller value than those observed for the initial catalogue.Cavers and Vasudevan (2013) interpreted the Markov chain of 32-states for five distinctly different zones to contain the basic combinatoric structure superimposed by the thumb-print of the undulatory structure of the recurrence weights.Since the earthquake sequencing is in general nonlinear and non-stationary, we contend that the time series representing the above Markov chain is also non-linear and non-stationary, and is conducive to an ensemble empirical mode decomposition (EEMD) procedure to understand its intrinsic mode functions (IMFs).The ensemble empirical model decomposition of the time series leads to nine intrinsic mode functions and a trend.Each one of the IMFs reveals the amplitude fluctuation of the state-to-state transitions.While there is a commonality in the relative dominance of the subduction-style earthquakes, represented by the top right corner grid of the relative weight matrices (Fig. 3), the presence or absence of certain state-to-state transitions in certain IMFs reveals the importance of integral multiples of the optimal time interval.
A simple observation of the first six or seven IMFs stresses the importance of multiple-zone approach to global seismicity problem in that the earthquake sequencing for the time period we considered has similar oscillatory behaviour of the state-to-state transition probabilities from the point of view of the amplitude scaling and the oscillating period.The growth and decay of oscillations in easily identifiable packets in each IMF following certain periodicity is an intrinsic signature of the role of multiple zones in earthquake sequencing.

Figure 1 .
Figure 1.A graph representation of earthquake sequencing with arcs (with weights w ij ) representing transitions between states.

Figure 2 .
Figure 2. (a) A time series of the state-to-state transition frequencies of the modified Markov chain model of the earthquake sequencing.The sampling time ( t) of 9 days is used.(b) The state-to-state transition frequencies of the modified Markov chain model of the earthquake sequencing.

Figure 3 .
Figure 3. Ensemble empirical mode decomposition of the time series.(a-i) Intrinsic mode functions from the first to the ninth; (j) intrinsic mode function of the trend; (k-s) state-to-state relative weight matrices for the intrinsic mode functions from the first to the ninth; (t) stateto-state relative weight matrix of the trend.Time steps and the corresponding calendar dates: 0 t -1 January 1982; 200 t -6 December 1986; 400 t -10 November 1991; 600 t -14 October 1996; 800 t -18 September 2001; 1000 t -23 August 2006; 1024 t -27 March 2007.We provide this information here to avoid any cluttering of the plots.

Figure 4 .
Figure 4. Hilbert-Huang amplitude spectrum of the intrinsic functions.

Figure 5 .
Figure 5. Representation of a point process (a-d) versus representation of a state-to-state transition (e and f).(Adapted fromThurner et al., 1997.)

Figure
Figure 6.(a, b)The Fano and Allan factor graphs, respectively, derived from the earthquake catalogue data using the approach ofTelesca et al. (2001Telesca et al. ( , 2009Telesca et al. ( , 2011) ) andTelesca and Lovallo (2008); (c, d) the Fano and Allan factor graphs, respectively, for the time series of the state-to-state transition frequencies of the modified Markov chain model of the earthquake sequencing.

Table 1 .
Tectonic zone identifier, tectonic zone and the number of earthquakes considered for M w > 5.6 and depth < 70 km from 1 January 1982 to 31 March 2008.