Predictability of geomagnetic series

The aim of this paper is to lead a practical, rational and rigorous approach concerning what can be done, based on the knowledge of magnetic series, in the field of prediction of the extreme geomagnetic events. We compare the magnetic vector differential at different locations computed with different resolutions, from an entire day to minutes. We study the classical correlations and the simplest possible prediction scheme to conclude a high level of predictability of the magnetic vector variation. The results obtained are far from a random guessing: the error diagrams are either comparable with earthquake prediction studies or outperform them when the minute sampling is used in accounting for hourly magnetic vector variation. We demonstrate how the magnetic extreme events can be predicted from the hourly value of the magnetic variation with a lead time of several hours. We compute the 2-D empirical distribution of consecutive values of the magnetic vector variation for the estimation of conditional probabilities of different types. The achieved results encourage further development of the approach to prediction of the extreme geomagnetic events.


Introduction
Many have tried for a long time to predict the so-called magnetic situation or magnetic activity, as characterized, for example, by magnetic indices.The objective of better understanding of the time evolution of the geomagnetic field is also of practical interest: the simplest example of application is the planning of an aeromagnetic survey, which requires a quiet magnetic situation to be accurate.
The interest for such type of prediction has, of course, been renewed and amplified since it has been realized that damages caused by big magnetic storms in power lines and Correspondence to: E. Bellanger (ebellan@ipgp.jussieu.fr)other installations (e.g.Kappenman, 1996), have to be mitigated.Space navigation, especially manned flights, and the revolution of communications via satellites also require forecasting of space weather (e.g.Hastings, 1995).
The full subject encompasses short-term, middle-term, long-term predictions, periodic or quasi-periodic cycles.Furthermore space weather forecasting uses additional space information independent of the magnetic field.The objective of the present paper may then appear limited.We will evaluate the predictability of some long and homogeneous magnetic time series by using simple tools whose efficiency can be easily tested.Our results, however limited, are established in a rigorous way and, for each prediction, we will clearly and explicitly say what it means, avoiding any vague statement.

The magnetic time series
A magnetic observatory provides recordings of three components of the geomagnetic field: X, horizontal northward, Y , horizontal eastward, and Z, vertical downward.The sampling rate and the accuracy, in absolute value, depend on the epoch and the observatory.We won't describe here in any detail the full set of these magnetic data, but rather refer to Bellanger et al. (2002b) and Bellanger et al. (2002a).We will just give the necessary information on the series analyzed in the present study.The quality of the series used has been carefully checked: series containing gaps or obvious steps due to base lines problems were rejected.This selection led to the retention of series that do not have the same length and do not cover the same timespan.The long-term control of absolute values (the so-called problem of base lines), which remains the most difficult task to achieve in an observatory, is, contrary to the steps mentioned above, of small influence here, due to the short-term differences considered (daily, at most).

Minute values
Minute values of X, Y and Z have been available in a number of observatories since the introduction of fluxgate variometers (at a date depending on the observatory, after 1970).In particular, the Intermagnet program has produced each year since 1992 a CD-ROM containing definitive minute values of several tens of observatories (80 in 2001); the minute value, given in tenths of nT, is obtained by applying a Gaussian filter of 19 coefficients to a set of 19 measurements centered on the given minute and sampled every 5 s (Trigg, and Coles, 1999).
We will use here the X, Y , and Z components of an 11year minute value time series from the Port-aux-Franc ¸ais observatory and 4 years of minute values from Chambon-la-Forêt (see Table 1).

Hourly means
Hourly mean values of X, Y , Z can be obtained from the World Data Center (WDC) of Copenhagen or directly from the observatories.More than twenty 3-component series are available, which have variable lengths and continuity (see, e.g.Bellanger et al., 2002b).Hourly means, centered on the half hour, are computed from minute data when and where available, by a simple arithmetic averaging.Before minute values were available, hourly means were scaled by hand from photographic magnetograms.We will use here the hourly means of the Chambon-la-Forêt observatory covering the period 1960-1995.

Daily means
From the hourly mean series of the Eskdalemuir observatory (Table 1), covering the 1914-1998 timespan, we build a series of daily means (daily values are obtained by a simple arithmetic averaging of the 24-hourly values of the day), from 1 January 1914 to 31 December 1998 (31 046 days; the ESK series is the longest and most continuous available).

The analysis of predictability
We will analyze the predictability of a magnetic series corresponding to the three sampling rates: daily, hourly and minute.We start the analysis with daily values to conclude with the most refined, minute sampling.

Daily first differences
Let E be any of X, Y and Z; k is the sequential number of the sample (k = 1, 2 . . .N), and the first difference and the daily rate of change R, as well as any function of R, is an invariant measure of geomagnetic field (does not depend on the choice of local coordinates X, Y and Z).The common line of analysis will be described in detail for the ESK series of daily means; the subsequent analysis of hourly and minute data being performed along exactly the same lines will be presented more succinctly.Note that the daily sampled series Ė and R obtained from esk daily means are measured in nT per day.After presenting briefly the histograms of Ė and R, we will analyze the predictability of the R series.

Histograms
As the magnetic data suggest, it is natural to use bins with exponentially increasing sizes to describe their distribution.Specifically, the bin boundaries we use to plot the histograms are x p+1 = ax p = a p x 0 , a > 1, with p being an integer number, so that in log-log scale they are equally spaced along the abscissa.Figure 1 riod 1914-1998, in bi-logarithmic scale (a = 2).The first three ones show a linear increase in the density up to values of 3-4 nT/day, over three orders of magnitude, and then a power law decay, with an exponent between 2 and 3, from 10 nT/day until the "extreme" events with values larger than 100 nT/day.The histogram of R is almost symmetrical about the vertical axis through the value of 6-7 nT/day.Although these histograms contain much information on the process under study, we won't comment on them any longer in this study of predictability.

Autocorrelation of R
The autocorrelation of the R(k) series can be illustrated by a few diagrams.Figure 2 displays the 2-D histogram of pairs of consecutive values of R, (R(k),R(k + 1)), with the bins defined by a = 10 0.05 .The horizontal axis is for the "today" value and the vertical one for the "tomorrow" value.The elongation of the cloud of points is indicative of a correlation in the time series.The histogram, when normalized to the total number of pairs (N − 1 = 31 045), delivers their empirical distribution, which can be used to determine various conditional probabilities.
Figure 3 displays, as an example, the graphs of two of the conditional probabilities: i.e. probability of having R greater than 50 nT/day the day after the value R is observed, and i.e. the probability that, given the value of R(k) today, its value will be larger tomorrow.The first graph shows that (i) there is practically no chance for R to reach 50 the day after its current level is less than 15, (ii) if its current level is 120 or more, there are more chances for R to stay above 50 the day after than to fall below this level (note, however, that the number of high value extreme events is small, so that the statistic is less robust in this range and that the probability of the extreme values eventually collapses to 0).The second graph also illustrates some kind of dynamical law of the system: if R < 8 today, chances are higher than 50% for R to be larger the day after; if R > 80, chances for R to be larger are limited to 20%-30%.As already mentioned, other condi- , and the empirical distribution of R (blue line), i.e. unconditional probability tional probabilities can be derived from the (R(k),R(k The autocorrelation function of R, where is average and t is the lag time.The Correlation is 0.54 after one day (t = 1), but falls down to 0.22 after two days.

A simple prediction scheme
Suppose we are interested in predicting today "extreme" events, defined by values of R > R 0 , for tomorrow.The simplest prediction scheme suggests issuing an alert for tomorrow if R > r 0 today; r 0 and R 0 are parameters.We call the occurrence of a value of R larger than R 0 an extreme event.We count a success if it happens on the alert day, and a failure-to-predict otherwise.
Figure 4 shows the error diagram, or so-called (n, τ ) diagram (Molchan, 1997), achieved by the described prediction scheme on Eskdalemuir data, for a large set of values of parameters r 0 and R 0 .n is the percentage of unpredicted (extreme) events, τ is the ratio of the time covered by alerts, i.e. of the ratio of the number of alert days to the total time interval considered (N days).The effectiveness of the prediction is characterized by the distance of the lower envelope of the set of points (n, τ ), from the random guess strategy curve, i.e. the segment of diagonal n + τ = 1.This diagonal segment connects the point corresponding to the optimistic strategy (no alert, and failure to predict any event) to the point corresponding to the pessimistic strategy (full time alert, and no failure to predict).In general, one tries to minimize some loss function γ (n, τ ) depending on preparedness problems and measures envisioned in response to the prediction.The point where an isoline of γ (red line in Fig. 4) touches determines both the minimal achievable loss and the optimal set of adjustable parameters (here r 0 and R 0 ) of the prediction algorithm.
Here, for illustration, we adopt the linear cost function γ = n + τ (note that γ = 100% at any point of the random guess diagonal).
To facilitate further comparison we define here the extreme events with R 0 between 98 and 99 percentiles of R, i.e. the top 1-2% in a data set, and r 0 between 0 and 98 percentile.In case of the Eskdalemuir daily data, 41 < R 0 < 55 nT/day and 0 < r 0 < 41 nT/day.Each of the 800 (= 20×40 , arbitrary sampling of R 0 , r 0 variations) pairs (R 0 , r 0 ) that are uniformly distributed in the above defined intervals is mapped on the error diagram according to the score achieved by the prediction scheme with these parameters, as seen in Fig. 4. The lower envelope of the mapping of the domain spanned by (R 0 , r 0 ) pairs is used to determine the optimal parameters for prediction.We draw level lines of the cost function γ (segments of a straight line of slope −1 in the case of γ = n + τ ); the value attached to the one which is tangent to is the minimum of the cost function γ .In the case of γ = n+τ (e.g.Fig. 4), this value can be read directly on either axis.
For Eskdalemuir (Fig. 4), the minimal value of γ = 43% is achieved (in such definition of the R extremes), when r 0 = 11.3 and R 0 = 51.5.Twenty-two percent of alert days, i.e. 18.6 out of 85 years, is required to predict 287 out of 362 extremes.The score is far from a random prediction and compares with the case of reproducible earthquake prediction, which came from the 10 years of real-time global testing of M8-MSc algorithm with γ = 34% ( Keilis-Borok et al., 2001).Of course, one may want to use a more restrictive definition of the extremes and issue a smaller number of alerts.For example, for R 0 = 120, and r 0 = 60, 25 out of 40 extreme events are predicted by issuing 254 days of alert, i.e. about 0.8% of the total number of days (γ = 38%).

Hourly first differences
We now consider hourly mean values of X, Y , Z at Chambon-la-Forêt observatory (Table 1) in the 1974-85 timespan, covering a full solar cycle (note that we retained 12 years for the study of CLF hourly means, which cover the solar cycle number 21, and 11 years for the study of PAF hourly variations from minute values (Sect.3.3), which span the cycle number 22).We define, exactly as above, Ẋ(k), Ẏ (k), Ż(k), k now being the number of the sequential hour and the dot meaning the difference E(k + 1) − E(k) (Eq.1).R(k) is given again by Eq. ( 2).However, now Ė and R are measured in nT per hour.
Figure 5  be larger than R(k), R(k) given.It appears that this probability stabilizes after 60 nT/hour, but eventually collapses for the most extreme values (> 300 nT/hour).
We have also computed the error diagram (Fig. 7) for the simplest prediction scheme aimed at the extreme events defined as above (Sect.3.1.3)in the range from 2 to 1% of the top values of R. The minimal γ = n + τ = 42% is reached when R 0 = 49.4 and r 0 = 17.1 nT/hour.If we take R 0 = 112 and r 0 = 24 nT/hour, then γ = 24%, which corresponds to 8.1% of alert and prediction of 249 out of 295 extreme events (i.e.84.1%).Again, the choice of parameters r 0 , R 0 depends on what we want to do with the prediction.

Hourly indices from minute values
The last two series considered are derived from minute values recorded at the Chambon-la-Forêt (France) and Port-aux-Franc ¸ais (Kerguelen Islands; see Table 1) observatories.The PAF data cover a full Solar cycle, i.e. 1985-1995 and the CLF one the period 1992-1995.
We compute the first difference (Eq.1), Ė(k), for each minute with the sequential number k (5.77 • 10 6 values at Port-aux-Franc ¸ais, 2.10 • 10 6 values at Chambon-la-Forêt).Accordingly, the unit measure is nT/min.From the minute differences relative to the hour interval with sequential number m, we compute the average: where  m.In this definition R * can be viewed as an extension of the mathematical total variation of a real-value function that characterizes the volatility of the process (in the sense used in financial markets).R * (m) can also be viewed as an invariant hourly magnetic activity index computed from minute data.We will process R * the same way as we did for the hourly first difference R(m) in the previous section.Figures 8 and 9 present the 2-D histograms (R * (m),R * (m + 1)) determined from 95 000 PAF and 35 000 CLF data points.They are clearly more elongated and narrower than the previous ones (Figs. 2 and 5), indicating better correlation, which implies better predictability.
The conditional probabilities shown in Figs. 10 and 11 were computed from the 2-D histograms (Figs. 8 and 9).One may conclude from the graphs that (i) there is practically no chance for R * observed at PAF to reach 50 nT/min in the next hour if its current level is less than 10 nT/min (Fig. 10 P > 50%) to record an increase in R * during the next hour, if its present level is 1 nT/min (Figs. 10 and 11,bottom).The level of 1 nT/min or less is observed at PAF 55% of time, while at CLF it happens 80% of time (see the empirical distributions of R * , P(R * (m) < R * ), plotted in the same figures).
The probability graph P(R * (m + 1) > R * (m)) for CLF (Fig. 11 bottom) is much steeper and collapses to nearly no chance at 10 nT/min; for PAF, it extends above 100 nT/min.As already said, the two graphs P(R * (m + 1) > R * (m)) and P(R * (m) < R * ) determine the dynamics of the system and can be used for a proper statistical simulation of magnetic indices at a given location.
Figures 12 and 13 illustrate the predictability of the extreme events defined by R * and how it differs at the two sites: the minimal value of n + τ is much smaller than in the previous two cases.Table 2 summarizes the results of predicting magnetic field extremes.Specifically, n + τ for the magnetic rate prediction at CLF is reduced by a factor of 2 when R * is used instead of R.This score is outperformed for another factor of 2 by the prediction of the R * extremes at PAF.The predictability of geomagnetic series by the simple scheme was comparable in the first two cases, whereas it is by far better in the latter two cases than the reproducible intermediate-term prediction of the largest earthquakes (Keilis-Borok et al., 2001).
Finally, we show how the magnetic extreme events can be predicted from the hourly values of the magnetic variation, R * , several hours in advance.We consider the 11-year series from PAF, define extreme events values as R 0 = 20 nT/min, and compute the error diagrams for the simple prediction scheme applied with a lead-time t = 1, 2, 3, and 4 h (Fig. 14).Results are summarized in Table 3. Naturally, the optimal value of the cost function γ = n + τ is growing with lead time, although at a decreasing rate: from a factor of 1.6 from 1 h to 2 h and 1.1 from 3 h to 4 h.Similarly, the optimal alert threshold decreases from 5.0 nT/min for a lead time of 1 h to 3.0 for 4 h.We also give in Table 3 a decomposition of γ = n + τ presenting the percentages of failures-to-predict n in the whole set of 823 extreme events and of the relative alarm time over the 11 years (= 95 474 h) considered.

Discussion and conclusion
We have investigated the predictability contained in geomagnetic data using unusually long series, which allow for a firm statistical analysis.Two samplings have been used: daily and hourly, with two different variables in the last case.The simple tools we use allowed us to quantify in different manners the predictability of the series.The simplest of them focus on what can be said about the next value of the variable knowing its current value, e.g.what can be said for tomorrow knowing the situation for today.We have chosen two examples of conditional probabilities; many other ones can be computed in the same way, making use of the statistically firm empirical distributions.In every case, the answer is clear and rigorously established.The error diagram of our simple prediction scheme delivers comprehensible, data supported answers, in terms of chances of failures-to-predict and percentage of alert time, to questions concerning costs and benefits, which can be formulated in various ways.The simple tools used are limited in their forecasting capacity.Nevertheless, the same techniques can be easily generalized at will: for example, without changing the structure of the algorithm, tomorrow can be replaced by the day after tomorrow and so on.The functional form for the γ cost function has to be determined for each practical situation.In this paper, for illustration, we assumed that the cost (e.g.financial cost) of a failure to predict was the same as the cost of a false alarm.It is likely that a failure to predict would have a larger cost since a failure to predict may lead to damages, whereas a false alarm may only induce a loss of profit.In such a case, a greater weight should be attributed to the percentage of unpredicted extreme events (n) in the cost function.
The classical autocorrelation functions of R and R * drop quickly, at least by a factor of 2 in the first two days for daily series and by a factor of 2 in the first 4-12 h for hourly se-ries.The correlations and predictability may depend on geomagnetic location: the first are higher at CLF than at PAF, whereas the simple prediction scheme is much more efficient at PAF than at CLF.The magnetic vector rate is more predictable when measured by the R * index.However, Bellanger et al. (2002b) andBellanger et al. (2002a) showed that R * is almost identical, within a constant factor, in all lowand mid-latitude observatories; R * at CLF, as considered in this paper, is thus expected to describe the predictability of magnetic series at most of the surface of the Earth.PAF, an observatory in the auroral zone, has been studied to allow comparison.Bellanger et al. (2002b,a) also showed that R * characterizes the variation of the external field and can thus be considered as an activity index.In the present paper, R * has been preferred to any other index because it is an invariant characteristic of the magnetic field activity whose time sampling can be easily changed by adjusting the width of the averaging window (one hour in this paper); but considering, for example, aa or am (3-hour range indices, see Mayaud, 1980) would have given similar results.
We did not search for geomagnetic precursors stricto sensu.The data suggests that it is dubious that a peculiar magnetic variation observed in quiet time might warn of an approaching sudden burst of activity, like a disastrous magnetic storm.On the other hand, recurrences of enhanced probability of magnetic activity induced by the Sun's rotation of about 27 days and the solar cycle of about 11 years are well known.
Hopes in predicting sudden bursts of activity rather rely upon the real-time observations of solar wind and coronal mass ejections (CMEs) (Joselyn, 1995).However, due to the complexity of the involved physical processes (e.g.Boaghe et al., 2001), no complete quantitative theory of the magnetospheric dynamics is available at the present time, and thus no fully reliable prediction of magnetic activity is possible.Moreover, the efficiency of forecasts of geomagnetic activity from solar and interplanetary conditions is not systematically estimated a posteriori (Thomson, 2000), although prospective and retrospective validation is applied in other fields of geophysics (Mulargia, 1997;Kossobokov et al., 1999).
Our aim was to lead a practical, rational and rigorous approach concerning what could be done, based on the knowledge of magnetic series, in the field of extreme geomagnetic activity events prediction.Space weather extreme conditions have an important financial impact on a wide domain of activities (see, e.g.Allen & Wilkinson, 1993;Maynard, 1995, for a summary).The prediction scheme studied here, via the adjustment of the γ function to each specific customer, depending on the cost of an alarm vs. the cost of a failure to predict an extreme event, provides a quantitative tool in a decision theory perspective (Matthews, 1997;Keilis-Borok et al., 2001;Thomson, 2000).Despite that this method is far from covering all the needs of the large variety of customers interested in space weather forecasting (Feynman and Gabriel, 2000), it gives at least a simple, statistically robust and quantitative way for short-term prediction of geomag-   3). 1 h lag in red, 2 h lag in green, 3 h lag in blue and 4 h lag in violet.
displays the distributions of the absolute values | Ẋ|, | Ẏ | and | Ż| and of R for the whole pe-

Fig. 5 .
Fig. 5. R for next hour vs R for current hour at CLF (computed from hourly means): 2-D histogram.The color indicates the number of pairs in a cell, see color bar.

Fig. 9 .
Fig. 9. R * for next hour vs R * for current hour computed from minute values at CLF: 2-D histogram.

Table 1 .
Location of observatories: geographic (geocentric) coordinates and corrected geomagnetic (CGM) coordinates in degrees D histogram.The color indicates the number of pairs in a cell, see color bar.

Table 2 .
Summary of the simplest prediction of R

Table 3 .
Optimal threshold for issuing an alarm r 0 , cost value γ = n + τ , failure ratio n and number of failures to predict (N f ), alarm ratio τ and total hours of alarm (N a ) for an extreme event threshold R 0 fixed to 20 nT/min and for a prediction with a lead time t The error diagram for PAF (hourly variation R * computed from minute values).The error diagram for CLF (hourly variation R * computed from minute values).Bellanger, E., Anad, F., Le Mouël, J.-L., and Mandea, M.: The irregular variations of the external geomagnetic field from Intermagnet data, Earth Planets Space, submitted, 2002a.Bellanger, E., Blanter, E. M., Le Mouël, J.-L., Mandea, M., and Shnirman, M. G.: On the geometry of the external geomagnetic irregular variations, J. Geophys.Res., 107, 10.1029/2001JA900112, 2002b.Boaghe, O. M., Balikhin, M. A., Billings, S. A., and Alleyne, H.: Identification of nonlinear processes in the magnetospheric dynamics and forecasting of D st index, J. Geophys.Res., Error diagram for PAF (hourly variation R * computed from minute values) up to a 4-hour prediction (see text and Table n Fig. 12. n Fig. 13.n Fig. 14.