CC BY
19
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 10, ES2005, doi:10.2205/2007ES000251, 2008
Prediction of extreme events: Fundamentals and prerequisites of verification
V. G. Kossobokov1,2 and A. A. Soloviev1
Received 28 January 2008; accepted 9 February 2008; published 29 February 2008.
[1] In many cases extreme events of different nature induce catastrophic consequences.
Therefore, in each case prediction of them is a long-living challenging problem of extremely high stakes. With a break-through in informatics many data relevant to catastrophic extremes became available for intensive search and testing of empirical “precursors”, as well as of conceptual hypotheses, thus, creating a fertile land for pattern recognition technique. Here we present the results of application of the same, perhaps, the simplest methodology to geophysical and socio-economical systems. Specifically, we (i) demonstrate the achievements of the on-going global monitoring of seismic activity aimed at prediction of the great and major earthquakes worldwide, which accommodates more than 15 years of real-time experience, and (ii) describe in more detail the quantitative experimentation in finding precursors of starts and ends of economic recessions, episodes of a sharp increase in the unemployment rate, and surges of homicides in a mega-city. Index Terms: 0500
Computational Geophysics; 1721 History of Geophysics: Nonlinear geophysics; 7223 Seismology: Earthquake interaction, forecasting, and prediction; 7290 Seismology: Computational seismology; KEYWORDS: seismology, earthquake prediction, model and theoretical seismicity.
Citation: Kossobokov, V. G. and A. A. Soloviev (2008), Prediction of extreme events: Fundamentals and prerequisites of verification,
Russ. J. Earth. Sci, 10, ES2005, doi:10.2205/2007ES000251.
1. Introduction
[2] Usually, prediction of extreme events is quite a difficult problem. By definition, an extreme event occurs rarely in a sequence of kindred phenomena that, usually, implies investigating a small sample of case-histories with a help of delicate statistical methods and data of different quality, collected in various conditions. Many of extreme events cluster and/or have self-similar distribution in space-time that, evidently, contradicts with a typically accepted simplified model of random occurrence. Such situation complicates additionally search for and definition of precursors, which could be used effectively in a prediction method. In the frames of objectivism’s viewpoint on probability it is not possible to give quantitative and/or probabilistic claims of the efficiency of a method for prediction of extreme events without a long series of its successes and failures-to-predict that, in turn, is impossible without its long enough testing by prediction determined in real time. Statistics of the ratio
1 International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Moscow, Russia
2Institute de Physique du Globe de Paris, Paris, France
Copyright 2008 by the Russian Journal of Earth Sciences.
ISSN: 1681-1208 (online)
of the number of failures to the total of successes and failures and the relative measure of the space-time volume of alarms, obtained during such testing, is necessary and sufficient for the assessment of reliability and potential of a method as of a prediction instrument, as well as it provides basic information for its improvement. Let us note that potential of usage is problem specific, i.e., it depends on a problem, and requires knowledge of a specific cost-and benefit function for the choice of an optimal strategy of prediction.
[3] These simple fundamentals are illustrated here with examples of the on-going prediction of extreme events in geophysical and socio-economical systems. Each system is considered as a complex hierarchical dissipative one and, evidently, possesses an important feature in common: persistent reoccurrence of extremes, i.e., abrupt overall changes, interpreted here as “critical transitions”. The following critical transitions are included: (i) large earthquakes in geophysical systems of the lithosphere blocks-and-faults, (ii) starts and ends of economic recessions, (iii) episodes of a sharp increase in the unemployment rate, (iv) surge of the homicides in socio-economic systems.
[4] These studies are based on a heuristic search of phenomena preceding critical transitions. We use the methodology of pattern recognition of infrequent events developed for studying rare phenomena of highly complex origin that, by their nature, limit the possibilities of using classical statistical or econometric methods. Our goal is to identify by an-
Figure 1. Possible outcomes of prediction.
alyzing the observable quantitative integrals and indicators the robust and unambiguously defined prediction algorithms of the “yes or no” variety. Given the values of integrals and indicators available to a given date, each of the algorithms provides unambiguous answer to the question whether a critical transition should be expected in the next time interval.
[5] Specifically, in terms of pattern recognition, an algorithm (a “recognition rule”) solves uniquely the following problem:
[6] given the time series of certain relevant indicators prior to a moment of time t,
[7] predict whether the critical transition will or will not occur during the time period (t, t + t).
[8] If the prediction is “yes”, the period (t, t + t) is the “period of alarm.”
[9] The possible outcomes of such a prediction are illustrated in Figure 1. Such “yes or no” prediction of specific extraordinary phenomena is different from predictions in a more traditional sense, i.e., extrapolation of a process in time, which is better supported by conceptual classical theories.
[10] The probabilistic component of such prediction is represented by the estimated probabilities of alarms on one side and failures to predict on the other. The probabilistic component is inevitable when a highly complex non-stationary process is considered, since the predictability of an originating non-linear dynamical system is limited in principle.
[11] If the algorithm is found and validated, it may be used in the two ways: (i) as a quantitative and reproducible description of phenomena premonitory to the critical transition that would provide empirical constraints for the theoretical modeling of the relevant process and (ii) as a practical tool complementing the existing methods of prediction of the critical transition.
[12] The methodology of the studies described below could be attributed to the so-called “technical” analysis, consisting of a heuristic search for phenomena preceding critical transitions. The alternative would be a “fundamental” analysis, focusing on “cause-and-effect” mechanisms leading to a critical transition under consideration. Regretfully, our knowledge of “cause-and-effect” mechanisms might be too concep-
tual and far from the nature of phenomenon under study. Therefore, the methodology used here is the pattern recognition of infrequent events developed by the artificial intelligence school of the Russian mathematician I. M. Gelfand [Gelfand et al., 1976] for the study of rare phenomena of highly complex origin. A distinctive feature of this methodology is the robustness of the analysis, which helps to overcome both the complexity of the process considered and the chronic imperfection of the data; in that aspect it is akin to exploratory data analysis, as developed by the statistics school of J. Tukey [Tukey, 1977]. Robust analysis - “a clear look at the whole” - is imperative in a study of any complex system [Gell-Mann, 1994]. The surest way not to predict such a system is to consider it in too fine detail [Crutchfield et al., 1986].
[13] The approach differs from but complements classical statistical and econometric methods such as regression analysis and ARIMA [Engle and McFadden, 1994]; see also Stock and Watson [1989], Klein and Niemira [1994], and Mostaghimi and Rezayat [1996] usually applied in socio-economical research. For comparison of this approach with the multiple regression analysis see Keilis-Borok et al. [2000] where a preliminary algorithm for the prediction of economic recessions is developed. Multiple regression analysis is not applied but the linear regression coefficients for single macroeconomic indicators are used to approximate the trends.
[14] The pattern recognition approach has been successfully applied to prediction in seismology and earthquake prediction [e.g. Gelfand et al., 1976; Keilis-Borok and Press, 1980; Keilis-Borok and Soloviev, 2003; Press and Allen, 1995; Press and Briggs, 1975], geological prospecting [e.g. Press and Briggs, 1977] as well as of the outcome of American elections [Keilis-Borok and Lichtman, 1993; Lichtman and Keilis-Borok, 1989] and in many other fields, as given in the references in these papers. Here the simplest version of such a methodology, called the “Hamming distance” [Gvishiani and Kosobokov, 1981; Keilis-Borok and Soloviev, 2003; Lichtman and Keilis-Borok, 1989; and refs. therein] is used. It is applied for classification of binary vectors into two classes on the basis of learning samples, which are used to determine a reference binary vector (kernel) with components that are typical for one class. The binary vectors under consideration are classified as belonging to this class if Hamming distance from them to the kernel is not greater than a certain threshold.
2. Geophysical Systems
[15] In this section we consider prediction of large earthquakes that occur in seismically active systems of the lithosphere blocks-and-faults.
2.1. Are Earthquakes Predictable?
[16] The temporal predictability of large earthquake occurrences requires a special comment on the recently re-
Table 1. Classification of earthquake predictions
Temporal, in years Spatial, in source zone size L
Long-term 10 Long-range Intermediate-term 1 Middle-range Short-term 0.01-0.1 Narrow Immediate 0.001 Exact Up to 100 5-10 2-3 1
vived discussions [Cyranoski, 2004; Geller et al., 1997; Wyss, 1997], (Nature Debates, 1999, http://www.nature.com/na-ture/debates/earthquake/equake_frameset.html). No current theory of dynamics of seismic activity can answer this question. Inevitably, a negative statement that asserts a non-trivial limitation on predictability is merely a conjecture. On the other hand, forward testing of a reproducible prediction method and, so far, in no other way, can unequivocally establish a certain degree of predictability of earthquakes. The results of the on-going real-time monitoring of the global seismic activity aimed at intermediate-term middle-range prediction of the largest earthquakes (http://www.mitp.ru) has proved [Kossobokov and Shebalin, 2003; Kossobokov et al., 1999] the high statistical significance of the two methods, algorithms M8 [Keilis-Borok and Kossobokov, 1990a] and MSc [Kossobokov et al., 1990], which short descriptions are given below, did confirm a positive statement on predictability of earthquakes. Furthermore, it appears that in some cases the inverse cascading of seismic activity to a catastrophe evolves through long-, intermediate-, short-, immediate-term and even nucleation [Ellsworth and Beroza, 1995] phases.
[17] Following common perception many investigators usually overlook spatial modes of predictions concentrating their efforts on predicting the “exact” fault segment ready to rupture (e.g., the Parkfield earthquake prediction experiment), which is by far a more difficult and might be an unsolvable problem. Being related to the rupture size L = L(M) of the incipient earthquake of magnitude M, such modes could be summarized in a classification of location of a source zone from a wider prediction ranges (Table 1).
[18] From a viewpoint of such a classification, the earthquake prediction problem is naturally approached by a hierarchical, step-by-step prediction technique, which accounts for multi-scale escalation of seismic activity to the main rupture [Keilis-Borok, 1990]. Table 1 disregards term-less pre-
dictions although identification of earthquake-prone areas, e.g., by pattern recognition methods [Gorshkov et al., 2003], deliver a zero-approximation for a target earthquake location. Moreover, the Gutenberg-Richter law suggests limiting magnitude range of prediction to about one unit, i.e., M0 < M < M0 + AM and AM < 1. Otherwise, the statistics would be related to dominating smallest earthquakes and, therefore, attributing it to much larger events is misleading.
[19] The on-going real-time monitoring of the global seismic activity aimed at intermediate-term middle-range prediction of the largest earthquakes has a 15-year history now [Healy et al., 1992; Kossobokov et al., 1999]. Tables 2 and 3 give the up-to-date summary of the prediction outcomes and prove certain predictability of the great and major earthquakes beyond any reasonable doubt (the achieved confidence is above 99%).
[20] It is notable that to drive the achieved confidence level below 95%, the real-time monitoring should fail to predict the next six M8.0+ or nineteen M7.5+ events in a row, which seems unlikely. The results require special comments in the following sections. Since the estimates presented in the tables use the most conservative measure of the alarm volume accounting for empirical distribution of epicenters, called measure ^ below, we describe it first, and then explain what stand behind M8 and MSc and their global and regional testing.
2.2. How to Measure Space Occupied by Seismic Activity?
[21] Are the results of the earthquake prediction experiment better than the random guessing or they are not? A statistical conclusion about that could be attributed in the following general way:
[22] Let T and S be the total time and territory considered; At is the territory covered by the alarms at time t; t x ^ is a measure on T x S (we consider here a direct product measure t x ^ reserving a general case of a time-space dependent measure v for the future more sophisticated null-hypotheses); N counts the total number of large earthquakes with M > M0 within T x S and n counts how many of them are predicted. The time-space occupied by alarms,
Table 2. Worldwide performance of earthquake prediction algorithms M8 and M8-MSc: Magnitude range M8.0+
Large earthquakes Percentage of alarms Confidence level, %
Test period Predicted by Total
M8 M8-MSc M8 M8-MSc M8 M8-MSc
1985-2007 12 9 17 32.93 16.78 99.83 99.93
1992-2007 10 7 15 29.17 14.54 99.71 99.70
Table 3. Worldwide performance of earthquake prediction algorithms M8 and M8-MSc: Magnitude range M7.5+
Large earthquakes Percentage of alarms Confidence level, %
Test period Predicted by Total
M8 M8-MSc M8 M8-MSc M8 M8-MSc
1985-2007 32 16 57 30.27 9.79 99.99 99.99
1992-2007 22 10 45 24.29 8.79 99.97 99.50
A = y At, in percentage to the total space-time considered T
equals
jd(T x ^)
A
P = --------------- .
I d(T x ^)
T XS
[23] By common definition the two dual levels of statistical significance and confidence of the prediction results equal to
a = 1 — B(n — 1, N, p)
and
1 — a = B(n — 1, N, p) ,
where B is the cumulative binomial distribution function. The smaller is the significance level a, the larger is the confidence level 1 — a and the higher is significance of predictions under testing.
[24] When testing temporal predictability of earthquakes it is natural to make the following choice of the product measure t x ^: the uniform measure t, which corresponds to the Poisson, random recurrence of earthquakes and the measure ^ proportional to spatial density of epicenters. Specifically, determine the measure ^ of an area proportional to the number of hypo- or epicenters of earthquakes from a sample catalog, for example, earthquakes above certain magnitude cutoff Mc. This empirical spatial measure of seismic distribution is by far more adequate than the literal measures of volume in km3 or territory in km2 for estimating statistical significance of the prediction results. Evidently, the literal measures of volume or territory equalize the areas of high and low seismic activity, at the extreme, the areas where earthquake happen and do not happen.
[25] The actual, empirical distribution of earthquake locations is the best present day knowledge estimate of where earthquakes may occur. The recipe of using the ^-measure and counting p is the following: Choose a sample catalog. Count how many events from the catalog are inside the volume or the territory considered; this will be your denominator. At a given time, count how many events from the catalog are inside the area of alarm; this will be your numerator. Integrate the ratio over the time of prediction experiment. This is the exact way of computing Percentage of alarms and Confidence level in Table 2 (where the catalog sampled all earthquakes of magnitude 4 or larger from the NEIC Global Hypocenter’s Data Base in 1963-1984).
[26] This simple recipe has a nice analogy, called Seismic Roulette, that justifies using statistical tools available since Blaise Pascal (1623-1662):
• Consider a roulette wheel with as many sectors as the number of events in a sample catalog, a sector per each event;
• Make your bet according to prediction: determine
which events are inside area of alarm, and put one chip in each of the corresponding sectors;
• Nature turns the wheel.
[27] If you play seismic roulette systematically, then you win and lose systematically. If the roulette is not perfect and you are smart enough to choose an effective strategy, then your wins will outscore loses! There is evident option of switching to “antipodal” strategy [Keilis-Borok and Soloviev, 2003; Molchan, 1994] when the losses outscore wins. The results of the global test of the algorithms M8 and MSc did confirm “imperfection” of Nature in recurrence of the great and major earthquakes and suggest using it for the benefit of the population exposed to seismic hazard.
2.3. The M8 and MSc Algorithms
[28] Both algorithms are reproducible earthquake prediction methods that satisfy the consensus definition [Allen et al., 1976] and make use of seismic activity reported in routine seismic catalogs. The M8 is applied first. It scans the territory in question for the areas in alarm (Figure 2), so-called Time of Increased Probability, TIP. The MSc is applied to reduce the area of alarm by analyzing dynamics at lower magnitude levels of seismic hierarchy. Sometimes, the data is enough to get a near-perfect outline of the incipient large earthquake. More often the catalog of earthquakes is exhausted already at the M8 analysis and the prediction remains in the middle range.
[29] The M8 intermediate-term earthquake prediction algorithm was designed by retroactive analysis of dynamics of seismic activity preceding the greatest, magnitude 8.0 or more, earthquakes worldwide, hence its name. Its prototype [Keilis-Borok and Kossobokov, 1984] and the original version [Keilis-Borok and Kossobokov, 1987] were tested retroactively at recorded epicenters of earthquakes of magnitude 8.0 or greater from 1857-1983. Figure 3 shows, as an example of the M8 prediction in the real time, the case-history of
Figure 2. General scheme of applying reproducible earthquake-prediction algorithm: Areas of investigation overlay seismic region; seismic sequences in each area gives reproducible description of the present state, which is then used to diagnose an alert, so-called time of increased probability, TIP.
the 4 June 2000, Ms8.0 Sumatra earthquake: The Andaman-Sumatra-Java segment of the global prediction map issued in January 2000 along with epicenters of the great main shock and its first aftershocks are given on the left, while, on the right, the figure depicts the space-time diagram of seismic activity in the circle of investigation (Test no. 34) with radius 667 km where the alarm was in progress when the great earthquake happened and below it presents the functions of algorithm M8 with their abnormal values marked by heavy dots. The arrows indicate the great shock, and small circles stand for smaller magnitude earthquakes used by the algorithm for determining the alarm. The distance along the seismic belt measured in kilometers from the center of the circle is plotted on the vertical axis. Time is plotted along the horizontal axis.
[30] The algorithm M8 uses traditional description of a dynamical system adding to a common phase space of rate (i.e. number of mainshocks, N) and rate differential (i.e., deviation of N from a longer-term average, L) the dimensionless
concentration (i.e., the average source size divided by the average distance between sources, Z) and a characteristic measure of clustering (i.e., maximum number of aftershocks, B). The analysis of seismic activity in one region may distinguish a number of magnitude ranges and deliver a hierarchy of predictions [Keilis-Borok and Kossobokov, 1990b].
[31] The algorithm recognizes criterion, defined by extreme values of the phase space coordinates, as a vicinity of the system singularity. When a trajectory enters the criterion, probability of extreme event increases to the level sufficient for effective provision of a catastrophic event (evidently, the ranges of extreme values of the M8 algorithm functions define the kernel used for classification of times by the Hamming distance). The exact definitions and computer code of the M8 algorithm are published [Healy et al., 1992; Keilis-Borok and Kossobokov, 1990a; Kossobokov, 1997].
[32] Retrospectively the standard version of the algorithm [Keilis-Borok and Kossobokov, 1990a] was applied to predict earthquakes with magnitudes from above 8.0 to 4.9 in
Figure 3. Global testing of algorithms M8 and MSc, M0 = 8.0 [Kossobokov and Shebalin, 2003]: The 4 June 2000 Sumatra earthquake. Observe the highlighted circular areas of alarm in the first approximation determined by algorithm M8 and the highlighted rectangular areas of alarm in the second approximation determined by algorithm MSc. A foreshock of magnitude 4.7 (diamond) occurred within a day in advance of the great shock.
a number of regions worldwide. Its modified versions apply also in regions of seismic activity lower than required by the original version [Bhatia et al., 1989; Gahalaut et al., 1992; Kossobokov et al., 1989; Latoussakis and Kossobokov, 1990; Peresan et al., 2005; Romachkova et al., 1998].
[33] The second approximation prediction method MSc [Kossobokov et al., 1990] was designed by retroactive analysis of the detailed regional seismic catalog prior to the Eureka earthquake (1980, M=7.2) near Cape Mendocino in California, hence its name, Mendocino Scenario, and an abbreviation. Qualitatively, the MSc algorithm outlines such an area of the territory of alarm where the activity, from the beginning of seismic inverse cascade recognized by the first approximation prediction algorithm (e.g. by M8), is continuously high and infrequently drops for a short time. Such an alternation of activity must have a sufficient temporal and/or spatial span. The phenomenon, which is used in the MSc algorithm, might reflect the second (possibly, shorter-term and, definitely, narrow-range) stage of the premonitory rise of seismic activity near the incipient source of the main shock. In reduction of territorial uncertainty of the M8 predictions, the MSc algorithm outperforms by at least a factor of 2 a few simple alternatives like the earthquake-prone cells in the area of alarm or the most active cells that contain certain part of the recent seismic activity [Kossobokov et al., 1990].
2.4. Testing Earthquake Predictions
[34] After prediction of the Spitak 1988 and Loma Prieta 1989 earthquakes in real-time J. H. Healy, V. G. Kossobokov, and J. W. Dewey designed a rigid test to evaluate the M8 algorithm [Healy et al., 1992]. Since 1991 each half-year the algorithm has been applied in a real time prediction mode to monitor seismic dynamics of the entire Circum Pacific (that is the reason for distinguishing the two periods of testing in Table 2: since the design of the algorithm in 1985, and since the formal publication of the settings for global monitoring in 1992). More extended testing, for all seismically active territories on Earth where seismic data is enough to run the standard version of algorithm M8 was carried on in parallel [Kossobokov and Khokhlov, 1993; Kossobokov et al., 1992, 1999]. Unfortunately, testing in seismic regions of the Former Soviet Union where the rescaling of the original M8 algorithm was tested first in 1986 on the “Earthquakes in the USSR” catalog aimed mostly at M6.5+ earthquakes were discontinued due to the collapse of the state and some of its seismological structures. The testing included Vrancea, Caucasus, Turkmen territories, Pamirs and Tien Shan. The reestablishment of seismic monitoring aimed at prediction of large magnitude earthquakes in Vrancea, Caucasus and
Central Asia looks feasible nowadays, specifically after the development of a recently proposed scheme for the spatial stabilization of the intermediate-term middle-range predictions [Kossobokov et al., 2002]. The scheme, named M8S, makes use of the multiple application of the M8 algorithm in a large number of objectively distributed circles of investigation and aims at elimination of spatially sporadic alarms. In fact, it appears to guarantee a more objective and reliable diagnosis of times of increased probability and, at the same time, is less restrictive to input seismic data. At the moment it is used for the real-time monitoring of the Italian territory being aimed at M6.5+, M6.0+, and M5.5+ earthquakes [Peresan et al., 2005].
[35] In the Global Test aimed at M8.0+ earthquakes the algorithms M8 and MSc are applied in 262 overlapping circles of investigation, of which 170 scan near-uniformly Circum-Pacific and its surroundings, 92 circles taken from Alpine-Himalayan Belt and Myanmar (25 in Mediterranean, 25 in Asia Minor and Iran, 28 in Pamirs-Hindukush, and 14 in Myanmar). These cover about 8090% of the major seismic belts of the Earth. The complete set of predictions in 1985-2007 could be viewed at http://mitp.ru/predictions.html, although the access to those in progress is restricted. In general, the alarms last for about five years, but could expire before or extend beyond this limit under unusual local changes of seismic regime. On average the M8 alarms cover less than one third of the whole seismic territory considered, while MSc reduces this area by another factor of two (Table 2). The probability gain in confirmed predictions depends on locality and varies from 2-3 in regions of extremely high activity, like Tonga-Kermadek, to 20-100 in regions where recurrence of the great earthquakes is much lower than average, like southern Sumatra or Tibet.
[36] Aimed at M7.5+ earthquakes the algorithms are applied in 180 circles, which in total cover about 75% of the major seismic belts. 147 of them represent seismic regions of Circum Pacific, while the remaining 33 ones compose of 15 from Mediterranean, 4 from Iran, 11 from Pamirs-Hindukush, and 3 from Myanmar. On average the M8 alarms cover about one quarter (24.29% in accordance with measure ^ since 1992, Table 3) of the whole seismic territory considered, while MSc reduces this area below 10%. For this magnitude range, certain decay in performance is observed in the recent years. There are indications that this could be inflicted by the changes either in the global seismic regime or in reporting the magnitudes or both: (i) most of the failures-to-predict occurred during the unusual rise of seismic energy release, have magnitude below 7.75 and are thrust or normal faulting [Kossobokov et al., 1999]; (ii) starting form 1993 the NEIC changed the procedures of the global database compilation, substituting MS from Pasadena and Berkley with values of MW from Harvard and USGS.
2.5. Can Mega Earthquakes Be Predicted?
[37] The statistics given in Table 2 do not include the recent mega-earthquakes in Indonesia that are much stronger than M8.0+ events. Specifically, the size of the 26 December 2004, MW9.3 (Ms8.8) off the west coast of Northern
Sumatra Great Asian, Sumatra-Andaman mega-thrust and its follower the 28 March 2005, MW8.7 (Ms8.4) Nias earthquake, brings them out of the list of target earthquakes of the Global Test. (Note that the most recent great earthquakes on 12 September 2007 MW8.4 and 7.9 (Ms8.5 and 8.1) that stroke southern Sumatra at 1110 and 2349 GMT, correspondingly, were successfully predicted in course the Global Test.)
[38] First of all, the linear dimension of the source of the first one is about 1000-1300 km, i.e., about the diameter of circles of investigation used in the Global Test of M8 to predict M8.0+ earthquakes. The linear dimension of the second one is above 450 km. The source length of the M8.0+ events in 1985-2003, usually accounts to about 150-300 km. Therefore, since the logic of the methodology suggests the proportions of investigation about 5-10 times larger than the target earthquake size, it would be naive and ambiguous to expect a success of the monitoring aimed either at M8.0+ or M7.5+ earthquakes in predicting the 26 December 2004 and 28 March 2005 events. According to the M8 algorithm predictions we were not expecting any M8.0+ or M7.5+ events in the Indian Ocean neither during the second half of 2004 nor in the first half of 2005 and, in fact, these did not happen.
[39] On the other hand, if on 1 July 2004 someone, enough ambiguous to extend application of the M8 algorithm into unexampled magnitude range aiming at M9.0+ earthquakes, then he or she would have diagnosed Time of Increased Probability in advance of the 2004 Sumatra-Andaman mega-trust event. The genuine M8 computer code run with the target earthquake magnitude threshold equal to 9.0 and the radius of CI’s increased to 3000 km determines the current alarm.
[40] In fact, this is a unique unexampled confirmation that the algorithm, designed for prediction of M8.0+ earthquakes and tested in many applications rescaled for prediction of smaller magnitude earthquakes (e.g., down to M5.5+ in Italy, http://www.mitp.ru/m8s/M8s_italy.html), is applicable for prediction of the mega-earthquakes of M9.0+ (and M8.5+). Of course, we are not that ambiguous to go from the first indication to a routine prediction, but feel the 26 December 2004 case history very important for general understanding of the methodology and the Problem of Earthquake Prediction, in general.
[41] What is the extent in space and time of the M8 algorithm TIPs for M9.0+? The answer is thought provoking: in the 25 years of retrospective analysis of available data there was one cluster of TIPs in 1984-1989 around western Mediterranean (a compact union of the eight out of the 262 circles of investigation) plus another one in 1994-1999 around Cascadia plate (a compact union of the five circles of investigation off the western coast of U.S.), which produce no M9.0+ event. The union of TIPs to date has global extent: it encountered the maximum of 145 circles of investigation in 2003, 124 - by the time of the 2004 Sumatra-Andaman mega-trust; the 47 circles of investigation in alarm to date cover about one half of the global seismic belts (49.09% in accordance with measure ^). Having in mind the evidence, which suggests clustered occurrence of seismic events including mega-earthquakes, we cannot reject such a possibility of further confirmations in the nearest future.
2.6. Implications
[42] The M8 and MSc algorithms, which test is presented here, make use of seismic activation and the growing correlation of earthquakes at the approach of the Big One. The predictions could be done on the basis of earthquake catalogs routinely available in the majority of seismic regions. There are evident limitations in performance. With more complete catalogs and, hopefully, with other relevant data the areas of alarm may be substantially reduced in the second and, perhaps, further approximations at the cost of additional failures-to-predict.
[43] The M8 and MSc algorithms are neither optimal nor unique. Together with other methods [Harte et al., 2003; Keilis-Borok and Rotwain, 1990; Keilis-Borok et al., 1988; Shebalin et al., 2003; Vorobieva, 1999; etc.] they hallmark a break-through in earthquake prediction research that leads from term-less assessment of seismic hazard to reliable intermediate-term alert of increased probability. The accuracy could be improved in course of a systematic monitoring of the alarm areas and by designing a new generation of earthquake prediction technique of higher accuracy.
[44] Thus, the approach has already demonstrated efficiency of the pattern recognition technique in solving the earthquake prediction problem at global and regional scales and form the basis of Quantitative Earthquake Prediction. The achievements of pattern recognition in the design of the reproducible algorithms predicting large earthquakes and the verified statistical validity of their predictions confirm the underlying paradigms:
• Seismic premonitory patterns exist;
• Formation of earthquake precursors at scale of years involves large size fault system;
• The phenomena are similar in a wide range of tectonic environment;
• The phenomena are universal being observed in other complex non-linear systems.
[45] Seismic Roulette is not perfect. Therefore, the existing reliable predictions of limited accuracy could be used in a knowledgeable way to the benefit of population living in seismic regions. The methodology linking them to optimal strategies for disaster management exists and is rather developed [Molchan, 2003]. The intermediate-term middle-range accuracy is quite enough for undertaking earthquake preparedness measures, which would prevent a considerable part of damage and human loss, although far from the total.
[46] The predictions also provide reliable empirical constrains for modeling earthquakes and earthquake sequences. The prediction results evidence that distributed seismic activity is a problem in statistical physics. They favor the hypothesis that earthquakes follow a general hierarchical process that proceeds via a sequence of inverse cascades to produce self-similar scaling (intermediate asymptotic), which then truncates at the largest scales bursting into direct cascades [Gabrielov et al., 1999].
3. Socio-Economic Systems
[47] In this section we consider prediction of starts and ends of economic recessions, episodes of a sharp increase in the unemployment rate, and surges of the homicides in a mega-city. Prediction is given by a discrete sequence of alarms (Figure 1). Its accuracy is captured by statistics of false alarms (including their total space-time) and failures to predict.
[48] It have been found for the five recessions in the USA since 1962 to 1996 that each of them is preceded by a specific pattern of 6 economic indexes, which are defined at the lowest (binary) level of resolution. This pattern was present during 6 to 14 month before each recession and at no other time, suggesting a hypothetical prediction algorithm. The algorithm is exceedingly robust: the retrospectively diagnosed alarms remain about the same after variation of its adjustable numerical parameters, and of other non-unique decisions, involved in its determination. Another algorithm has been formulated for predicting the end of an American economic recession by means of analysis of the same macroeconomic indicators within the recession period. It indicates up to 6 months long time interval, during which recession will end. First application of the algorithm to out of sample data (not used in its development) is successful: it predicted that the last recession started in April 2001 would end between July and December 2001 and that recession indeed ended in November 2001.
[49] A specific “premonitory” pattern of three macroeconomic indicators that may be used for algorithmic prediction of FAUs has been found for unemployment in France between 1962 and 1997. Among seven FAUs identified within these years six are preceded within 12 months by this pattern that appears at no other time. The application of this algorithm to Germany, Italy and the USA yields similar results. The first advance prediction, for the USA for early 2000, has been successful.
[50] The analysis of statistics of several types of crime in Los Angeles over the period 1975-2002 focused on how these statistics change before a sharp and lasting rise (“a surge”) of the homicide rate. The goal was to find an algorithm for predicting such a surge by monitoring the rates of different crimes. The results may be summarized as follows: episodes of a rise of burglaries and assaults simultaneously occur 4 to 11 months before a homicide surge, while robberies decline. Later on, closer to the rise in homicides, robberies start to rise. These changes are given unambiguous and quantitative definitions, which are used to formulate a hypothetical algorithm for the prediction of homicide surges. The retrospective analysis shows that this algorithm is applicable through all the years considered despite substantial changes both in socio-economic conditions and in the counting of crimes. Moreover, it gives satisfactory results for the prediction of homicide surges in New York City as well. Sensitivity tests show that predictions are stable to variations of the adjustable elements of the algorithm.
[51] Decisive validation of these findings requires experimentation in advance prediction, for which these studies set up a base. Particularly encouraging for this further research
Table 4. American Economic Recessions, 1960-2003
# Peaks Troughs
1 April 1960 February 1961
2 December 1969 November 1970
3 November 1973 March 1975
4 January 1980 July 1980
5 July 1981 November 1982
6 July 1990 March 1991
7 March 2001 November 2001
Note: Peak is the last month before a recession, and trough is the
last month of a recession (a recession ends in this month).
is the wealth of yet untapped possibilities: we have used so far only a small part of the data and mathematical models that are currently available and that are relevant to dynamics of complex socio-economic systems.
3.1. Prediction Targets
Figure 5. Target of prediction - the Start of the Homicide Surge (“SHR”); schematic definition. Gray bar marks the period of homicide surge.
[52] Recessions. The first months of the recessions and the first months after them are considered as moments of critical transitions. These months are given by the National Bureau of Economic Research (NBER). Seven recessions occurred from January 1960 to April 2002 are listed in Table 4 according to the NBER data. When starts of the recessions are predicted the targets are the first months after the peaks. When ends of the recessions are predicted the targets are the first months after the troughs.
[53] Unemployment. A specific phenomenon in the dynamics of unemployment: a sharp increase in the rate of unemployment growth is considered as the prediction tar-
Time
Figure 4. Fast acceleration of unemployment (FAU): schematic definition. Thin line - monthly unemployment; with seasonal quasiperiodic variations. Thick line - monthly unemployment, with seasonal variations smoothed away. The arrow indicates a FAU - the sharp bend of the smoothed curve. The moment of a FAU is the target of prediction.
get. Qualitatively, this phenomenon is illustrated in Figure 4. The thin line is the monthly number of unemployed u(t), including seasonal variations. After smoothing u(t) to eliminate such variations, the function U(t) is obtained. The prediction target is the starting month of a strong and lasting increase in U(t). An example is the turning point indicated by the arrow in Figure 5. We call this target by the acronym FAU, for “Fast Acceleration of Unemployment.”
[54] Homicide statistics. A specific phenomenon in crime dynamics: a large and lasting increase in the homicide rate is considered as the prediction target. Qualitatively, this phenomenon is illustrated in Figure 5. The prediction target is the starting month of a large and lasting increase in the smoothed monthly homicide rate. It is indicated by the arrow in Figure 4. We call it by the acronym SHS, for “Start of the Homicide Surge.”
3.2. Common Notations
[55] Socio-economic systems are described by monthly series of indexes, e.g. industrial production, long- and shortterm interest rates, statistics of different crime types etc. Let f (m) is one of such series.
[56] Common notation here and below is Wf ^ -, p^ -
the local linear least-squares regression of a function f (m) within the sliding time window (q, p):
Wf ^-, p^ = Kf (q, p)l + Bf (q, p), q < l < p . (1)
It is assumed here that m, l, q, and p are integers that stand for consecutive numbers of months during a time period under consideration.
[57] To expose premonitory (relative a critical transition) changes in the behaviour of f (m) two following functions are introduced, depicting the trend by a linear least-squares regression (1) of f (m).
Kf
= Kf (m — s, m) .
(2)
This function that is a regression coefficient defined as in (1) approximates the trend of f (m) in a time window of s
months long, (m—s, m). The value of Kf (may be used for determination of a precursory pattern that appears in month m since it is attributed to the end of the time window where it is determined; accordingly it does not depend on information on the future (after month m).
q
Rf m = f(m) — Wf m, m — 1
q
(3)
This function depicts deviation of an index from its longterm extrapolation. Here the linear least-squares regression (1) is determined on a time interval (q, m — 1), and it is assumed that this interval is rather long, i.e. m — q is large.
[58] Let g(m) be function (2) or (3) or source series f (m) itself. If g(m) demonstrates premonitory behavior then a robust quantitative definition of this is given as follows. The values of g(m) are defined on the lowest level of resolution, distinguishing only the values above and below a threshold Tg(Q). It is defined as a percentile of a level Q, that is, by the condition that g(m) exceeds Tg (Q) during Q% of the months considered.
3.3. Start of Recession
[59] The time period from January 1959 to April 1996 has been initially considered [Keilis-Borok et al., 2000]. Six recessions occurred during this period; they are identified in Table 4 by the last month before a recession (“peak”) and the last month of a recession (“trough”). The time series, consisting of monthly values of the six indexes, are analyzed. These indexes that were already known [Stock and Watson, 1989, 1993], as correlated with the approach of a recession are listed below (abbreviations are the same, as in [Stock and Watson, 1993]):
[60] IP. Industrial Production, total: index of real (constant dollars, dimensionless) output in the entire economy. This represents mainly manufacturing because of the difficulties in measuring the quantity of output in services (services include travel agents, banking, etc.).
[61] LHELL. Index of “help wanted” advertising. This is put together by a private publishing company that measures the amount of job advertising (column-inches) in a number of major newspapers.
[62] LUINC. Average weekly number of people claiming unemployment insurance.
[63] INVMTQ. Total inventories in manufacturing and trade, in real dollars. Includes intermediate inventories (for example held by manufacturers, ready to be sent to retailers) and final goods inventories (goods on shelves in stores).
[64] FYGM3. Interest rate on 90 day U.S. treasury bills at an annual rate (in percent).
[65] G10FF. Difference between interest rate on 10 year U.S. Treasury bond, and federal funds interest rate, on annual basis.
[66] It has been found looking at the dependence of the indexes on time that trough the whole period the indexes IP, and INVMTQ have an evident upward trend; the indexes LHELL, LUINC, FYGM3, and G10FF have no such trend.
[67] Time periods. The time of each recession and 5 months after it were eliminated, since the behaviour of the indexes may have some special features during a recession and its aftermath; the interval 5 month is chosen, because subsequent recessions are officially considered as different ones, only if they are separated by at least six months. At this stage the period before the first recession (since it starts too close to the beginning of the data set) and the period after the sixth one were excluded also. Thus the following 5 periods were used for the search of premonitory phenomena.
Wi: August 1961-December 1969 (101 months);
W2: May 1971-November 1973 (31 months);
W3: September 1975-January 1980 (53 months);
W4: January 1981-July 1981 (7 months);
W5: May 1983-July 1990 (87 months).
The combination of these periods W = W1 U W2 U W3 U W4 U W5 contains 279 months. In all calculations the unit of time is a year and duration of a month is 1/12. Months are identified by their sequential numbers, i = 1, 2, . . . .
[68] Hypothesis. The index G10FF (Figure 6) seems to become unusually low before the recessions. Premonitory behavior of other indexes is not so clear but their temporal trend is changing, when a recession is approaching. To expose these changes functions Kf (—) (2) and
Rf (m) = Rf ^m ^ (3) with q being the first month after
the end (the through) of the previous recession are used.
[69] For further analysis besides the original index G10FF, the following 5 functions, which seem to have different prevailing values close to an incipient recession and far from it have been selected:
Functions: RIP(m), RINVMTQ ^, klhell
K L
m\ rf 10 > R
(m).
Abbreviations: IPR, INVR, LHK5, LUK10, FYG3R.
[70] The behaviour of these functions and of the index G10FF suggests a hypothesis that the recessions considered are often preceded by relatively large values of LUK10 and FYG3R and relatively small values of IPR, INVR, LHK5, and G10FF.
[71] Discretization. Following to Section 2.2. thresholds Tg (Q) have been determined for these six time series. The relevant values of Q are given in Table 5. Now the values g(m) are coded on the lowest level of resolution, 0 or 1, discriminating only the values g(m) < Tg (Q) and g(m) > Tg(Q). It will be convenient to give the same notation, say 1, to the values, which became more frequent, when a recession is approaching. Accordingly, 1 is specified
Figure 6. Index G10FF, the threshold of discretization is shown by a horizontal line, shaded vertical bars indicate recessions.
for the values g(m) > Tg(Q) in the case of the functions LUK10 and FYG3R, and for the values g(m) < Tg(Q) in the case of the four other time series.
[72] After the discretization the description of situation in each month has been reduced to a binary vector with 6 components. Each component has been defined in such a way, that (if the hypothesis is correct) the values “1” would became more frequent when a recession is approaching. Accordingly, the description of pre-recession situations would be close to the vector (1,1,1,1,1,1), which is called the kernel. Let D(m) be the number of zeros in a code of a month, a
- its Hamming’s distance from the kernel. One can assume that the approach of a recession is recognized by the small values of D(m). A priori this is not clear, in spite of the way the zeros are defined; for example, premonitory changes of the functions may appear not simultaneously, even if the above hypothesis is correct. But it has been found that the values D(m) < 2 are confined to 6 to 14 months, preceding each recession. The change of D(m) with time suggests the following prediction algorithm: an alarm is declared for
Table 5. Values of Q
g(m) Precursory value Q, %
IPR Low 75.0
INVR Low 25.0
G10FF Low 90.0
LHK5 Low 66.7
LUK10 High 16.7
FYG3R High 25.0
three months after each month with D(m) < 2 (regardless of whether this month belongs or not to an alarm which has been already declared). Three months of alarm after the month when D(m) < 2 are introduced for the following reasons: a premonitory phenomenon does not necessarily appear just before a recession; and it is preferable to merge the alarms, which are close to each other.
[73] The alarms thus defined are shown in Figure 7 by black bars. Note that each alarm extends also to first 2 or 3 months of a recession. One can see that all 5 recessions were preceded by continuous alarms. The longest one lasted 13 months, one alarm lasted 10 months, and 3 alarms
- 5 months. There were no false alarms. Total duration of alarms, 38 months, is 13.6% of the time covered by the analysis (W set). No recessions happened since March 1991; application of this algorithm to the subsequent years, September 1991-April 1996, does not give a false alarm.
[74] Advance prediction covered subsequent 11 years, up to now. The algorithm detected an alarm on May 2001, one month later than it started according to NBER (Table 4). It is too early to evaluate the rates of the errors. Note that in practice detection of the alarm can be delayed by the delay in updating the indexes.
3.4. End of a Recession
[75] The macroeconomic indicators listed above (Section 3.3.) were analyzed in order to find phenomena preceding the end of an American economic recession.
[76] The data concerning the first 6 recessions from Table 4 were analyzed initially. Thus the following 6 periods were
I'I'I'I'I'I'I'I'I'I'I
59 61 63 65 67 69
73
71
'I'I'I'I'I'r
75 77 79 81 83 85 87 89
w
91
93
'I'I'I'I
95 97
Figure 7. Alarms (shown by black bars) and recessions (shown by grey bars).
considered (for each recession a month before it (the peak) and its last month (the trough) are included).
W1: April 1960-February 1961 (11 months);
W2: December 1969-November 1970 (12 months);
W3: November 1973-March 1975 (17 months);
W4: January 1980-July 1980 (7 months);
W5: July 1981-November 1982 (17 months);
W6: July 1990-March 1991 (9 months).
The total duration of these periods W = W1 U W2 U W3 U W4 U W5 U W6 is 73 months.
[77] The original index G10FF and 5 functions RIP(m),
RINVMTQ(m), KLHELL(m) , KLUINC(m) , RFYGM3(m)
were used to formulate the algorithm for prediction of the recession end. Abbreviation the functions are the same as
given in Section 3.3. When function Rf (m) =f m) <3>
is calculated for periods W2 — W6 q is the first month after the end (the through) of the previous recession; when it is calculated for period W1 q is January 1960.
[78] Discretization has been made for these six time series g(m) by means of thresholds Tg(Q) that have been determined considering all months included in the set W. The relevant values of Q are given in Table 6. The values g(m) are coded on the lowest level of resolution, 0 or 1, discriminating only the values g(m) < Tg(Q) and g(m) > Tg (Q). The same notation “1”” has been given to the values, which became more frequent, when the recession end is approaching. Accordingly, “1” is specified for the values g(m) > Tg(Q) in the case of G10FF and LUK10, and for the values g(m) < Tg(Q) in the case of the four other time series.
[79] Hypothetical prediction algorithm. After discretization the monthly description of situation during a recession is reduced to a binary vector with 6 components and the “ideal” situation prior to the end of a recession, when all indicators are precursory, is the vector (1,1,1,1,1,1), which is called the kernel. Let D(m) be the number of zeros in a
Table 6. Precursory trends and values of Q
g(m) Precursory values Q, %
IP Low 75.0
INVR Low 50.0
G10FF High 33.3
LHK5 Low 75.0
LUK10 High 50.0
FYG3R Low 50.0
code of a month - its Hamming distance from the kernel. The approach of the end of a recession may be recognized by the small values of D(m). This suggests the following algorithm: the precursory pattern appears if for three consecutive months D(m) < 3 and the recession end is expected during an interval of three months long after appearing this pattern. The alarms (continuous intervals of this kind) obtained by this algorithm for the recessions considered are shown in Figure 8. The end of each of the 6 recessions is preceded by an alarm and there are no false alarms. The total duration of alarms in all 6 recessions is 16 months, which is 22% of total duration of the set W (73 months).
Figure 8. Results of application of the algorithm for prediction of the end of a recession.
Table 7. FAUs in France, 1965-1997
Time, year: month January 1970 January 1974 September 1977 July 1980 July 1983 May 1990 September 1995
“Magnitude” F* 7.6 22.6 5.3 15.7 9.2 20.3 9.4
Duration (me — m*), 22 24 16 20 19 21 20
months
[80] The end of the last recession. The last recession (#7 in Table 4; Figure 8) starts in April 2001. According to the score D(m) the alarm for the end of that recession starts in July 2001. In 2003 the NBER made a conclusion that it did end in November 2001.
3.5. Fast Acceleration of Unemployment
[81] Each episode of a sharp increase in the unemployment rate is called here Fast Acceleration of Unemployment, FAU. The study [Keilis-Borok et al., 2005] described below used databases issued by the Organization for Economic Cooperation and Development [OECD, 1997] and the International Monetary Fund [IMF, 1997]. Past FAUs were identified by an analysis of the monthly statistics of unemployment. To explore the predictability of FAUs the monthly indicators listed below were analyzed:
[82] 1. IP: Industrial production index, composed of weighted production levels in numerous sectors of the economy, in % relative to the index for 1990.
[83] 2. L: Long-term interest rate on 10-year government bonds, in %.
[84] 3. S: Short-term interest rate on 3-month bills, in %.
[85] For France the data sources are sufficiently complete for the time period between January 1965 and May 1997, and it is this period that is considered here. The analogues of these indicators for the USA have been successfully used in research on predicting American economic recessions [Keilis-Borok et al., 2000; Stock and Watson, 1993; Section 3.3.].
[86] The definition of FAU described above qualitatively (see Figure 4) is formalized here. Let u(m) is the monthly number of unemployed including seasonal variations (m = 1,2,...). FAUs are defined as follows. First, smoothing out the seasonal variation of u(m) a function U (m) =
m
W “I -------, m + 6 ) is obtained - a value of the regression
m — 6
(1) over the time interval (m — 6, m + 6). Next, a function F (—) = KU (m+s, m) — KU(m, m — s) is determined - the
Table 8. Trends and thresholds
Indicator Premonitory trend s Q, %
IP: Industrial production index Upward 12 50%
L: Interest rate, long-term bonds Upward 12 33%
S: Interest rate, short-term bills Upward 12 25%
difference between the linear trends in regression (1) of U(m) within s subsequent months and s preceding months. This function with s = 24 months is used as a coarse measure of unemployment acceleration. Finally, the FAUs are defined by the local maxima of F(m) exceeding a certain threshold F. The time m* and the height F* of such a maximum are, respectively, the time and the magnitude of a FAU. Acceleration ends in a month me of the subsequent local minimum of F (m).
[87] Monthly unemployment in France and the function F(m) for the time period considered, from January 1965 through May 1997, are shown in Figure 9. One may see 10 that the threshold 4 identifies obviously outstanding peaks of F(m). Seven FAUs identified by the condition F* >F=4 are listed in Table 7. As we see, each such FAU is the beginning of a long unemployment rise, lasting 16 to 24 months. Since it is determined after a strong smoothing of the unemployment rate, the meaningful accuracy of prediction may hardly be better than about 2 months. Three “major” FAUs, marked in bold in Table 7, are distinctly larger than the others.
[88] This definition of the FAUs is applicable only in retrospect, two years after a FAU occurs in order to ensure a reliable identification of past FAUs.
[89] The “premonitory” trends of the indicators, which
tend to occur more frequently as a FAU approaches are explored here. The trends are approximated by the function Kf ) (2) with f replaced by the symbol of an indicator.
[90] Discretization has been made for each of these functions g(m) by means of thresholds Tg (Q) that have been determined using all months under consideration.
[91] Premonitory trends. Values of thresholds Tg(Q) and empirically determined premonitory trends of the indicators are summarized in Table 8. The premonitory behavior of the indicators has a transparent qualitative explanation.
[92] The values g(m) are coded on the lowest level of resolution, 0 or 1, discriminating only the values g(m) < Tg (Q) and g(m) > Tg (Q) and the description of the unemploymentrelevant situation is reduced to a monthly time series of a binary vector with 9 components, as is usual in the pattern recognition of infrequent events. For convenience, the same code, 1, is given to the “premonitory” trend of each indicator, regardless of whether it is an upward or a downward one.
[93] It is considered here how the approach of a FAU is reflected in the collective behavior of the indicators. Its simplest description is function D(m) - the number of nonpremonitory indicators for the month m. If the identification of premonitory trends is correct then the value of D(m)
Figure 9. Unemployment in France. Top: Monthly unemployment, thousands of people. Thin line: u(m), data from the OECD database; note the seasonal variations. Thick line: U(m), data smoothed over one year. Bottom: Determination of FAUs. F(m) shows the change in the linear trend of unemployment U(m). FAUs are attributed to the local maxima of F(m) exceeding threshold F = 4.0 shown by solid horizontal line. The thick vertical lines show moments of the FAUs.
should decrease as a FAU approaches. By definition D(m) is the number of zeros in the binary code of the situation. This is the so-called “Hamming distance” between that code and the code of the “perfect” premonitory situation, when all the components are equal to 1, that is, all the trends are all premonitory.
[94] Since only three indicators IP, L, and S are considered the value of D(m) may vary from 0 to 3. Change of D(m) through the time considered is juxtaposed with FAUs in Figure 10. One can see that the minimal value D(m) = 0 appears within 1 to 12 months before a FAU and at no other time. Data on Figure 10 suggest the following hypothetical prediction algorithm: An alarm is declared for 6 months after each month with D(m) = 0 (regardless of whether this month belongs or not to an already determined alarm). A waiting period of 6 months is introduced because in three cases (1977, 1980, and 1995) the premonitory pattern does not appear right before a FAU. Results of prediction are
shown in Figure 11. One can see that this algorithm predicts 6 out of 7 FAUs, including all three major ones.
[95] The algorithm was also applied to the data on monthly unemployment rates for the U.S. civilian labor force, as given by USDL. Unlike Europe, unemployment in USA had no general trend during the years considered. One can see this in Figure 12. The FAUs are the times when unemployment started to rise, that are the local minima of the unemployment rate. They are formally defined as follows. Let R(m) be the smoothed monthly unemployment rate in a month m. Then R(m) has the local minima in a month m* if for j = 1, 2, 3,4 R(m* — j) > R(m*) and R(m* + j) > R(m*). Seven such minima are identified within the period 19601999 in August 1962 (9), March 1967 (3), February 1969
(28), July 1973 (24), May 1979 (19), March 1981 (21), and May 1989 (38). The duration of the unemployment rise is given in brackets after the corresponding months m* , which are the targets of our prediction.
Figure 10. Collective performance of premonitory trends. premonitory trends at month m. Vertical lines show FAUs. declared for 6 months after each month when D(m) = 0.
Function D(m) is the number of nonAlarms (shown by dark gray bars) are
Ltll______
Figure 11. Retrospective prediction of FAUs in France: alarms for FAUs (shown by black bars) and periods of the unemployment growth (shown by grey bars) FAUs. Checkered bars indicate the times, for which data on economic indicators were unavailable.
[96] Application of the algorithm. Indicators IP, S, and L have the following American equivalents (see Section 2.3.). For IP - “industrial production, total”, IP; for S
- interest rate on 90-day U.S. treasury bills, at an annual rate (in percent), FYGM3; for L - interest rate on 10-year U.S. treasury bonds, at an annual rate (in percent); FYGT10.
[97] Alarms and FAUs are juxtaposed in Figure 13. One can see that 4 out of 7 FAUs are captured by alarms; three FAUs, in 1962, 1969, and 1981, are missed; and there
are three false alarms, in 1968, 1983, and 1994. The alarms within the periods of unemployment growth are not regarded as false ones. Total count of errors for the USA is worse that for France, though the result is better than random.
[98] First advance prediction. An advance prediction by the suggested algorithm is shown in Figure 13. It was found by analysing the data for USA up to December 2000 that D(m) = 0 during the four months, from February to May 2000. Accordingly, the algorithm declared the alarm for
12
2----------------------------------:------ ------------- I I-----------—-------: ■ -:-------i.
1960 1965 1970 1975 1980 1985 1990 1995
Figure 12. Unemployment rates in the U.S. Thin line: r(m), original data. Thick line: R(m), data after smoothing out the seasonal variations. The thick vertical lines show the moments when unemployment started to rise (local minima of smoothed unemployment rate).
Figure 13. Prediction of FAUs in the USA. Notations are the same as in Figure 11.
the period from February to November 2000 [Keilis-Borok et al., 2001]. The data on unemployment rates for the subsequent period, up to July 2002, confirmed that prediction: a FAU materialized in July 2000. Obviously, that confirmation is most encouraging, but much longer experimentation is necessary to validate the algorithm.
[99] Numerous other warnings of a coming rise of unemployment in the USA did appear during the first months of 2000, and even in the popular media. The particular feature of the prediction discussed here, however, is that a formal unambiguous algorithm obtained it and that it indicates a specific time interval when the unemployment will start to rise.
[100] The second alarm has been declared from May 2006 to April 2007 (Figure 14). Now (in December 2007) it is impossible to determine confidently that there is FAU during this period because the last month for that the data on the unemployment rate are available is October 2007. When the data for November 2007-January 2008 will be available we shall be able to decide is there FAU in March 2007 or not.
3.6. Homicide Statistics
[101] Dynamics of crimes reflects important aspects of sustainability of our society and the risk of its destabilization -a prelude to a disaster. Here, a prominent feature of crime dynamics - surge of the homicides in a mega-city is considered [Keilis-Borok et al., 2003]. The study integrates the professional expertise of the police officers and of the scientists working on pattern recognition of infrequent events.
[102] In this study statistics of several types of crime in Los Angeles over the period 1975-2002 is analysed. The analysis focuses on how these statistics change before a sharp and lasting rise (“a surge”) of the homicide rate. The goal is to find an algorithm for predicting such a surge by monitoring the rates of different crimes.
[103] The following data sources are used:
[104] (i) The National Archive of Criminal Justice Data, placed on the web site (NACJD: http://www.icpsr.umich. edu/NACJD/index.html). Carlson [1998] gives its description. This site contains data for the years 1975-1993.
Figure 14. Advance prediction of FAUs in the USA. Unemployment rate in the USA, 1997-2002: thin curve shows original data; thick curve shows the rates with seasonal variations smoothed out. The gray bars show the alarm periods, defined by analysis of macroeconomic indicators. The black vertical line shows the actual start of the unemployment rise, as defined by an analysis of monthly unemployment rates for the U.S. civilian labor force.
Figure 15. Total monthly number of homicides in Los Angeles city, 1975-1993. Data are taken from (NACJD: http://www.icpsr.umich.edu/NACJD/index.html), Carlson [1998]. Thin curve - original time series, h(m), per 3,000,000 inhabitants. Thick curve - smoothed series H(m), with seasonal variations eliminated. Vertical lines show the targets of prediction - episodes of SHS. Gray bars are the periods of homicide surge. Checkered bars are the alarms declared by the hypothetical prediction algorithm.
[105] (ii) Data bank of the Los Angeles Police Department (LAPD Information Technology Division); it contains similar data for the years 1990 - May 2001.
[106] Here the data for 1975-1993 as taken from (NACJD: http://www.icpsr.umich.edu/NACJD/index.html), Carlson [1998] are analyzed. Let h(m), m =1, 2 ... , be the time series of the monthly number of all homicides. Figure 15 shows the plot of h(m) in Los Angeles, per 3,000,000 inhabitants of the city. To identify the episodes of SHS (see Figure 5 above) the seasonal variations, which are clearly seen in Figure 15, are smoothed out by replacing h(m) with its linear
least square regression (1): H(m) = Wh ( —m—, m + 6).
\m — 6 /
Since H(m) is defined on the time interval (m — 6, m + 6),
it depends on the future. Thus, it is admissible to define prediction targets (but not precursors).
[107] The function H(m) is shown in Figure 15 by the thick curve. Three time periods of a lasting homicide rise are clearly seen: 1977-1980, 1988-1992 and a relatively shorter period 1985-1986. The starting months of these periods: April 1977, March 1985, and August 1988 are chosen as prediction targets. They are marked in Figure 15 by vertical lines.
[108] Here the monthly data on seven types of crimes out of the 13 types listed in Table 9 are analyzed to look for “premonitory” trends of each crime that tend to appear more frequently as an SHS approaches. Prediction itself is based
Table 9. Types of crime considered (after Carlson [1998]; abbreviations are indicated in brackets)
Homicide Robberies Assaults Burglaries
• All (H) • All (Rob) • All (A)* • Unlawful not forcible entry (UNFE)
• With firearms (FRob) • With firearms (FA)
• With knife or cutting • With knife or cutting • Attempted forcible entry (AFE)*
instrument (KCIR) instrument (KCIA)
• With other dangerous • With other dangerous
weapon (ODWR) weapon (ODWA)*
• Strong-arm robberies • Aggravated injury
(SAR)* assaults (AIA)*
Note: * Analyzed in sensitivity tests only.
on the collective behavior of these trends, as analyzed below. Orientation on a set of precursors has been found to be rather successful in prediction research: an ensemble of “imprecise” precursors usually gives better predictions than a single “precise” precursor [e.g. Keilis-Borok and Rotwain, 1990; Zaliapin et al., 2003].
[109] Observation. According to police experience, the crimes considered here often rise before an SHS.
[110] To smooth out seasonal variations, we replace the plot f (m) of each type of crime is replaced by its regression
(1): F(m) = Wf (—^ —, m ). Regression is done over
m — 12
the prior 12 months and does not depend on the future, so that it can be used for prediction. These plots exhibit two consecutive patterns:
[111] (i) First, one can see a simultaneous escalation of burglaries and assaults within several (4 to 11) months before an SHS; at the same time robberies are declining.
[112] (ii) Later on, closer to an SHS, one can see, albeit not so clearly, a simultaneous escalation of different kinds of robberies.
[113] The first pattern is formally defined and explored here. To quantify the above observation the function Kf
(2) is used to approximate the trends of the crimes where f identifies the type of crime. Next, following the pattern recognition approach, the trends (the values of Kf ) are discretized on the lowest level of resolution: a binary one distinguishes only the trends above and below a threshold Tg (Q) where g denotes a relevant function.
[114] According to the above observations, it is expected that “premonitory” trends lay above the respective thresholds for assaults and burglaries, while they lay below these thresholds for robberies. One can see this in Figure 16, showing the functions Kf (m — 12, m) for 7 crime types. For convenience, the same code, 1, is given to the “premonitory” trend of each crime, regardless of whether it is above or below the threshold of discretization. The seven monthly crime statistics considered here are thus reduced to a binary vector with 7 components.
[115] The values of Q used for discretization are given in Table 10.
[116] The simplest description of the collective behavior of the trends is D(m) - the number of non-premonitory trends at a given month m. If the identification of premonitory trends is correct then D(m) should be low in the proximity of an SHS. By definition D(m) is the number of zeros in the binary code of the monthly situation. This is the so-called “Hamming distance” between that code and the code of the “pure” premonitory situation, {1,1,1,1,1,1,1} when all seven trends listed in Table 10 are premonitory.
[117] Figure 17 shows the change of D(m) with time. The value of D(m) may vary from 0 to 7 but the minimal observed value is 1. That value appears within 4 to 11 months before an SHS and at no other time. An examination of the temporal change of D(m): An alarm is declared for 9 months each time when D(m) < A for two consecutive months (regardless of whether these two months belong or not to an already declared alarm).
[118] The condition D(m) < A means, by definition, that A or less trends are not premonitory at the month m. A count of D(m) (Figure 17) suggests to take A = 1. A waiting period of 9 months is introduced because the premonitory trends do not appear right before an SHS. The requirement that this condition holds for two months in a row makes prediction more reliable and reduces the total duration of alarms.
[119] The alarms obtained by this algorithm are shown in Figure 17 by the checkered bars. The total duration of these alarms is 30 months, representing 14 percent of all months considered. In real prediction such a score would be considered quite satisfactory.
[120] The algorithm has been tested by application to “out of sample” data not used in its development. Such tests are always necessary to validate and/or improve a prediction algorithm. Such a test is possible since the algorithm is self-adaptive: the thresholds Tg (Q) are not fixed but are adapted to crime statistics, as the percentile of a level Q.
[121] Los Angeles, 1994—2002. So far we used the data source [Carlson, 1998], (NACJD: http://www.icpsr.umkh. edu/NACJD/index.html) covering the years 1975-1993 was used. To extend the analysis past 1993, the data of the LAPD Information Technology Division, covering the time period from January 1990 to May 2002 have been involved. It has been found by comparing the data for the overlapping three years that they are reasonably close, particularly after smoothing.
[122] Figure 18 shows the homicide rates through the whole period from 1975 to May 2002. Two SHS episodes are identified in the later period 1994-2001. They are indicated in Figure 16 by dashed vertical lines. The first episode is captured by an alarm, which starts in the month of SHS without a lead time. The second episode is missed in that an alarm has started two months after it. That error has to be put on the record; nevertheless the prediction remains informative: during these two months the monthly homicide number rose by only a few percent, giving no indication that a lasting homicide surge has started.
[123] New York City. Figure 19 shows the monthly total homicide rates in New York City per 7 million inhabitants of the city. Two SHS episodes (February 1978 and February 1985) are identified. The prediction algorithm gives two alarms, as shown in Figure 19 by checkered bars. One of them predicts the second SHS, while the first one is
Table 10. Premonitory trends for selected crime types
# Crime type Premonitory trend s Q, %
1 Rob Below threshold 12 66.7
2 FRob ii 11 12 66.7
3 KCIR ii 11 12 50.0
4 ODWR ii 11 12 87.5
5 FA Above threshold 12 50.0
6 KCIA ii 11 12 50.0
7 UNFE ii 11 12 50.0
Figure 16. The regression coefficients Kf(m — 12, m) for seven crime types. See the definition (2) in Section 3.2. and notations in Table 9. Original data are taken from (NACJD: http://www.icpsr.umich.edu/NACJD/index.html), Carlson [1998]. Horizontal lines and arrows show respectively discretization thresholds and premonitory trends in accordance with Table 10. Vertical lines show episodes of SHS. Gray bars indicate months when D(m) < 1.
Figure 17. Homicide surges and alarms determined by the prediction algorithm. Start of a homicide surge is shown by the vertical line. Function D(m) is the number of crime statistics not showing premonitory trends at a month m. Alarms (shown by checkered bars) are declared for 9 months, when D(m) < 1 during two consecutive month. Adjustable parameters correspond to variant 10 given in Table 7.
missed. Another alarm is considered as a false one. Though the failure to predict and a false alarm are disappointing, the results as a whole appear to be useful: one of the two SHS is captured by alarms lasting together 21 months, amounting to 10 percent of the time interval considered.
3.7. Discussion of Socio-Economical Implications
[124] If the crime dynamics is considered then on the practical side, the results given above enhance the capability to identify a situation that is “ripe” for homicide surges and, accordingly, to escalate the crime prevention measures. In a broader scheme of things, a surge of crime is one of potential ripple effects of natural disasters. Accordingly the risk of a natural disaster is higher in such a situation. The approach applied here - a heuristic “technical” analysis - is not competing with but complementary to the cause-and-effect “fundamental” analysis. The cause that triggered a specific homicide surge is usually known, at least in retrospect. This might be, for example, a rise in drug use, a rise in unemployment, a natural disaster etc. However, that does not render predictions considered in this study redundant. On the contrary, this approach might predict an unstable situation when a homicide surge might be triggered, thus enhancing the reliability of cause-and-effect predictions. Among available data that can be incorporated in the analysis are other types of crimes [Bursik et al., 1990], economic and demographic indicators [Messner, 1983] and the territorial distribution of crimes. It seems worthwhile to try the same approach with other targets of prediction - e.g. surges of all violent crimes; and to other areas, e.g. separate Bureaus of the city of Los Angeles, or to other major cities. In a broader scheme of things, this analysis discriminates
Figure 18. Performance of prediction algorithm through 1975-2002. Data from [Carlson, 1998], (NACJD: http://www.icpsr.umich.edu/NACJD/index.html) for 1975-1993 have been used to develop the algorithm. It was than applied to the data from the Data Bank of the Los Angeles Police Department (LAPD Information Technology Divison) for subsequent 9 years. Notations are the same as in Figure 15. Dashed vertical lines indicate SHS episodes that occurred after 1993.
Figure
tions a
19.
the
Application of the prediction Figure 15. Data
same as in
http://www.icpsr.umich.edu/NACJD/index.html).
inhabitants.
algorithm to New York City. Nota-
are taken from [Carlson, 1998], (NACJD: Homicide statistics is shown per 7,000,000 of
stable situations from unstable, where the risk of different disasters is higher. At the same time it would be important to set up an experiment in advance prediction of homicide surges in Los Angeles using the algorithm hypothesized here. Successes and errors will both provide for evaluation of this algorithm and for developing a better one.
[125] The alarms, obtained here, are rather durable. In a second approximation they may be used as a background for the search of the shorter-term ones. Note, that the shorter alarms are not necessarily more useful for practical purposes [e.g. Molchan, 1994].
[126] On theoretical side if confirmed by advance prediction these findings would expand the known limits of predictability of socio-economic systems. They and the realm of robust behavior patterns in macroeconomics provide heuristic constraints for macroeconomic models. Besides econometric models, as reviewed by Engle and McFadden [1994], this refers also to models of non-linear systems such as developed in statistical mechanics [Allegre et al., 1995; Blanter et al., 1997; Burridge and Knopoff, 1967; Gabrielov et al., 2000; Newman et al., 1994; Zaliapin et al., 2003].
[127] The economic integration of the European Union is increasing, along with the general globalization of the economy. This will not necessarily render the analysis irrelevant, since the algorithms are very robust. In any case, economic integration would make it easier to develop new algorithms
of this kind, with an even higher level of averaging of the processes considered.
[128] Similarly, the results may not become irrelevant due to some drastic change of the mechanisms controlling the quickly accelerating modern economy. The premonitory patterns considered here probably reflect some type of scenario of transition to critical phenomena (one that is common for many mechanisms in the case of non-linear dynamics). An indirect confirmation is the uniform performance of the prediction algorithm over the last 35 years.
[129] The results described above suggest the following directions for further research: (i) here the indicators that are sufficient for prediction have been identified, but not all indicators potentially useful for that purpose; (ii) with less robust predictions, premonitory patterns may be considered in the context of the accelerator hypothesis and, more generally, the cointegration of each indicator separately and in combinations, using Dickey-Fuller tests and Granger causality concepts [Ericsson, 1998; Ericsson and MacKinnon, 1999; Watson, 1994]; (iii) in the search for subsequent approximations to prediction it would be interesting to explore, e.g., the changes of unemployment during the FAUs; also, it may be preferable for purely computational reasons to consider the reciprocal of unemployment as used in previous studies and the prediction of a decrease of unemployment obviously deserves a parallel study.
4. Final Remarks
[130] The examples of the on-going prediction of extreme events in geophysical and socio-economical systems demonstrate the efficiency, confirmed and potential, of the algorithms designed by a common pattern recognition approach for solving itchy challenging problems of practical interest. The necessary prerequisites of verification of a pattern recognition solution are the on-line availability of data used in computation for relevant integrals and indicators, as well as a setup of an unambiguous test experiment. The first is feasible in many branches of Science now, while the second needs heuristic search of potential solutions and special efforts of durable monitoring of the phenomenon in question.
[131] The achieved experience in the straight forward practical approach to earthquake prediction problem provided already a unique collection of successes and failures that permit their systematic analysis and further development of the methodology. Obviously, the progress in Quantitative Earthquake Prediction will require more data, novel pioneering studies, and verification of arising hypotheses on correlations between the occurrence of seismic events and observable phenomena.
[132] Similar to earthquake prediction studies decisive validation of the socio-economic findings described above requires experimentation in advance prediction. Particularly encouraging for further research is the wealth of yet untapped possibilities, since only a small part of the data and mathematical models that are currently available in the fields and that are evidently relevant to considered critical transitions have been used so far.
References
Allegre, C. J., J.-L. Le Mouel, C. Ha Duyen, and C. Narteau (1995), Scaling organization of fracture tectonics (SOFT) and earthquake mechanism, Phys Earth Planet. Inter, 92, 215, doi:10.1016/0031-9201(95)03033-0.
Allen, C. R., (Chaiman), W. Edwards, W. J. Hall, L. Knopoff, C. B. Raleigh, C. H. Savit, M. N. Toksoz, and R. H. Turner (1976), Predicting earthquakes: A scientific and technical evaluation, with implications for society, Panel on Earthquake Prediction of the Committee on Seismology, Assembly of Mathematical and Physical Sciences, National Research Council, p. 62, U.S. National Academy of Sciences, Washington, D.C..
Bhatia, S. C., S. V. Chalam, V. K. Gaur, V. I. Keilis-Borok, and V. G. Kossobokov (1989), On intermediate-term presic-tion of strong earthquakes in the Himalayan arc region using pattern recognition algorithm M8, Proc. Indian Acad. Sci., Earth Planet. Sci., 98(1), 111.
Blanter, E. M., M. Shnirman, J.-L. Le Mouel, and C. Allegre (1997), Scaling laws in blocks dynamics and dynamic selforganized criticality, Phys. Earth Planet. Inter., 99, 295, doi:10.1016/S0031-9201(96)03195-0.
Burridge, R., and L. Knopoff (1967), Model and theoretical seismicity, Bull. Seismol. Soc. Am., 57, 341.
Bursik, R. J., Jr., H. G. Grasmick, and M. B. Chamlin (1990), The effect of longitudional arrest patterns on the development of robbery trends at the neighborhood level, Criminology, 28, 431, doi:10.1111/j.1745-9125.1990.tb01333.x.
Carlson, S. M. (1998), Uniform Crime Reports: Monthly
Weapon-specific Crime and Arrest Time Series, 1975—1993 (Na-
tional, State, and 12-City Data), ICPSR 6792, Inter-university Consortium for Political and Social Research, p. P.O. Box 1248, Ann Arbor, Michigan 48106.
Crutchfield, J. P., J. D. Farmer, N. H. Packard, and R. S. Shaw (1986), Chaos, Sci. Am., 255, 46.
Cyranoski, D. (2004), A seismic shift in thinking, Nature, 431, 1032, doi:10.1038/4311032a.
Ellsworth, W. L., and G. C. Beroza (1995), Seismic evidence for an earthquake nucleation phase, Science, 268, 851, doi:10.1126/science.268.5212.851.
Engle, R. F., and D. L. McFadden, Eds., (1994), Handbook of Econometrics, 4, 1078 pp., North-Holland, Amsterdam.
Ericsson, N. R. (1998), Empirical modeling of money demand, Empirical Economics, 23, 295.
Ericsson, N. R., and J. G. MacKinnon (1999), Distributions of error correction tests for cointegration, International Finance Discussion, p. Papers 655, Board of Governors of the Federal Reserve System, http://ideas.repec.org/a/ect/emjrnl/ v5y2002i2p285-318.html.
Gabrielov, A., V. Keilis-Borok, I. Zaliapin, and W. I. Newman (2000), Critical transitions in colliding cascades, Phys. Rev.
E, 62, 237, doi:10.1103/PhysRevE.62.237.
Gabrielov, A., W. I. Newman, and D. L. Turcotte (1999), An exactly soluble hierarchical clustering model: Inverse cascades, self-similarity, and scaling, Phys. Rev. E, 60, 5293, doi:10.1103/ PhysRevE.60.5293.
Gahalaut, V. K., I. Kuznetsov, V. Kossobokov, A. Gabrielov, and V. Keilis-Borok (1992), Application of pattern recognition algorithm in the seismic belts of the Indian convergent plate margins — M8 algorithm, Proc. Indian Acad. Sci. (Earth Planet. Sci.), 101(3), 239.
Gelfand, I. M., Sh. A. Guberman, V. Keilis-Borok, L. Knopoff,
F. Press, I. Ranzman, I. M. Rotwain, and A. M. Sadovsky (1976), Pattern recognition applied to earthquake epicenters in California, Phys. Earth Planet. Inter., 11, 227, doi:10.1016/0031-9201(76)90067-4.
Geller, R. J., D. D. Jackson, Y. Y. Kagan, and F. Mulargia
(1997), Earthquakes cannot be predicted, Science, 275, 1616, doi:10.1126/science.275.5306.1616.
Gell-Mann, M. (1994), The Quark and the Jaguar: Adventures in the Simple and the Complex, 392 pp., Freeman and Company, New York.
Gorshkov, A., V. Kossobokov, and A. Soloviev (2003), 6.
Recognition of earthquake-prone areas, in Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, edited by V. I. Keilis-Borok and A. A. Soloviev, p. 141, Springer, Heidelberg.
Gvishiani, A. D., and V. G. Kosobokov (1981), On foundations of the pattern recognition results applied to earthquake-prone areas, Izv. Phys. Solid Earth (in Russian), 2, 21.
Harte, D., D.-F. Li, M. Vreede, and D. Vere-Jones (2003),
Quantifying the M8 prediction algorithm: Reduction to a single critical variable and stability results, New Zealand J. Geol. Geophys., 46, 141.
Healy, J. H., V. G. Kossobokov, and J. W. Dewey (1992), A test to evaluate the earthquake prediction algorithm, M8, Open File Rep. 92-401 (with 6 Appendices), p. 23, Geol. Surv. of U.S., Reston, VA.
International Monetary Fund, (IMF) (1997), International
Financial Statistics, CD-ROM pp., World Bank, Washington, DC.
Keilis-Borok, V. I. (1990), The lithosphere of the Earth as a nonlinear system with implications for earthquake prediction, Rev. Geophys., 28, 19, doi:10.1029/RG028i001p00019.
Keilis-Borok, V. I., and V. G. Kossobokov (1984), A complex of long-term precursors for the strongest earthquakes ofthe world, Proc. 27th Geological Congress, 61, p. 56, Nauka, Moscow.
Keilis-Borok, V. I., and V. G. Kossobokov (1987), Periods of high probability of occurrence of the world’s strongest earthquakes, Computational Seismology, Allerton Press Inc., 19, 45.
Keilis-Borok, V. I., and V. G. Kossobokov (1990a), Premonitory activation of earthquake flow: Algorithm M8, Phys. Earth Planet. Int., 61, 73, doi:10.1016/0031-9201(90)90096-G.
Keilis-Borok, V. I., and V. G. Kossobokov (1990b), Times of increased probability of strong earthquakes (M > 7.5), Diagnosed by algorithm MS in Japan and adjacent territories, J. Geophys. Res., 95(BS), 12,413.
Keilis-Borok, V. I., and A. J. Lichtman (1993), The self-
organization of American society in presidential and senatorial elections, in Limits of Predictability, edited by Yu. A. Kravtsov, p. 223, Springer-Verlag, Berlin-Heidelberg.
Keilis-Borok, V. I., and F. Press (19S0), On seismological applications of pattern recognition, in Source Mechanism and Earthquake Prediction Applications, edited by C. J. Allegre, p. 51, Editions du Centre national de la recherche scientifique, Paris.
Keilis-Borok, V. I., and I. М. Rotwain (1990), Diagnosis of
time of increased probability of strong earthquakes in different regions of the world: algorithm CN, Phys. Earth Planet. Int., 61, 57, doi:10.1016/0031-9201(90)90095-F.
Keilis-Borok, V. I., and A. A. Soloviev, Eds., (2003), Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, 337 pp., Springer-Verlag, Berlin-Heidelberg.
Keilis-Borok, V. I., L. Knopoff, I. М. Rotwain, and C. R. Allen (19SS), Intermediate-term prediction of occurrence times of strong earthquakes, Nature, 335(6192), 690, doi:10.103S/ 335690a0.
Keilis-Borok, V., J. H. Stock, A. Soloviev, and P. Mikhalev (2000), Pre-recession pattern of six economic indicators in the USA, J. Forecast, 19, 65, doi:l0.l002/(SICI)l099-l3lX(20000l) l9:l<65::AID-FOR73O>3.O.CO;2-U.
Keilis-Borok, V. I., A. Soloviev, C. Allegre, A. Sobolevskii, and М. Intriligator (2001), Dynamics of Macroeconomic Indicators before the Rize of Unemployment in Western Europe and the USA, Sixth Workshop on Non-Linear Dynamics and Earthquake Prediction, 15—27 October 2001, H4.SMR/1330-11, p. 36, ICTP, Trieste.
Keilis-Borok, V. I., D. Gascon, A. A. Soloviev, М. Intriligator, R. Pichardo, and F. E. Winberg (2003), On predictability of homicide surges in megacities, in Risk Science and Sustainability, edited by T. Beer, A. Ismail-Zadeh (NATO Science Series. II. Mathematics, Physics and Chemistry, vol. 112), p. 91, Kluwer Academic Publishers, Dordrecht-Boston-London.
Keilis-Borok, V. I., A. Soloviev, C. Allegre, A. Sobolevskii, and М. D. Intriligator (2005), Patterns of macroeconomic indicators preceding the unemployment rise in Western Europe and the USA, Pattern Recognition, 38(3), 423, doi:l0.l0l6/j.patcog.
2004.0S.005.
Klein, Ph. A., and М. P. Niemira (1994), Forecasting Financial and Economic Cycles, 542 pp., Wiley, New York.
Kossobokov, V., and P. Shebalin (2003), 4. Earthquake Prediction, in Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, edited by V. I. Keilis-Borok, and A. A. Soloviev, p. 141, Springer, Heidelberg.
Kossobokov, V. G., and A. V. Khokhlov (1993), An experimental intermediate-term prediction of earthquakes in real time: Verification of the MS algorithm, in Mathematical Modeling of Seismotectonic Processes in the Lithosphere Oriented on the Problem of Earthquake Prediction, Issue 1 (in Russian), p. 53, MITP RAN, Moscow.
Kossobokov, V. G., B. K. Rastogi, and V. K. Gaur (19S9), On self-similarity of premonitory patterns in the regions of natural and induced seismicity, Proc. Indian Ac. Sci., Earth and Planetary Sciences, 98(4), 309.
Kossobokov, V. G., J. H. Healy, V. I. Keilis-Borok, J. Dewey, and A. V. Khokhlov (1992), Testing an intermediate-term earthquake prediction algorithm: A real-time test design and the results of retrospection, Doklady RAN (in Russian), 325(1), 46.
Kossobokov, V. G., L. Romashkova, G. Panza, and A. Peresan (2002), Stabilizing intermediate-term medium-range earthquake predictions, J. Seismology Earthquake Engineering, 4(2— 3), 11.
Kossobokov, V. G., V. I. Keilis-Borok, and S. W. Smith (1990), Localization of intermediate-term earthquake prediction,, J. Geophys. Res., 95(B12), 19,763.
Kossobokov, V. G., V. I. Keilis-Borok, L. L. Romashkova, and
J. H. Healy (1999), Testing earthquake prediction algorithms: Statistically significant real-time prediction of the largest earthquakes in the Circum-Pacific, 1992—1997, Phys. Earth Planet. Int., 111(3-4), 187.
Kossobokov, V. G. (1997), User Manual for M8, in Algorithms for earthquake statistics and prediction, edited by J. H. Healy, V. I. Keilis-Borok, and W. H. K. Lee, IASPEI Software Library, vol. 6, Seismol. Soc. Am., p. 167, El Cerrito, CA.
Latoussakis, J., and V. G. Kossobokov (1990), Intermediate term earthquake prediction in the area of Greece: Application of the Algorithm M8, Pure Appl. Geophys., 134(2), 261, doi:10.1007/BF00877001.
Lichtman, A. J., and V. I. Keilis-Borok (1989), Aggregate-level analysis and prediction of midterm senatorial elections in the United States, 1974-1986, Proc. Natl. Acad. Sci. USA, 86(24), 10,176, doi:10.1073/pnas.86.24.10176.
Messner, S. F. (1983), Regional differences in the economic correlates of the urban homicide rate, Criminology, 21, 477, doi:10.1111/j.1745-9125.1983.tb00275.x.
Molchan, G. M. (1994), Models for optimization of earthquake prediction, in Computational Seismology and Geodynamics, edited by D. K. Chowdhury, vol. 1, p. 1, Am Geophys Un, Washington, DC.
Molchan, G. M. (2003), 5. Earthquake prediction strategies: A
theoretical analysis, in Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, edited by V. I. Keilis-Borok, and A. A. Soloviev, p. 141, Springer, Heidelberg.
Mostaghimi, M., and F. Rezayat (1996), Probability forecast of a downturn in U. S. economy using classical statistical theory, Empirical Economics, 21, 255, doi:10.1007/BF01175973.
Newman, W., A. Gabrielov, and D. L. Turcotte, (Eds.), (1994), Nonlinear Dynamics and Predictability of Geophysical Phenomena, http://prola.aps.org/toc/PRE/v61/i4 pp., Am Geo-phys Un., Int Un Geodesy Geophys.
OECD, (1997), Main Economic Indicators: Historical Statistics 1960-1996, CD-ROM pp., OECD, Paris.
Peresan, A., V. Kossobokov, L. Romashkova, and G. F. Panza (2005), Intermediate-term middle-range earthquake predictions in Italy: A review, Earth-Sience Reviews, 69(1-2), 97, doi:10.1016/j.earscirev.2004.07.005.
Press, F., and C. Allen (1995), Patterns of seismic release in the southern California region, J. Geophys. Res., 100(B4), 6421, doi:10.1029/95JB00316.
Press, F., and P. Briggs (1975), Chandler wobble, earthquakes, rotation and geomagnetic changes, Nature (London), 256, 270, doi:10.1038/256270a0.
Press, F., and P. Briggs (1977), Pattern recognition applied to uranium prospecting, Nature (London), 268, 125, doi:10.1038/268125a0.
Romachkova, L. L., V. Kossobokov, G. Panza, and G. Costa
(1998), Intermediate-term prediction of earthquakes in Italy: Algorithm M8, Pure Appl. Geophys., 152, 37, doi:10.1007/s000 240050140.
Shebalin, P., V. Keilis-Borok, I. Zaliapin, S. Uyeda, T. Nagao, and N. Tsybin (2003), Short-term premonitory rise of the earthquake correlation range: 18 case histories for the largest earthquakes in Japan (M > 7) and 3 in California (M > 7.4),
in IUGG2003, June 30-July 11, 2003, p. Abstracts, A.184, Sapporo, Japan.
Stock, J. H., and M. W. Watson (1989), New indexes of leading and coincident economic indicators, in NBER Macroeconomics Annual, p. 351, University of Chicago Press, Chicago.
Stock, J. H., and M. W. Watson (1993), A procedure for predicting recessions with leading indicators, in Business Cycles, Indicators, and Forecasting (NBER Studies in Business Cycles, Vol.28), edited by J. H. Stock and M. W. Watson, p. 95, University of Chicago Press, Chicago.
Tukey, J. W. (1977), Exploratory Data Analysis, Addison-Wesley Series in Behavioral Science: Quantitative Methods, 688 pp., Addison-Wesley, Reading (Mass).
Vorobieva, I. A. (1999), Prediction of a subsequent large earthquake, Phys. Earth Planet. Inter., 111, 197, doi:10.1016/S0031 -9201(98)00160-5.
Watson, M. W. (1994), Vector autoregressions and cointegration, Chapter 47, in Handbook of Econometrics, vol. IV, edited by R. F. Engle and D. L. McFadden, p. 2843, North-Holland, Amsterdam.
Wyss, M. (1997), Cannot earthquakes be predicted?, Science, 278, 487, doi:10.1126/science.278.5337.487.
Zaliapin, I., V. Keilis-Borok, and M. Ghil(2003), A Boolean delay model of colliding cascades. II: Prediction of critical transitions, J. Stat. Phys., 111(3-4), 839, doi:10.1023/A:1022802432 590.
V. G. Kossobokov, International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, 117997, Moscow-997, Profsoyuznaya street, S4/32, Russia; Institute de Physique du Globe de Paris, 4 Place Jussieu, 75252 Paris, Cedex 05, France (volodya@mitp.ru)
A. A. Soloviev, International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, 117997, Moscow-997, Profsoyuznaya street, S4/32, Russia (soloviev@mitp.ru)
