Retrospective Evaluation of the Five-Year and Ten-Year CSEP-Italy Earthquake Forecasts

The CSI catalog spans the time period from 1981 until 2002 and reports local magnitudes, in agreement with the BSI magnitudes that will be used during the prospective evaluation of forecasts. Schorlemmer et al. [2010a] found a clear change in earthquake numbers per year in 1984 due to the numerous network changes in the early 1980s and therefore recommend using the CSI data from 1 July 1984 onwards. For the retrospective evaluation, we selected earthquakes with local magnitudes $M_{L}\geq 4.95$ from 1985 until the end of 2002. To mimic the durations of the prospective experiments, we selected three non-overlapping five-year periods (1998-2002, 1993-1997, 1988-1992). To test the robustness of the results, we also used the entire 18-year span of reliable data from 1985 until 2002. We selected shocks as test data if they occurred within the CSEP-Italy testing region [see Schorlemmer et al., 2010a].

The CPTI catalog covers the period from 1901 until 2006 and is based on both instrumental and historical observations (for details, see Schorlemmer et al. [2010a]). The catalog lists moment magnitudes that were estimated either from macroseismic data or calculated using a linear regression relationship between surface wave, body wave or local magnitudes. Because the prospective experiment will use local magnitudes, we converted the moment magnitudes to local magnitudes using the same regression equation that was used to convert the original local magnitudes to moment magnitudes for the creation of the CPTI catalog [Gruppo di lavoro MPS, 2004; Schorlemmer et al., 2010a]:

Schorlemmer et al. [2010a] estimated a conservative completeness magnitude of $M_{L}=4.5$, so that one could justify using the entire period from 1901 until 2006 for the retrospective evaluation. However, we focused mainly on the data since the 1950s because it seems to be of higher quality [Schorlemmer et al., 2010a]. We divided the period into non-overlapping ten-year periods to mimic the duration of the prospective experiment, but we also evaluated the forecasts on a 57-year time span from 1950 until 2006 and on the 106-year period from 1901 until 2006. As for the CSI catalog, we only selected shocks within the testing region. Some quakes, mostly during the early part of the CPTI catalog, were not assigned depths. We included these earthquakes as observations within the testing region because it is very unlikely that they were deeper than 30 km (see also Schorlemmer et al. [2010a]).

In this section, we consider the distribution of the number of observed events in the five- and ten-year periods relevant for the Italian forecasts. Analysis of this empirical distribution can test the assumption (made by all time-independent forecasts) that the Poisson distribution approximates well the observed variation in the number of events in each cell and in the entire testing region. (CSEP-Italy participants decided to forecast all earthquakes and not only so-called mainshocks – see section 4.4).

The Poisson distribution is defined by its discrete probability mass function:

where $n=0,1,2,...,$ and $\lambda$ is the rate parameter, the only parameter needed to define the distribution. The expected value and the variance $\sigma^{2}_{POI}$ of the Poisson distribution are both equal to $\lambda$.

Because the span of time over which the CSI catalog is reliable is so short, we used the CPTI catalog for the seismicity rate analysis. The sample variance of the distribution of the number of observed earthquakes in the CPTI catalog over non-overlapping five-year periods from 1907 until 2006 (inclusive) is $\sigma^{2}_{5yr}=23.73$, while the sample mean is $\mu_{5yr}=8.55$. For non-overlapping ten-year periods of the CPTI catalog, the sample variance is $\sigma^{2}_{10yr}=64.54$, while the sample mean is $\mu_{10yr}=17.10$. Because the sample variance is so much larger than the sample mean on the five- and ten-year timescales, it is clear that the seismicity rate varies more widely than expected by a Poisson distribution.

Figure 1 shows the number of observed earthquakes in each of the twenty non-overlapping five-year intervals, along with the empirical cumulative distribution function. The Poisson distribution with $\lambda=\mu_{5yr}=8.55$ and its $95\%$ confidence bounds are also shown. One should expect that one in twenty data points falls outside the $95\%$ confidence interval, but we observe four, and one of these lies outside the $99.99\%$ quantile.

We compared the goodness of fit of the Poisson distribution with that of a negative binomial distribution (NBD), motivated by studies that suggest its use based on empirical and theoretical considerations [Vere-Jones, 1970; Kagan, 1973; Jackson and Kagan, 1999; Kagan, 2010; Schorlemmer et al., 2010b; Werner et al., 2010a].The discrete negative binomial probability mass function is

where $n=0,1,2,...,$ $\Gamma$ is the gamma function, $0\leq\nu\leq 1$, and $\tau>0$. The average of the NBD is given by

while the variance is given by

implying that $\sigma^{2}_{NBD}\geq\sigma^{2}_{POI}$. Kagan [2010] discusses different parameterizations of the NBD. For simplicity, we used the above parameterization and maximum likelihood parameter value estimation. We found $\tau_{5yr}=6.49$ and $\nu_{5yr}=0.43$, with $95\%$ confidence bounds given by $[-0.39,13.37]$ and $[0.17,0.70]$, respectively. The large uncertainties reflect the small sample size of twenty. For the ten-year intervals, we estimated $\tau_{10yr}=9.24$ and $\nu_{10yr}=0.35$, with $95\%$ confidence bounds given by $[-2.74,21.22]$ and $[0.05,0.65]$, respectively. Figure 1 shows the $95\%$ confidence bounds of the fitted NBD in the number of observed events (left panel), and the NBD cumulative distribution function (right panel).

Because the NBD is characterized by two parameters while the Poisson has only one, we used the Akaike Information Criterion (AIC) [Akaike, 1974] to compare the fits:

where $L$ is the likelihood of the data given the fitted distribution and $p$ is the number of free parameters. For the five-year and ten-year intervals, the NBD fit the data better than the Poisson distribution, despite the penalty for the extra parameter: for the five-year intervals, $AIC_{NBD}=117.32$ and $AIC_{POI}=126.20$, while for the ten-year intervals, $AIC_{NBD}=70.05$ and $AIC_{POI}=77.56$. To test the robustness of the better fit of the NBD over the Poisson distribution, we also checked the distribution of the number of events in one-year, two-year and three-year intervals of both catalogs. In all cases, the NBD fit better than the Poisson distribution, despite the penalty for an extra parameter.

Several previous studies showed that the distribution of the number of earthquakes in any finite time period is not well approximated by a Poisson distribution and is better fit by an NBD [Kagan, 1973; Jackson and Kagan, 1999; Schorlemmer et al., 2010b; Werner et al., 2010a] or a heavy-tailed distribution [Saichev and Sornette, 2006]. The implications for the CSEP-Italy experiment, and indeed for all CSEP experiments to date, are important.

The only time-independent point process is the Poisson process [Daley and Vere-Jones, 2003]. Therefore, a non-Poissonian distribution of the number of earthquakes in a finite time-period implies that, if a point process can model earthquakes well, this process must be time-dependent (although there might be other, non-point-process classes of models that are time-independent and generate non-Poissonian distributions). Therefore, the Poisson point process representation is inadequate, even on five- or ten-year timescales for large $m\geq 4.95$ earthquakes in Italy, because the rate variability of time-independent Poisson forecasts is too small, and they will fail more often than expected. As a result, the agreement of CSEP-Italy participants to use a Poisson distribution should be viewed as problematic for time-independent models because no Poisson distribution that their model could produce will ever pass the tests with the expected $95\%$ confidence rate. On the other hand, time-dependent models with varying rate might generate an NBD over a longer time period (the empirical NBD can even be used as a constraint on the model).

Solutions to this problem have been discussed previously. The traditional approach has been to filter the data via declustering (deletion) of so-called aftershocks (as done, for instance, in the RELM mainshock experiment [Field, 2007; Schorlemmer et al., 2007]). However, the term “aftershock” is model-dependent and can only be applied retrospectively. A more objective approach is to forecast all earthquakes, allowing for time-dependence and non-Poissonian variability. In theory, each model could predict its own distribution for each space-time-magnitude bin [Werner and Sornette, 2008], and future predictability experiments should consider allowing modelers to provide such a comprehensive forecast (see also section 7).

A third, ad-hoc solution [see Werner et al., 2010a] is more practical for time-independent models in the current context. Based on an empirical estimate of the observed variability of past earthquake numbers, one can reinterpret the original Poisson forecasts of time-independent models to create forecasts that are characterized by an NBD. One can perform all tests (defined below in section 5) using the original Poisson forecasts, and repeat the tests with so-called NBD forecasts.

We created NBD forecasts for the total number of observed events by using each forecast’s mean and an imposed variance identical for all models, which we estimated either directly from the CPTI catalog or from extrapolation assuming that the observed number of events are uncorrelated. Appendix A describes the process in detail. Because the resulting NBD forecasts are tested on the same data that were used to estimate the variance, one should expect that the NBD forecasts perform well. The broader NBD results in less specificity, but also fewer unforeseen observations. We will re-examine this ad-hoc solution in the discussion in section 7.

In Figure 2, we show the results of the N, L, S, and M-tests applied to the time-independent forecasts for the most recent five-year target period from 1998-2003 of the CSI catalog. We discuss each of the test results below. As a summary of all the results we present here and below, Tables 2 and 3 list all the tests that the forecasts fail for each of the considered target periods of the CSI and CPTI catalog, respectively.

In Figure 4, we summarize the results of the N, conditional L, S and M-tests for the time-independent models and five non-overlapping ten-year target periods of the CPTI catalog. These results mimic the prospective ten-year experiment and help gauge the variability of the results. The online material that accompanies this article (available at http://mercalli.ethz.ch/~mwerner/CSEP_ITALY/ESUPP/) provides additional figures of the forecasts, maps of their likelihood ratios against a uniform forecast, and concentration diagrams [Rong and Jackson, 2002; Kagan, 2009] for the entire CPTI data set from 1901 until 2006. Because the figures are based on the longest target period, which we consider explicitly in section 6.3, they include all earthquakes observed during the ten-year target periods and provide an informative visual presentation of the results.

The long-term forecasts submitted for the CSEP-Italy experiment were calculated for five-year and ten-year periods. Because the forecasts are time-independent and characterized by Poisson uncertainty, one can test suitably scaled versions of the forecasts on longer time periods: 18 years (the duration of the reliable part of the entire CSI catalog, from 1985 through 2002), 57 years (the duration of the most reliable part of the CPTI catalog, from 1950 through 2006), and 106 years (the entire CPTI catalog). In this section, we present the results of testing these scaled forecasts. The online material presents further figures of the forecasts, likelihood ratios and concentration diagrams based on the 106-year target period.

The test results of the 18-year period from 1985 to 2002 of the CSI catalog are shown in Figure 5. Twenty-three earthquakes occurred during this period. The N-test results reveal the same features already observed previously: a group of models overpredicts the number of earthquakes (Akinci-et-al.Hazgridx, Chan-et-al.HzaTI, Meletti-et-al.MPS04, Schorlemmer-Wiemer.ALM, Zechar-Jordan.CPTI, Zechar-Jordan.Hybrid). While the confidence bounds of the negative binomial distribution remain substantially wider than the bounds based on the Poisson distribution, there are only two forecasts for which the test results are ambiguous (Akinci-et-al.Hazgridx and Nanjo-et-al.RI). The ALM forecasts and the Meletti-et-al.MPS04 forecast fail the conditional L-test and the S-test, with Schorlemmer-Wiemer.ALM scoring less than a uniform spatial forecast. The failures are due to the earthquakes we discussed previously that occur either in background bins or in locations with low expected rates.

Increasing the duration of the retrospective tests to the 57 most recent years of the CPTI catalog (1950-2007) yields 83 earthquakes and leads to similar results but with greater statistical significance (Figure 6). In addition to the rejections mentioned in the preceding paragraph, the N-test now unequivocally rejects the forecasts Akinci-et-al.Hazgridx and Nanjo-et-al.RI, even when the confidence bounds of an NBD are considered. The conditional L-test rejects the forecast Meletti-et-al.MPS04 because of a likelihood score of negative infinity (discussed in section 6.2.2). The S-test results show that the forecast Nanjo-et-al.RI can be rejected in addition to the ALM forecasts and Meletti-et-al.MPS04. No forecasts can be rejected by the M-test, despite 57 years of data.

The longest period over which we evaluated the scaled forecasts was 106 years, spanning the full duration of the CPTI08 catalog and containing 183 earthquakes (Figure 7, see the online material for maps of the forecasts, likelihood ratios and concentration diagrams). The N-test results now show a clear separation between the group of forecasts that consistently overpredict, the forecast Nanjo-et-al.RI, which underpredicts, and the forecasts that cannot be rejected by assuming confidence bounds based on either a Poisson or a negative binomial distribution. Application of the conditional L-test additionally rejects the forecasts Nanjo-et-al.RI and Werner-et-al.CSI, while the S-test now also fails Werner-et-al.CSI and Zechar-Jordan.CSI.

Interestingly, four forecasts fail the M-test: Akinci-et-al.Hazgridx, Chan-et-al.HzaTI, Meletti-et-al.MPS04 and Nanjo-et-al.RI. In Figure 8, we compare the observed with their predicted magnitude distributions. For reference, we added a pure Gutenberg-Richter (GR) distribution with b-value equal to one, which passes the M-test. The magnitude distributions predicted by Akinci-et-al.Hazgridx and Nanjo-et-al.RI are close to exponential, but with b-values larger than one. As a result, large earthquakes are less likely, and the forecasts are penalized for the occurrence of three $m>7$ earthquakes. The magnitude distribution of the forecast Chan-et-al.HzaTI seems to reflect its non-parametric kernel estimation method (see section 2) and also underpredicts the rate of large shocks. Finally, the magnitude distribution of Meletti-et-al.MPS04 is non-monotonic: several characteristic magnitude bulges can be seen. However, the largest events occur between the bulges, for which the forecast is penalized.

The assumption of Poisson rate variability in the CSEP-Italy experiments (as well as other CSEP experiments, including RELM [Field, 2007; Schorlemmer et al., 2007]) has certain advantages. In particular, this is a simplifying assumption: because the Poisson distribution is defined by a single parameter, the forecasts do not require a complete probability distribution in each bin. Moreover, Poisson variability has often been used as a reference model against which to compare time-varying forecasts, and it yields an intuitive understanding.

Despite these advantages, however, this assumption is questionable, and the method of forcing each forecast to be characterized by the same uncertainty is not the only solution [see also the discussion by Zechar et al., 2010]. Werner and Sornette [2008] remarked that most forecast models generate their own likelihood distribution, and this distribution depends on the particular assumptions of the model; moreover, there is no reason to force every model to use the same form of likelihood distribution. The effect of this forcing is likely stronger for time-dependent, e.g. daily forecasts [Lombardi and Marzocchi, 2010], and it is difficult to judge (without the help of modelers) the quality of approximating each model-dependent distribution by a Poisson distribution. On the other hand, one can check whether or not the Poisson assumption is appropriate with respect to observations. In section 4.3, we showed that the target earthquake rate distribution is approximated better by an NBD than by a Poisson distribution. Therefore, time-independent forecasts that predict a Poisson rate variability necessarily fail more often than expected at $95\%$ confidence because the observed distribution differs from the model distribution. To improve time-independent forecasts, the (non-Poissonian and potentially negative binomial) marginal rate distribution over long timescales needs to be estimated. However, the parameter values of the rate NBD change as a function of the temporal and spatial boundaries of the study region over available observational periods [Kagan, 2010]. Whether a stable asymptotic limit exists (loosely speaking, whether seismic rates are stationary) remains an open question. For time-dependent models, on the other hand, several classes exist which are capable of producing a rate NBD over finite time periods including branching processes [Kagan, 2010] and Poisson processes with a stochastic rate parameter distributed according to the Gamma distribution.

Despite this criticism, it is unlikely that the Poisson distribution would be replaced by a model-dependent distribution that is substantially different, particularly for long-term models. Therefore, the p-values of the test statistics used in the N, L, S and M-tests might be biased towards lower values, but they do provide rough estimates. Nevertheless, one should bear in mind that a quantile score that is outside the $95\%$ confidence bounds of the Poisson distribution may be within the acceptable range if a model-dependent distribution were used. As an illustration, and to explore the effect of the Poisson assumption in these experiments, we created a set of modified forecasts with rate variability estimated from the observed history. The width of the $95\%$ confidence interval of the total rate forecast increased, in certain cases substantially. Several forecasts are rejected if a Poisson variability is assumed, while they pass the test under the assumption of an NBD. Overall, however, the p-values (quantile scores) of the test statistics based on the Poisson approximation often give good approximate values. Only in borderline cases did the Poisson assumption lead to (potentially) false rejections of forecasts.

The modified forecasts based on an NBD are not an entirely satisfactory solution to the problem. First, the model distribution in each bin should arise naturally from a model’s hypotheses, rather than an empirical adjustment made by those evaluating the forecast. Second, even if a negative binomial distribution adequately represents the distribution of the total number of observed events in an entire testing region, one should specify the parameter values for each bin to make the non-Poisson forecasts amenable also to the L-, S and M-tests. Therefore, future experiments should allow forecasts that are not characterized by Poisson rate uncertainty.

More generally, future experiments might consider other forecast formats and additional model classes. For example, stochastic point process models provide a continuous likelihood function which can characterize conditional dependence in time, magnitude and space (and focal mechanisms, etc.). As a result, full likelihood-based inference for point processes and tools for model-diagnostics are applicable to this class of models [e.g. Ogata, 1999; Daley and Vere-Jones, 2003; Schoenberg, 2003]. However, when considering new classes of forecasts, one may keep in mind that a major success of the RELM and CSEP experiments was the homogenization of forecast formats to facilitate comparative testing.

We explored results from the N, L, S and M-tests in this study because they are the “staple” CSEP tests. Other metrics for evaluating forecasts should certainly be considered, especially with regard to alarm-based tests [e.g. Molchan and Keilis-Borok, 2008; Zechar and Jordan, 2008] and further conditional likelihood tests [Zechar et al., 2010]. Overall, the N, L, S and M-tests are intuitive and relatively easy to interpret. However, we demonstrated a weakness in the L-test and replaced it with a conditional L-test that better assesses the quality of a forecast [see also Werner et al., 2010a]. Among the metrics, the S-test results were the most helpful in tracking down the weak features of forecasts, because the biggest differences between time-independent models lie in their spatial forecasts.

The M-test results were largely uninformative. Because the magnitude distributions considered here were so similar, this result is not surprising; indeed, it is in accordance with the statistical power exploration of Zechar et al. [2010]. No forecast could be rejected for target periods ranging from 5 to 57 years. Different tests, such as a traditional Kolmogorov-Smirnov (KS) test, should be compared with the current likelihood-based M-test, particularly in terms of statistical power.

The current status quo in CSEP experiments is to reject a forecast if it fails a single test at $95\%$ confidence. As we discussed above, the actual p-values provide a more meaningful assessment than a simple binary yes/no statement because the assumed confidence bounds may not accurately represent the model uncertainty. Furthermore, as the suite of tests grows, we should be concerned with joint confidence bounds of the ensemble of tests, rather than the individual significance levels of each test. Joint confidence bounds can be obtained from model simulations. A global confidence bound for the multiple tests can then be established. A similar question will arise when forecasts from the same model are tested within nested regions, as will be the case when considering the performance of a model’s forecast for Italy with that for the entire globe.

Finally, future experiments may consider developing tests that address particular characteristics of a forecast [see also the discussion by Zechar et al., 2010]. For example, a forecast might be a reflection of the hypothesis that the magnitude distribution varies as a function of tectonic setting. In this context, an M-test conditioned on the spatial distribution of observed earthquakes would provide a sharper test.

A summary of all results can be found in Tables 2 and 3. The Poisson N-test is possibly the strictest test within the present context, because none of the forecasts pass every N-test of the different periods. On the other hand, five forecasts pass all the N-tests with confidence bounds based on a negative binomial distribution (Gulia-Wiemer.ALM, Gulia-Wiemer.HALM, Werner-et-al.CSI, Werner-et-al.Hybrid and Zechar-Jordan.CSI). As we mentioned, several modelers indicated to us that their forecasts were calibrated on the moment magnitude scale. As a result, it is difficult to interpret their overpredictions beyond the obvious statement that the forecasts were poorly calibrated. The forecast Nanjo-et-al.RI is the only forecast that expects substantially fewer earthquakes than the observed sample mean, although the forecast fails the NBD N-test only for the longest of the considered target periods. Those forecasts that expect a number of shocks equal to the sample mean over their calibration period predict the number of quakes well, as should be expected.

With one important exception, the conditional L-test results largely reflect the S-test results, because the predicted magnitude distributions were consistent with observations from all but the 106-year target period. The exception concerns the occurrence of an earthquake in a space-magnitude bin in which an earthquake should have been impossible according to the forecast: the 1968 $M_{L}=6.39$ Belice earthquake happened in a spatial cell in which the forecast Meletti-et-al.MPS04 set a maximum magnitude of $M_{L}=6.25$. The discrepancy might be explained by the wrong magnitude conversion that the authors adopted and/or it may suggest that the model’s assumptions regarding the spatial variation of maximum magnitudes may need to be revised. However, if we had tested the forecast against the moment magnitude of the Belice earthquake ($M_{W}6.33$, according to the CPTI catalog), the forecast would have still failed, thus pointing towards the latter explanation.

The S-test results provided the most insight into the weaknesses of the forecasts. Only five forecasts pass all S-tests (Akinci-et-al.Hazgridx, Chan-et-al.HzaTI, Werner-et-al.Hybrid, Zechar-Jordan.CPTI and Zechar-Jordan.Hybrid). These forecasts fit the spatial distribution of the CSI and CPTI catalogs well, although they might overfit and perform poorly in the future. The models are also among the simplest, especially when compared to the forecast Meletti-et-al.MPS04. However, the forecasts Werner-et-al.CSI and Zechar-Jordan.CSI, which were calibrated on CSI data, cannot adequately forecast the spatial locations of earthquakes during the time period before the CSI data begins. This might indicate that the models are not smooth enough and do not anticipate sufficiently quiet regions becoming active.

The ALM group of forecasts (Gulia-Wiemer.ALM, Gulia-Wiemer.HALM and Schorlemmer-Wiemer.ALM) consistently fail the S-tests, and often perform worse than a uniform forecast, because isolated earthquakes occur in extremely low-probability “background” bins that cover roughly 50$\%$ of the region. We could not identify a common characteristic among the earthquakes that occurred in background bins. The incurred likelihood losses cannot be compensated by the gains achieved by adequately forecasting the majority of earthquakes. The results suggest that the ALM forecasts are overly optimistic in ruling out earthquakes in their background bins, i.e. the models are not smooth enough.

The forecast Meletti-et-al.MPS04 also often fails the S-test because of a minority of earthquakes that occur in low-probability regions. Almost all earthquakes that incur likelihood losses are located offshore. But while the forecast performs substantially better onshore, a few surprising onshore earthquake locations remain. Poor performance of a forecast for offshore earthquakes potentially raises the problem of the “weight” of each earthquake in the testing procedure. Specifically, if a model is intended for the practical purpose of seismic hazard assessment, then a rejection of its forecast due to offshore earthquakes may not have the same importance as a rejection due to earthquakes in regions of higher exposure and/or vulnerability.

Eight forecasts pass all M-tests (Gulia-Wiemer.ALM, Gulia-Wiemer.HALM, Schorlemmer-Wiemer.ALM, Werner-et-al.CSI, Werner-et-al.Hybrid, Zechar-Jordan.CPTI, Zechar-Jordan.CSI and Zechar-Jordan.Hybrid). Five of them are based on a simple Gutenberg-Richter distribution with uniform b-value equal to one (Werner-et-al.CSI, Werner-et-al.Hybrid, Zechar-Jordan.CPTI, Zechar-Jordan.CSI and Zechar-Jordan.Hybrid). This suggests that the hypothesis of a universally applicable, uniform Gutenberg-Richter distribution with b-value equal to one [e.g. Bird and Kagan, 2004] cannot be ruled out for the region of Italy.

Four forecasts fail the M-test during the 1901-2007 target period of the CPTI catalog. The magnitude distributions of the forecasts Akinci-et-al.Hazgridx, Nanjo-et-al.RI, Chan-et-al.HzaTI and Meletti-et-al.MPS04 do not adequately forecast the largest magnitudes and the three observed $M_{L}>7$, in particular. In the case of the Akinci-et-al.Hazgridx and Nanjo-et-al.RI forecasts, the reason seems to be a b-value of the Gutenberg-Richter distribution that is too large. The non-parametric estimate of Chan-et-al.HzaTI also decays too quickly. The magnitude distribution of Meletti-et-al.MPS04 reveals several characteristic magnitude values of elevated rates, but earthquakes also occur between them in extremely low-probability bins. However, these results should be interpreted cautiously because the same magnitude forecasts pass the 1950-2007 period, and because the greater uncertainty of the data prior to 1950 arguably influences the results.

The initial submission deadline for long-term earthquake forecasts for CSEP-Italy was 1 July 2009. Because the formal experiment was not intended to start until 1 August 2009, there was a brief period for initial analysis and quality control of the submitted forecasts. We provided a quick summary of the features of the forecasts and preliminary results of a retrospective evaluation to the modelers during this period. As a result, six of the eighteen time-independent and time-dependent long-term forecasts were modified and resubmitted before the final deadline on 1 August 2009. This initial quality control period was therefore extremely useful, and future experiments might consider expanding and formalizing the initial quality control period.

The short one-month period was, however, too short to evaluate the forecasts retrospectively in the detail we present here. During the course of this study, the problem of the wrong magnitude scaling was discovered. Because at least 4 of the 18 forecasts are affected, a second round of submissions was solicited for 1 November 2009, and 15 revisions (and 2 new forecasts) were submitted. This suggests that the feedback provided to modelers based on the present study was useful and informative. The task of converting even a relatively simple hypothesis into a testable, probabilistic earthquake forecast should not be underestimated, and we suggest that future experiments include some form of retrospective testing prior to final submission.

The retrospective evaluation also showed that at least the time-independent forecasts can be evaluated in a meaningful manner and that useful information about the models can be extracted. Such information is critical for the development of better forecasts and for the evaluation of the underlying hypotheses of earthquake occurrence.

At the same time, retrospective evaluation cannot replace the prospective tests with zero degrees of freedom. Given the relative robustness of the results from the retrospective evaluation, we anticipate that the prospective experiment will provide further useful and more definite information about the quality of the forecasts. Most importantly, if the second round of forecast submissions contains significantly improved forecasts with fewer technical errors, we expect to see real progress in our understanding of earthquake predictability.