Statistical Significance of Earth’s Electric and Magnetic Field Variations Preceding Earthquakes in Greece and Japan Revisited

By analyzing the seismicity in a new time domain, termed natural time, we recently found that the change of the entropy under time reversal (Physica A 2018, 506, 625–634) and the relevant complexity measures (Entropy 2018, 20, 477) exhibit pronounced variations before the occurrence of the M8.2 earthquake in Mexico on 7 September 2017. Here, the statistical significance of precursory phenomena associated with other physical properties and in particular the anomalous variations observed in the Earth’s electric and magnetic fields before earthquakes in different regions of the world and in particular in Greece since 1980s and Japan during 2001–2010 are revisited (the latter, i.e., the magnetic field variations are alternatively termed ultra low frequency (ULF) seismo-magnetic phenomena). Along these lines we employ modern statistical tools like the event coincidence analysis and the receiver operating characteristics technique. We find that these precursory variations are far beyond chance and in addition their lead times fully agree with the experimental findings in Greece since the 1980s.


Introduction
Since the 1980s continuous monitoring of the variations of the electric field (in the frequency range ≤1 Hz) of the Earth has started in Greece. This led to the conclusion that there exist transient changes of the Earth's electric field, termed Seismic Electric Signals (SES), preceding earthquakes (EQs) [1,2]. These changes have been classified into single SES (meaning that one transient change appears before an EQ) or SES activity which corresponds to the case when several transient changes are recorded within a short time [3]. The lead time ∆t in the former case is ∆t ≤ 11 days (d) while in the latter case ∆t is appreciably longer around a few months or so with maximum ≈ 5 1 2 months [4] and minimum ≈ a few weeks [3,5] (e.g., 22 d see [3]). SES are believed [5] to be generated when the gradually increasing stress in the EQ preparation area reaches a critical value σ cr and then the existing electric dipoles that have been formed due to point defects (e.g., [6]) exhibit a cooperative orientation, thus leading to the emission of a transient electric signal. This SES generation model is considered to be unique [7] among other models in that SES would be generated spontaneously during the gradual increase of stress without requiring any sudden change of stress such as micro fracturing. This model is confirmed in the following sense: Upon analysing the time series of the pulses comprising an SES activity in a new time domain, termed natural time χ, it has been shown [8,9] that these pulses exhibit infinitely ranged temporal correlations which is characteristic of criticality. For a time series comprising N events, we define an index for the occurrence of the k-th event by χ k = k/N, which we term natural EQs with M J MA ≥ 7.6 that occurred in Japan during this period. It was later shown [35] by means of Monte-Carlo calculation as well as by the receiver operating characteristics (ROC) technique that the probability to obtain the above result by chance is of the order of 10 −5 . In addition, Varotsos et al. [32] showed that the two phenomena (initiation of SES activity and minimum of the order parameter fluctuations) are also linked closely in space. This opened the window for a reliable estimation of the epicentral area of an impending EQ, as subsequently confirmed by the estimation [36] of the epicentral area for all shallow major mainshocks of magnitude 7.6 or larger that occurred in Japan during the aforementioned period 1984-2011.
It is the scope of this paper to investigate the statistical significance of the aforementioned precursory phenomena and in particular the anomalous variations observed in the Earth's electric and magnetic field before EQs. Along these lines we will study two independent datasets, one from Greece and another one from Japan. Specifically, the dataset from Greece comprises all the SES reported during the almost 2 1 2 year period from 1 April 1987 until 30 November 1989 that possibly preceded EQs within the area N 41. 5 36.0 E 26.0 19.0 (see Figure 1). The selection of this early dataset is made for two reasons: First, the statistical significance of the data set has been previously studied in other works (e.g., [37]) showing that SES predictions are far beyond chance but with conventional techniques. Second, during this early period, predictions based on SES were issued when the expected EQ magnitude was 5 or larger, while later, i.e., after mid-1990s the European Advisory Committee for EQ prediction of the Council of Europe recommended (see p. 102 of [38]) that predictions should be issued only for large magnitudes, i.e., with local magnitudes M L ≥ 5.5 and hence Ms(ATH) ≥ 6.0, thus the number of predictions based on SES became considerably smaller.  [3] in Greece during the period from 1 April 1987 to 30 November 1989 are depicted by the red triangles while the blue stars correspond to the epicentres of the earthquakes (EQs) with local earthquake magnitude reported by GI-NOA (ML(ATH)) ≥ 4.5. The dates of occurrence and their magnitudes are shown in Figure 2.
As for the dataset from Japan, it comprises the geomagnetic data and local EQs registered in Kakioka (KAK) station in Japan during the 10 year period 2001-2010, which has been recently studied by Han et al. [29]. They utilized Molchan's error diagram [39] to evaluate whether geomagnetic changes contain precursory information beyond chance in a strikingly similar fashion with the procedure Entropy 2018, 20, 561 4 of 17 followed in [40,41] which also employed Molchan's diagram to show that the predictions issued on the basis of SES during the period 1987-1995 were appreciably better than random score (see Figure 2 of [40]). Han et al. [29] concluded that their magnetic field precursory variations perform better than chance when the lead time ∆ is around one week and the alarm window L is less than four d or ∆ is 13-14 d and L is less than one week. Interestingly, the former value almost coincides with the lead time reported for the single SES in Greece and the quantity ∆ + L of the latter agrees only with the minimum (≈22 d) of the lead time of SES activities observed both in Greece and Japan (cf. in these cases ∆t may be appreciably longer, i.e., around a few months or so). We report for example two such cases from Japan (cf. several similar cases in Greece are reported in [5]): First, the magnetic (and electric) field change before the volcanic-seismic swarm activity in 2000 in the Izu Island region started more than two months before [14,15]. (This lead time was attributed by Han et al. [29] to the volcanic-seismic nature in the Izu Island activity, but they did not consider that similarly long lead times had been recorded repeatedly in Greece for tectonic EQs since the early 1990s [5,42].) The same holds for the magnetic field anomaly observed before the M9 Tohoku EQ since it appeared (e.g., [43,44]) during the week around 5 January 2011 and the EQ occurred in 11 March 2011 (cf. on 5 January 2011 natural time analysis has revealed [34] that the minimum of the order parameter of seismicity has been observed which, as mentioned, may allow an estimation of both the time of occurrence [45] and the epicentral location [46] of the forthcoming mainshock). Thus, in short, Han et al. [29] claim that the maximum lead time is ∆ + L = 22 d which however does not agree either with longer lead times reported in Greece or with the aforementioned established examples in Japan for the cases of Izu and Tohoku. The present work sheds light on this disagreement.
In the present study, we solely employ modern statistical tools like the event coincidence analysis (ECA) [47] and the receiver operating characteristics (ROC) technique [48]. We clarify in advance that here we restrict ourselves to the probabilities to obtain by chance the results referring solely to the occurrence times of the events while the relevant calculation for achieving all the parameters of an impending EQ (time, epicenter, and magnitude) should also consider the probabilities to obtain the epicentral area and the magnitude of the impending EQs, e.g., [49,50].

Event Coincidence Analysis (ECA)
ECA has been designed [47] to quantify the strength, directionality and time lag of the statistical relations between event series. According to Donges et al. [47]: the method was introduced in a less general setting to study possible statistical interrelationships between nonlinear regime shifts in African peleoclimate during the past 5 million years and events in hominin evolution such as the appearance and disappearance of species [51]. ECA considers [47] , and ECA is based on counting (possibly lagged) coincidences between events of different types. In this context [52], B type events are considered as possibly influencing the timings of A type events, and not vice versa. The assumption to be quantified and tested by ECA is that the events in B precede the events in A. This is made by introducing an instantaneous coincidence if two events with timings t A i and t B j (t A i ≥ t B j ) are closer in time than a coincidence interval ∆T, that means and generalizing it to a lagged coincidence if Entropy 2018, 20, 561 5 of 17 holds, where τ ≥ 0 is the time lag parameter. In order to quantify the strength of the interrelations between the two event series, two variants of the coincidence rate addressing B type events as either precursors or trigger have been introduced [47]: The precursor coincidence rate where Θ(x) is the Heaviside unit step function (equal to 0 for x ≤ 0 and 1 for x > 0) and 1 I (x) is the indicator function for the interval I (equal to 1 for x ∈ I and 0 otherwise), measures the fraction of A type events that are preceded at least by one B type event (i.e., multiple B type events within ∆T are counted once) and the trigger coincidence rate measures the fraction of B type events that are followed by at least one A type event (i.e., multiple A type events within ∆T are counted once). The distinction between precursor and trigger coincidence rates allows the introduction of a notion of directionality while the parameter τ explicitly takes into account lagged interrelations between B type and A type events. Based on these two coincidence rates and appropriate assumptions (for example for the inter-event time distributions [47]) various statistical tests can be made to examine whether B type events are precursors to A type events for a risk enhancement test [53] or whether B type events are triggers for A type events (trigger test [53]). Typical examples are the climate-related disasters as risk enhancement factor for armed-conflicts in ethnically fractionalized countries [53] or the role of flood events as triggers of epidemic outbreaks [47].
Here, we employed the function CC.eca.es of the CoinCalc package [52] for R [54] that implements ECA and selected the (default) Poisson test. The reason behind this selection, although it is known that EQs appear in sequences due to aftershocks, was that the magnitude range of the EQs (that define A type events) in most of the studies here barely exceeds one magnitude unit while according to Båth law [4,[55][56][57][58] the difference in magnitude between the mainshock and its largest aftershock is approximately 1.2 magnitude units. In other words, in the data presented in Section 3 below, it is improbable that aftershocks are used for the determination of A type events (cf. only in Section 3.2 and when we consider all the 50 EQs of  Table 2 of [59]).

ROC Analysis
ROC analysis is a modern technique [48] that depicts the trade off between hit rates and false alarm rates in a conceptually simple way. It has been recently applied [35,[60][61][62][63][64][65][66][67] for estimating the predictability of various complex systems and as such it might be useful and complementary to ECA, which was mentioned in the previous subsection.
Suppose that we have in total N events out of which P are strong and important to be predicted and Q = N − P of which are not as strong (for example in the case of EQs one may consider as strong an EQ of magnitude M exceeding some threshold M thres , i.e., M ≥ M thres , which might be also called "strong event", while smaller EQs may be considered as "non-events"). Suppose also that before the occurrence of each event either strong or not, a quantity labeled is estimated based on some prediction algorithm. If ≥ t , where t is some threshold value, we assume that the forthcoming event will be of the P class and if < t a Q class event is expected. Upon the occurrence of the forthcoming event, four situations may appear: (a) If ≥ t and the event is P class we have a successful prediction called True Positive (TP); (b) if < t and the event is Q class we have a successful prediction called True ROC analysis is based on a graph in which we plot the True Positive rate (TPr), or hit rate (H), which is the number |TP| of TPs divided by P: versus the False Positive rate (FPr), or false alarm rate (F), which is the number |FP| of FPs divided by Q: A ROC curve is obtained as we vary t and plot TPr as a function of FPr. Thus, each prediction scheme (e.g., selection of t ) may be considered as a point in the ROC curve. If the predictor is decided at random, the ROC curve will scatter close to the diagonal (which will result if both P and Q tend to infinity) but for real problems it fluctuates as P and Q vary. The statistical significance of a ROC curve depends on the area above the horizontal axis (TPr = 0) and under the curve [68]. By assuming an appropriate set of ellipses, called k-ellipses, we may estimate [66] the p-value to obtain by chance a point on the ROC. Here, we employed the computer code VISROC [66] that allows the visualization of this p-value to obtain by chance a point on the ROC graph and hence this enables us to perform statistical tests on whether a prediction scheme outperforms chance.

Greece
In Table 1 we give all the SES reported [3,5,69] during the period from 1 April 1987 until 30 November 1989. The dates of all these SES are plotted in Figure 2 as vertical red lines along with the local magnitude ML(ATH) (≥4.5) of the EQs that occurred within the area N 41. 5 36.0 E 26.0 19.0 reported by GI-NOA (for their epicentres see Figure 1). By applying ECA between the event series corresponding to the occurrence dates of these EQs (event time series A, see Section 2.1) and the event series of the SES dates of Table 1 (event time series B, see Section 2.1), we find (see Figure 3) that for τ = 3 d and ∆T = 4 d the precursor p-value to obtain by chance the real interconnection between the event time series A and B is 2.7% (while the corresponding precursor coincidence rate, which equals the hit rate, is 26.5%). In other words, the SES are statistically significant precursors to Ms(ATH) ≥ 5.0 EQs. If we increase ∆T to 6 d, we obtain a p-value 4.5% with r p (∆T = 6 d, τ = 3 d) = 32.4% also pointing to statistical significance. Upon focusing on the cases when r p exceeds r t (see red squares in Figure 3c), we also found that SES are statistically significant precursors to the EQs depicted in Figure 2 in the following three cases: First, the p-value results in 4.9% when using τ = 18 d and ∆T = 6 d with r p (∆T = 6 d, τ = 18 d) = 32.3%; second, the p-value results in 3.4% when using τ = 43 d and ∆T = 4 d with r p (∆T = 4 d, τ = 43 d) = 26.5%; third, the p-value results in 3.7% when using τ = 58 d and ∆T = 4 d with r p (∆T = 4 d, τ = 58 d) = 26.5%. Additional SES activities verifying these three cases can be found in [5,42,70,71]. By summarizing, we found that SES are statistical significant precursors to EQs for the following four distinct time periods [τ, τ + ∆T]: 3 to 9 d, 18 to 24 d, 43 to 47 d and 58 to 62 d which are strikingly reminiscent of the finding [3][4][5]42] that single SES have a ∆t ≤ 11 d while for SES activities the minimum lead time ∆t is around a few weeks and the maximum 5 1 2 months, as mentioned in the introduction. Since in a finite sample, as studied here, there is always a chance of false positive significance in a certain fraction of individual tests, we also repeated ECA 30 times with the SES dates randomly chosen. Out of these 30 cases studied, only three of them led to smaller than six and larger than zero statistically significant pairs (∆T, τ) for which r p exceeds r t . Moreover, in all these three cases the obtained τ values included consecutive τ's (e.g., in one such case the τ-values 10, 11, and 49 were obtained). In the real case studied above, however, there were no consecutive τ values which points to the existence of a physical mechanism rather than a chancy correlation. In other words, the probability of obtaining by chance the four distinct periods [τ, τ + ∆T]: 3 to 9 d, 18 to 24 d, 43 to 47 d and 58 to 62 d is below 1/30 ≈ 3.3%.
Based on these results, we also applied the ROC technique by dividing the study period into consecutive weeks (i.e., a time interval corresponding to the larger two ∆T values of the previous paragraph). Based on whether a certain week included a ML(ATH) ≥ 4.5 EQ (as depicted in Figure 2) or not, we decided if this week had to be predicted (strong event, see Section 2.2) or not. This resulted in 28 weeks including such EQs out of the 140 weeks in total. The SES dates were also segmented (coarse grained) into weeks for which an alarm was issued or not. This procedure has led to 27 alarm weeks. The operation points for example of the two cases: (∆T = 6 d, τ = 3 d) and (∆T = 6 d, τ = 18 d) out of the four cases identified by ECA and discussed in the previous paragraph are shown by the circle and the square, respectively, in the ROC depicted in Figure 4. The p-values to obtain these results by chance are 2.0% and 0.5%, respectively, which are statistically significant.

Japan
Peng et al. [59] defined the energy enhancement parameter P, as the ratio of the observed energy (Z) at Kakioka (KAK) station in the vertical component of the magnetic field in the ULF region, i.e., at around [28] 0.01 Hz, over its modelled value (Z * ) estimated on the basis of the measurements made at a remote (≈1000 km away) station [59]. This parameter, given in Figure 2a of Han et al. [29], is plotted for the 10-year period 2001-2010 in Figure 5 together with the EQ magnitude M J MA for 50 sizeable EQs (see Table 2 of [59]) that occurred within 100 km from KAK, which according to [29,59] could have been preceded by energy enhancements in the ULF variation of the vertical component of the geomagnetic field. Here, following [29] we will assume that when P exceeds a threshold value P * an alarm for an impending EQ is issued. For example, in Figure 5 the red vertical lines correspond to the dates for which P exceeds the threshold value P * = 3, i.e., when P > 3. These dates as a first approximation have been selected as a coarse-grained average over three-day periods from Figure 2a of Han et al. [29] and lead to the 34 alarm dates given in Table 2.
We start with ECA for the 10-year period 2001-2010 (3652 d) and consider the event series of the 50 EQs shown in Figure 5 as event time series A and the event series comprising the 34 alarm dates shown in Table 2 as event series B. We investigated again the cases for which the precursor coincidence rate exceeds the trigger coincidence rate (see the red square in Figure 6), since we want to perform a risk enhancement test. Such an analysis resulted in one statistically significant solution (∆T = 1 d, τ = 32 d) with p-value 1.2%. The corresponding hit rate is r p (∆T = 1 d, τ = 32 d) = 6%. The colours indicate (a) the precursor p-value, (b) the precursor coincidence rate, and (c) the ratio of the precursor over the trigger coincidence rate. These values were obtained by applying CoinCalc [52] as discussed in subsection 2.1. The red squares indicate the statistically significant cases of (∆T, τ) that apart from having a precursory p-value smaller than 5%, they also exhibit a precursory coincidence rate exceeding the trigger coincidence rate.
Based on these results, we also applied the ROC technique by dividing the study period into (b) the precursor coincidence rate; and (c) the ratio of the precursor over the trigger coincidence rate. These values were obtained by applying CoinCalc [52] as discussed in Section 2.1. The red squares indicate the statistically significant cases of (∆T, τ) that apart from having a precursory p-value smaller than 5%, they also exhibit a precursory coincidence rate exceeding the trigger coincidence rate.  Tables 1, 2 and  4 of [5], see also [69]). The sites of the measuring stations are shown in Figure 1.     Tables 1, 2 and  4 of [5], see also [69]). The sites of the measuring stations are shown in Figure 1.  [59]. This parameter, given in Fig.2a of Han et al. [29], is plotted for 250 the ten year period 2001-2010 in Figure 5 together with the EQ magnitude M J MA for 50 sizeable EQs

251
(see Table 2 of [59]) that occurred within 100km from KAK, which according to [29,59]  field. Here, following [29] we will assume that when P exceeds a threshold value P * an alarm for an 254 impending EQ is issued. For example, in Figure 5 the red vertical lines correspond to the dates for 255 which P exceeds the threshold value P * = 3, i.e., when P> 3. These dates as a first approximation have 256 been selected as a coarse-grained average over three days periods from Fig.2a of Han et al. [29] and 257 lead to the 34 alarms dates given in Table 2.  Table 2 of [59] that occurred within 100km from KAK. The red vertical lines correspond to the dates (see Table 2) when P exceeds the threshold value P * = 3, i.e., when P> 3.
We start with ECA for the ten year period 2001-2010 (3652d) and consider the event series of 259 the 50 EQs shown in Figure 5 as event time series A and the event series comprising the 34 alarm 260 dates shown in Table 2 as event series B. We investigated again the cases for which the precursor 261 coincidence rate exceeds the trigger coincidence rate (see the red square in Figure 6), since we want 262 to perform a risk enhancement test. Such an analysis resulted in one statistically significant solution 263 (∆T = 1d, τ = 32d) with p-value 1.2%. The corresponding hit rate is r p (∆T = 1d, τ = 32d) = 6%.
264 Figure 5. The energy enhancement parameter P (green) as reported in Figure 2a of Han et al. [29] together the EQ magnitude M J MA (blue candlesticks ending at circles, right scale) reported by Japan Meteorological Agency (JMA) for the 50 EQs of Table 2 of [59] that occurred within 100 km from Kakioka (KAK). The red vertical lines correspond to the dates (see Table 2) when P exceeds the threshold value P * = 3, i.e., when P > 3. Hereafter, we focus on the large magnitude EQs (see Figure 7) searching for statistically significant cases with the highest hit rate and the smallest p-value. We first investigate the EQs with M J MA ≥ 5.5 (six EQs in total). Such a statistically significant solution was obtained for (∆T = 13 d, τ = 90 d) with a p-value 2.4% and r p (∆T = 13 d, τ = 90 d) = 50%. Secondly, we investigate the EQs with M J MA ≥ 5.8 (four EQs in total). Such a statistically significant solution was again obtained for (∆T = 13 d, τ = 90 d) with a p-value 0.6% and r p (∆T = 13 d, τ = 90 d) = 75%.
We now turn to ROC analysis. Since the results of the previous paragraph point to the importance of ∆T = 13 d and hence a two-week period, we divided the study period of 3652 d into 261 consecutive 14 d periods. Based on whether a certain two-week period included a M J MA ≥ 5.5 EQ or not, we decided if this period had to be predicted (strong event, see Section 2.2) or not. This resulted in six periods including such EQs out of the 261 in total. The alarm dates were also segmented (coarse grained) into two-week periods for which an alarm was issued or not. This procedure led to 29 alarms, which when considering that three alarms were followed by at least one EQ, we obtain an FPr F = 10.2% while the hit rate is H = 50% as mentioned. The operation point of the pair (∆T = 13 d, τ = 90 d) discussed in the previous paragraph is shown by the circle in the ROC depicted in Figure 8a. The p-value to obtain this result by chance is 0.1%. By applying the same procedure to predict M J MA ≥ 5.8 EQs, we obtain the operation point indicated by the circle in the ROC analysis shown in Figure 8b. The p-value to obtain this result by chance is again found to be 0.1%. p−value Figure 6. ECA for the ten year period 2001-2010 in Japan when studying the magnetic field data shown in Figure 5. The colours indicate (a) the precursor p-value and (b) the precursor coincidence rate which were obtained by applying CoinCalc [52]. The red square indicates the single statistically significant combination (∆T = 1d, τ = 32d) that apart from having a precursory p-value smaller than 5%, it also exhibits a precursory coincidence rate exceeding the trigger coincidence rate. Figure 6. ECA for the 10-year period 2001-2010 in Japan when studying the magnetic field data shown in Figure 5. The colours indicate (a) the precursor p-value and (b) the precursor coincidence rate which were obtained by applying CoinCalc [52]. The red square indicates the single statistically significant combination (∆T = 1 d, τ = 32 d) that apart from having a precursory p-value smaller than 5%, it also exhibits a precursory coincidence rate exceeding the trigger coincidence rate.   Figure 5 and focusing on the prediction of large magnitude EQs. The colours indicate in (a) and (c) the precursor p-value and in (b) and (d) the precursor coincidence rate which were obtained by applying CoinCalc [52]. Panels (a),(b) correspond to M J MA ≥ 5.5 while (c),(d) to M J MA ≥ 5.8. The red squares indicate the statistically significant cases of (∆T, τ) that apart from having a precursory p-value smaller than 5%, they also exhibit a precursory coincidence rate exceeding the trigger coincidence rate. CoinCalc [52]. Panels (a,b) correspond to M J MA ≥ 5.5 while (c,d) to M J MA ≥ 5.8. The red squares indicate the statistically significant cases of (∆T, τ) that apart from having a precursory p-value smaller than 5%, they also exhibit a precursory coincidence rate exceeding the trigger coincidence rate.

283
We first comment on the results obtained in subsection 3. be interesting to study the existing geoelectric records and vary the threshold for the identification of 300 SES activities so that the operation points in the ROC diagram of Fig.4 can be extended to ROC curves.

301
This way an optimal operation can be selected. Such a study, however, falls beyond the scope of the 302 present paper.

303
As for the results of ECA in Japan, it is striking that they do not reveal any of the two time periods,

Discussion
We first comment on the results obtained in Section 3.1 in which by means of ECA, we found that SES recorded since 1980s in Greece are statistically significant precursors to EQs for four distinct time periods 3 to 9 d, 18 [4,24,38] by an SES activity on 30 April 1995; for the period 58 to 62 d, the M6.6 Leonidion EQ on 6 January 2008 that was preceded [4,70] by an SES activity on 7 November 2007. We recall that during the almost 40-year period from the 1980s to 2018 in Greece, several large magnitude EQs 5.5 < M < 7 occurred and their preceding SES activities had lead times up to 5 1 2 months [3][4][5]42,71]. As mentioned in the introduction, since the mid-1990s SES predictions are issued when the expected magnitude is Ms(ATH) ≥ 6.0. A study concerning such predictions for the period 2001 to 2011 can be found in Table 7.1 of [4]. The results that are summarized in Section 7.3 of [4] show that 15 out of the 16 EQs with Ms(ATH) ≥ 6.0 where preceded by SES. This points to the direction that if we consider higher magnitude thresholds, better results are expected. In the same context, it would be interesting to study the existing geoelectric records and vary the threshold for the identification of SES activities so that the operation points in the ROC diagram of Figure 4 can be extended to ROC curves. This way an optimal operation can be selected. Such a study, however, falls beyond the scope of the present paper.
As for the results of ECA in Japan, it is striking that they do not reveal any of the two time periods, i.e., (∆ = 1 week, L < 4 d) and (∆ = 13-14 d, L < 1 week), deduced by Han et al. [29] from the analysis of the Japanese magnetic field data during 2001-2010. To the contrary, the application of ECA to Japanese data resulted in a lead time of around one month, i.e., 32 to 33 d, and in particular for the six large magnitude EQs with M J MA > 5.5 during 2001-2010 it led to even longer lead times around 3 1 2 months. Quite interestingly, these lead times identified by ECA agree with experimental results obtained in Greece, but have not been identified by Han et al. [29] (for example the M6.5 Andravida EQ in Greece on 8 June 2008, was preceded [4,70] by an SES activity that initiated on 29 February 2008). These new findings could be attributed to the superiority of the modern statistical tools employed in the present study.
We now comment on the point made by Han et al. [29] that a lead time of around one week can be established as an optimal strategy for earthquake prediction. As we have shown in Section 3.2, such a lead time did not result from the 10-year (2001-2010) magnetic field data analysed by Han et al. [29]. A procedure to achieve a time window of around one week for the occurrence of an impending EQ is to apply the following procedure based on natural time analysis: Just after the initiation of the SES activity, we analyse in natural time the EQs that occur in the candidate epicentral area (estimated either on the basis of the SES properties [70] or with the procedure developed in [36,46]) and compute the κ 1 value upon the occurrence of each EQ. Once this value approaches 0.070, e.g., see [4,9,70,72], the mainshock occurs within one week or so. Such an application for the M9 Tohoku EQ that occurred in 11 March 2011 in Japan has already appeared in [46]. Concerning this mega-earthquake, an interesting application of natural time analysis to ULF magnetic field variations can be found in [73].
Finally, we note that the two example cases presented here correspond to the most seismic active regions of the world (Japan) and of Europe (Greece). As such, they constitute representative examples which may find useful application to other seismic prone regions. Of course, the establishment and operation of geoelectric and/or geomagnetic measurement networks together with similar statistical studies will provide the ultimate test.

Summary and Conclusions
Here, we employed two modern statistical tools, i.e., ECA and ROC, to revisit the statistical significance of the anomalous changes of the Earth's electric and magnetic field as precursory to EQs. Two independent datasets (first, electric field data from Greece since the 1980s and second, magnetic field data from Japan during 2001-2010 which is alternatively termed ULF seismo-magnetic phenomena) were studied and we found that they exhibit precursory information far beyond chance. In addition, our results show that in both datasets the corresponding lead times are comparable to those determined experimentally in Greece since the 1980s. Quite interestingly, upon restricting ourselves to magnetic field changes before major EQs in Japan, we find that ECA succeeds in identifying their correct lead times, which are around 3 1 2 months, in agreement with comparable experimental values reported in Greece. On the other hand, previous conventional calculations for the analysis of Japanese data by other studies, e.g., [29], obtained appreciably shorter lead times (around one week and two weeks).
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: