FINANCE , STATISTICS, FORECASTING. ACCURACY TESTING AND ERROR MATRICES

Contents lists available at ScienceDirect

New Astronomy

journal homepage: www.elsevier.com/locate/newast

Forecasting space weather over short horizons: Revised and updated estimates

Gordon Reikard Statistics Department, U.S. Cellular, 8410 West Brywn Mawr, Chicago, IL 60631, United States

A R T I C L E I N F O

Keywords: Space weather Forecasting Frequency domain models Time series models

A B S T R A C T

Space weather reflects multiple causes. There is a clear influence for the sun on the near-earth environment. Solar activity shows evidence of chaotic properties, implying that prediction may be limited beyond short horizons. At the same time, geomagnetic activity also reflects the rotation of the earth's core, and local currents in the ionosphere. The combination of influences means that geomagnetic indexes behave like multifractals, exhibiting nonlinear variability, with intermittent outliers. This study tests a range of models: regressions, neural networks, and a frequency domain algorithm. Forecasting tests are run for sunspots and irradiance from 1820 onward, for the Aa geomagnetic index from 1868 onward, and the Am index from 1959 onward, over horizons of 1–7 days. For irradiance and sunspots, persistence actually does better over short horizons. None of the other models really dominate. For the geomagnetic indexes, the persistence method does badly, while the neural net also shows large errors. The remaining models all achieve about the same level of accuracy. The errors are in the range of 48% at 1 day, and 54% at all later horizons. Additional tests are run over horizons of 1–4 weeks. At 1 week, the best models reduce the error to about 35%. Over horizons of four weeks, the model errors increase. The findings are somewhat pessimistic. Over short horizons, geomagnetic activity exhibits so much random variation that the forecast errors are extremely high. Over slightly longer horizons, there is some improvement from estimating in the frequency domain, but not a great deal. Including solar activity in the models does not yield any improvement in accuracy.

1. Introduction

Forecasts for space weather are used in several areas, ranging from electric power to navigation and communications systems. The hor- izons are generally short, from as little as a few hours to as long as several days. Geomagnetic series show evidence of being complex processes. There is a clear influence for solar activity. Geomagnetic indexes exhibit cycles at 27 days, corresponding to the rotation of the sun, and the 11–22 year Schwabe and Hale cycles, corresponding to reversals in the sun's magnetic field, although they are slightly de- phased (Russell and McPherron, 1978). High magnetic activity during the twentieth century corresponded with a period in which solar ac- tivity was at its highest level in 11 millennia (Solanki et al. 2004, Frohlich, 2009). At the same time, there are periods in which the geomagnetic series are not closely correlated with the solar cycle, no- tably the 1970s (Rangarajan and Iyemori, 1997).

Geomagnetic activity also reflects the rotation of the earth's core, and local currents in the ionosphere (Jordan, 1979; Buffett, 2000; Weiss, 2002). The solid inner core, with a radius of about 1220 km, rotates within a liquid outer core, extending out to 3400 km. The

circulation of the outer core is driven by the transfer of heat from the inner core, to the core-mantle boundary. In effect, the magnetic field can be modeled as a dynamo: the change in the magnetic field depends on the diffusivity and the velocity of the fluid. Changes in magnetism generate an electric field, while electric currents give rise to magnetic fields. The magnetic field in the outer core has been estimated as roughly 50 times more intense than at the surface (Buffett, 2010).

The combination of these two processes – solar activity and the internal dynamo – implies that the geomagnetic series may be difficult to predict. The interaction two or more stochastic processes implies that geomagnetic activity should have multifractal properties (Lovejoy and Schertzer, 2013). The solar influence in and of itself may be difficult to forecast beyond very short horizons, since solar activity shows evidence of chaos (Hansen and Willson, 1997; Solanki and Krivova, 2011; Feynman and Ruzmaikin, 2011).

The literature on space weather forecasting includes several methods, ranging from regressions to neural networks and artificial intelligence. The approach tested here includes both frequency and time domain methods, and using solar activity as an input in models for geomagnetic indexes. The findings revise and update earlier results in

https://doi.org/10.1016/j.newast.2018.01.009 Received 30 November 2017; Received in revised form 8 January 2018; Accepted 17 January 2018

E-mail address: [email protected].

New Astronomy 62 (2018) 62–69

Available online 31 January 2018 1384-1076/ © 2018 Elsevier B.V. All rights reserved.

T

Reikard (2011, 2013), which were later determined to have been biased by programming errors. The programming errors were not present in and did not affect the results in Reikard (2015).

Section 2 reviews the data. Section 3 sets out the forecasting models. Section 4 runs forecasting experiments. Further analysis is presented in Section 5, and Section 6 concludes.

2. The data

Sunspots have been directly observed since the seventeenth century, and daily sunspot numbers have been reconstructed by interpolating gaps in the actual observations. The sunspots data was obtained from the long-term solar observations website, http://www.sidc.be/silso/ datafiles, maintained by the Royal Observatory of Belgium, and sup- ported by the International Council for Science. The dataset used was the revised series developed in Clette et al., (2014), starting in 1818.

Since late 1978, total solar irradiance has been directly measured by satellites, and is available in the ACRIM database (ACRIM, 2014; Willson and Mordvinov, 2003; Scafetta and Willson, 2014). Daily ir- radiance has been reconstructed to as far back as 1610 (Krivova et al., 2010; Vieira et al., 2011). The values for irradiance were downloaded from a website maintained by a consortium of laboratories and uni- versities in several countries: http://www.mps.mpg.de/projects/sun- climate/data.html. The data set used here is the SATIRE-H series for total solar irradiance, and is updated to 2017 using the ACRIM series. Irradiance is in watts per meter squared (W/m2).

Figs. 1 and 2 show 30-year cross sections of the sunspot and irra- diance data sets. Both of them show several characteristic features, the 27-day cycle, the 11–22 year cycles, and an upward trend in the early part of the twentieth century. Both also show nonlinear variability, with irregular outliers at intermittent intervals.

Two geomagnetic data sets were used. The Aa index, which begins in 1868, is based on the magnetic activity measured at two antipodal stations, Canberra, Australia, and Hartland, England (Mayaud, 1972). The index is the average of the northern and southern values of mag- netic activity, weighted to account for differences in the latitudes of the two stations and local induction effects. The Aa activity has the ad- vantage of being the longest record of geomagnetic activity.

The Am and Ap indexes were developed from the K index, a 3-hour measure of activity relative to an estimated quiet day. It was introduced at the beginning of the 1940s (Mayaud, 1980; Menvielle and Berthelier,

1991), and extrapolated backwards using data from selected magnetic observatories. The Ap index, beginning in 1932, uses 13 observatories. The Am index, beginning in 1959, spans 23 locations, 13 northern and 10 southern stations, arranged in groups representing longitudinal sectors (Menvielle and Berthelier, 1991; McPherron, 1995; Mayaud, 1980). Because the Am index spans a wider range of observation points, it is used in the following tests. The data prior to 2011 were down- loaded from the National Geophysical Data Center (NGDC, 2014). Data from 2011 onward were obtained from the International Service of Geomagnetic Indexes (http://isgi.unistra.fr/data_download.php). The data run through May 31, 2017.

Figs. 3 and 4 show short cross-sections of the Aa and Am indexes. Both show high degrees of nonlinear variability. During periods of high solar activity, the data is more volatile than at the trough of the solar cycle.

3. The forecasting models

Reviews of the literature are available in McPherron (1995) and Lundstedt (2005). Many of the earlier studies focused on a single class of models. Regression-based methods were used in Baker et al., (1990), Blanchard and McPherron (1993, 1995), Papitashvili et al., (1998), and Codrescu et al., (2004). One limitation of regressions is that they tend to capture the central tendency of the data more effectively than the turning points. In part for this reason, many studies have preferred neural networks. References include Lundstedt (1996), Gleisner and Lundstedt (1997), Wu and Lundstedt (1997a,b), O'Brien and McPherron (2000), Tulunay et al., (2005), Martin et al., (2005), Vandegriff et al., (2005), and Wang et al., (2008).

Other artificial intelligence techniques, including genetic algo- rithms, have been used in McPherron (1993), Wintoft and Lundstedt (1998), Orfila et al., (2002), Mirmomeni et al., (2006, 2007), Vorotnikov et al., (2008). Mirmomeni et al., (2007) use co-evolving models. Mirmomeni et al., (2010, 2011) and Gholipour et al., (2007) combine neuro-fuzzy techniques with spectral methods (see also Kalhor et al., 2011, 2012). Cao and Cao (2006) combine neural nets and wa- velets.

Despite the popularity of more advanced techniques, regressions with time-varying coefficients have often proven very effective in dealing with nonlinear data (Bunn, 2004; Granger, 2008). Let Yt denote geomagnetic activity, ln denote natural logs, ω denote a coefficient, the

1960 1962 1964 1966 1968 1970 1972 1974 1976 1978 1980 1982 1984 1986 1988 0

50

100

150

200

250

300

Fig. 1. Daily sunspot number, January 1, 1960 through December 31, 1989.

G. Reikard New Astronomy 62 (2018) 62–69

63

subscript t denote time variation and ε denote a random error. The lags at 1–3 days capture proximate dependence in the time series. The lag at 27 days captures the rotation of the sun.

= + + + + +

∼

− − − −ω ω ω ω ω

σ

lnY lnY lnY lnY lnY ɛ ɛ

P (0, ) t 0t 1t t 1 2t t 2 3t t 3 4t t 27 t t

t 2 (1)

where P is the probability distribution, assumed non-Gaussian, and σt2 is the residual variance. This specification was preferred by the Akaike (1973) information criterion. This model can be modified by using fewer autoregressive lags, and including exogenous terms like irradiance and sunspots.

Frequency domain models are attractive because of their ability to capture cyclicality and long memory (Fox and Taqqu, 1986). The fre- quency domain method used here involves Fourier transformation of the moving average representation, which is familiar from the well- known ARIMA class of models (Box and Jenkins, 1976). Using standard

notation, let ϕ(L) be the autoregressive operator, represented as a polynomial in the backshift operator: ϕ(L) = 1 − ϕ1 L − ⋅⋅⋅ − ϕpL

p. Let θ(L) be the moving average operator: θ(L) = 1 + θ1 L + ⋅⋅⋅ + θqLq. Let the superscript ξ denote the order of differencing (equivalently, the order of integration). The model is then of the form:

− = +μ θ ϕ(1 L) lnY [ (L)/ (L) ] ɛξ t t t t (2)

where μ is the intercept. The moving average representation, denoted β, is defined as the ratio of the moving average and autoregressive poly- nomials:

=β θ ϕ(L) ɛ [ (L)/ (L) ] ɛt t t t (3)

Following Koopmans (1974), the spectral density (F) can be ex- pressed as a function of the z-transform:

= −β β σF (z) (z )Yt 1 t 2 (4)

The z-transform converts the discrete values of a time series into a

1960 1962 1964 1966 1968 1970 1972 1974 1976 1978 1980 1982 1984 1986 1988 1363

1364

1365

1366

1367

1368

1369

Fig. 2. Total solar irradiance, reconstructed from sunspots and satellite observations, in watts per meter squared. January 1, 1960 through December 31, 1989.

1960 1962 1964 1966 1968 1970 1972 1974 1976 1978 1980 1982 1984 1986 1988 0

50

100

150

200

250

300

350

400

Fig. 3. The Aa geomagnetic index, January 1, 1960 through December 31, 1989.

G. Reikard New Astronomy 62 (2018) 62–69

64

frequency domain representation, and is given by: z = Λ exp (iφ) = Λ(cos φ + i sin φ), where Λ is the magnitude of z, i is the imaginary unit, and φ is the complex argument, in radians. Let γ denote a one-sided polynomial in positive powers of z. Then the log spectral density can be expressed as:

= + + −ln γ γ γF (z) (z )Yt 1 0 (5)

Taking antilogs, and combining (4) with (5):

= − −β β σ γ γ γ(z) (z ) exp [ (z) ] exp [ (z ) ] exp ( )1 t 2 1 0 (6)

Eq. (6) provides a frequency domain estimate of β, which can be projected outside the range of the series. The forecast in the time do- main is then recovered by inverse Fourier transformation. For other frequency domain techniques, see Lotsukov et al. (2001), Steinhilber and Beer (2013).

In many fields, outliers are treated as anomalous and smoothed out (Chen and Liu, 1993). However, in space weather forecasting, it is often necessary to predict the outlying values, since solar storms can disrupt the earth's magnetic field. Accordingly, it is of interest to use models which take explicit account of the outliers, using state transition terms (Hamilton, 1989, 1990). One way to do this is to create a binary state variable, denoted St, such that St = 1 when magnetic activity is in a peak state and 0 otherwise. The state term can enter the equation ad- ditively or multiplicatively, or both. Let It, denote irradiance. The equation for the daily data is then of the form:

= + + + + +

+ + +

− − − −

−

ω ω ω ω ω ω

ω ω

lnY lnY lnY lnY lnY lnI

S S lnY ɛ t 0t 1t t 1 2t t 2 3t t 3 4t t 27 5t t

6t t 7t t t 1 t (7)

The state must of course be predicted. One way to do this is a lo- gistic regression. The model is of the form:

= + + + +− −ω ω ω ωS lnY lnI S ɛt 0t 1t t 1 2t t 3t t 1 t (8)

This relates the state to recent values of the time series, prior values of the state, and solar activity.

It is also possible to set up models using neural networks. The system architecture used here is fairly standard. The neural net consists of a multilayer perceptron, with one direct connection and three hidden layers, trained using a backpropagation algorithm. Additional layers were tested, but these increased solution time with no gain in accuracy. The net was trained by epoch, rather than by example, since this was found to result in more accurate predictions. The relative merits of

regressions and neural nets are difficult to assess in advance. Neural nets have the advantage that they can capture complex patterns, but can be vulnerable to over-fitting noise.

4. Forecasting experiments

One issue in running forecasting experiments is the length of the interval to be covered. Although it would be possible to begin the tests for irradiance as early as 1610, during the Maunder minimum of 1645–1715 and the subsequent Dalton minimum of 1790–1820, the sun was quiescent, and highly predictable, which would bias the error downward. Instead, the forecasting tests for irradiance and sunspots begin in 1820. The tests for the two geomagnetic indexes begin at the inception points of the data sets, in 1868 and 1959 respectively.

In forecasting experiments, a common practice is to use a persis- tence model, i.e., setting the predicted value equal to the most recent observations, as a benchmark, and this is done here. The model fore- casts were run iteratively. In other words, for each observation in the series, the model was estimated over a given range of prior values, and then forecasted. At the next step, the model was re-estimated over the most recent value, while moving the first observation in the data set one period ahead, through the entire sample. All the forecasting is dynamic, and all predictions are true out-of-sample forecasts. In the tests for horizons beyond one day, the forecasts skip the intervening days. In other words, for t + 2, the forecast error is calculated only for t + 2; t + 1 is omitted.

The time-varying coefficients were estimated using moving win- dows. This was found to be preferable to a Kalman filter (Kalman, 1960). Initial experiments demonstrated that when an un- restricted Kalman filter was used, the coefficients exhibited too much randomness, which tended to reduce forecast accuracy. While it would be possible to impose restrictions on the Kalman filter, the moving window has the advantage that changes in the width can be used to control the rate of coefficient variation (Rossi and Inoue, 2012). Tests were run for a range of window widths. Over shorter widths, the coefficients become too volatile, reducing predictive accuracy. As the width expands, the coefficients stabilize, and forecast accuracy in- creases. At extremely wide widths, on the order of more than a year, the coefficients become too inertial, and forecast accuracy declines. The preferred window width was 270 days, encompassing 10 full solar ro- tations.

1960 1962 1964 1966 1968 1970 1972 1974 1976 1978 1980 1982 1984 1986 1988 0

50

100

150

200

250

300

350

400

Fig. 4. The Am geomagnetic index, January 1, 1960 through December 31, 1989.

G. Reikard New Astronomy 62 (2018) 62–69

65

The first set of tests is for irradiance and sunspots. The models were trained over the first 1500 observations, and then run forward from this point. The measure of forecast accuracy in this instance is the mean absolute error. The errors are reported both for all observations, and for the outliers. The error in units is preferable to the percent error in these data sets. In the case of total solar irradiance, the base value is so high that the percent error will be very low. By comparison, the number of sunspots can fall to zero at the trough of the 11-year cycle, making the percent error indeterminate. Even when the low values are non-zero, the percent error will tend to be higher at the cycle trough and lower at the cycle peak.

Model 1 is the persistence forecast. Model 2 uses a regression with lags over the first three days, and at 27 days. Model 3 uses the state transition terms. Model 4 is the neural net with the same inputs as in the regression. Model 5 is the frequency domain algorithm. In the state transition model, the high state was defined as values above 200 sun- spots, or irradiance above 1367 W/m 2. In the models for the state, the logistic regression and neural net were estimated over all prior ob- servations.

The findings are reported in Table 1. For irradiance, the persistence forecast consistently shows the lowest errors. None of the other models are particularly effective. The state transition model is only useful at the 1-day horizon. Beyond this, it is nearly impossible to predict the state correctly, so the model is no better than the alternatives. The neural network consistently yields the highest errors. The frequency domain method achieves roughly the same degree of accuracy as the time do- main models.

For sunspots, the persistence forecast is better only over the first three days. At horizons beyond four days, the other models do better. The neural net again shows the least accurate forecast. The other three methods are nearly tied, all better than persistence, but none really dominating the others.

Table 2 reports the results for the geomagnetic indexes. The per- sistence method fails badly, consistently producing the least accurate forecasts, with errors over 50% at the 1-day horizon, rising to over 70% at all horizons beyond a day. The neural net again does poorly. The input and bias weights vary so much over the moving windows that the model over-fits the noise. The remaining models all achieve about the same level of accuracy. However, the forecasts are not particularly accurate. The errors are in the range of 48% at 1 day, and 54% at all later horizons. Again, neither the time or frequency domain models are necessarily superior to the others. The state transition model is not particularly effective because of the difficulty in predicting the states. If the states could be known perfectly, the forecast error would decline by 1–2%. However, since the states cannot be predicted accurately beyond

one day, incorporating the transition terms in the model fails to im- prove the forecast.

The next set of experiments is for horizons of 1–4 weeks. To predict over these longer intervals, the data is converted to weekly values. Table 3 reports the findings for solar activity. The regressions and fre- quency domain algorithm both achieve lower errors than the persis- tence model, at all horizons. The neural network consistently produces higher errors.

Table 4 reports the findings for the geomagnetic indexes. The per- sistence errors are still very high, over 45% at the 1 week horizon. Surprisingly, the persistence errors decrease as the horizon increases, falling to 38% at four weeks. At the 1-week horizon, the regressions and frequency domain algorithm reduce the forecast errors to the 35% range. It should be noted that this is much more accurate than the forecasts from the same models using daily data. Again, including ir- radiance and sunspots fails to achieve any improvement in accuracy. Instead, the errors are marginally worse. The neural net is still less accurate than the regressions, but is better than the persistence method. As the horizon increases, the forecast error gradually increases. At four

Table 1 Forecast accuracy, solar activity.

Statistics are the mean absolute error.

Model Horizon (days)

1 2 3 4 5 6 7

Part 1: Irradiance Persistence 0.081 0.132 0.176 0.212 0.237 0.257 0.270 Regression on lags 0.104 0.175 0.230 0.270 0.300 0.319 0.341 Regression, state

transition terms 0.096 0.160 0.207 0.263 0.297 0.303 0.329

Neural network 0.114 0.217 0.363 0.368 0.371 0.372 0.388 Frequency domain 0.090 0.159 0.233 0.265 0.286 0.302 0.324 Part 2: Sunspots Persistence 13.10 19.21 24.45 29.01 32.56 35.78 38.54 Regression on lags 14.84 21.25 25.47 28.36 30.35 31.67 32.53 Regression, state

transition terms 15.08 21.51 25.71 28.53 30.47 31.78 32.65

Neural network 16.39 35.34 36.37 36.57 36.57 43.37 43.57 Frequency domain 14.94 21.06 25.28 26.25 30.05 31.59 32.52

Table 2 Forecast accuracy, geomagnetic indexes.

Statistics are the mean absolute percent error.

Model Horizon (days)

1 2 3 4 5 6 7

Part 1: Aa Index Persistence 55.24 72.24 78.51 80.85 81.14 81.82 81.96 Regression on lags 48.11 54.59 55.45 55.35 55.57 55.62 55.64 Regression, lags and

sunspots 48.24 54.65 55.77 55.91 55.96 56.01 56.04

Regression, lags and irradiance

48.23 54.66 55.76 55.86 55.91 55.99 56.01

Regression, state transition terms

48.13 54.68 55.56 55.65 55.68 55.72 55.75

Neural network 58.25 65.15 66.75 67.82 69.92 70.25 71.35 Frequency domain 48.12 54.57 55.68 55.77 55.64 55.83 55.84 Part 2: Am Index Persistence 54.83 72.41 79.18 81.62 82.02 81.65 81.51 Regression on lags 47.59 54.45 55.46 55.57 55.61 55.65 55.69 Regression, lags and

sunspots 47.73 54.61 55.67 55.81 55.81 55.69 55.94

Regression, lags and irradiance

47.71 54.67 55.66 55.84 55.98 56.02 55.91

Regression, state transition terms

47.63 54.51 55.68 55.81 55.79 55.81 55.83

Neural network 57.15 65.48 66.37 67.19 69.78 70.15 71.27 Frequency domain 47.58 54.71 55.67 55.43 55.81 55.79 55.88

Table 3 Forecast accuracy, solar activity.

Statistics are the mean absolute error.

Model Horizon (weeks)

1 2 3 4

Part 1: Irradiance Persistence 0.203 0.246 0.230 0.220 Regression on lags 0.203 0.225 0.233 0.239 Regression, state transition terms 0.202 0.227 0.239 0.241 Neural network 0.276 0.275 0.275 0.276 Frequency domain 0.206 0.217 0.230 0.237 Part 2: Sunspots Persistence 30.07 37.97 34.01 31.45 Regression on lags 27.70 30.98 32.09 33.04 Regression, state transition terms 27.68 30.99 32.11 33.07 Neural network 36.66 38.62 40.53 43.06 Frequency domain 27.59 30.30 30.89 30.95

G. Reikard New Astronomy 62 (2018) 62–69

66

weeks, the neural net error is higher than the persistence error, while the regression errors remain lower. The frequency domain model is slightly better than time domain methods at this horizon.

5. Further analysis

Some of the findings are counterintuitive. The geomagnetic indexes show some relationship to solar activity, but including measures of solar activity in the equations for space weather results in higher errors. To assess the relationship between geomagnetism and solar activity, the indexes were regressed on irradiance and sunspots, in logs. Table 5 presents the estimated regressions. When the regressions were run in levels, they showed evidence of high degrees of serial correlation. For this reason, they were re-estimated in fractional differences, a standard technique for dealing with serial dependence. The elasticities are all significant, and in some instances large, although this is due mainly to differences in the scale of the data. However, the share of the variance explained by solar activity is quite small. The R-bar-square is only 0.365 for the Aa index at a daily resolution, and 0.354 for the Am index. At this resolution, solar activity is found to account for at most only about one-third of the variation in the geomagnetic indexes. At the weekly resolution, the R-bar-square is even lower, 0.234 to 0.154.

The neural network performs remarkably poorly on these data sets. One reason is the moving window. In many prior studies, the net was trained over a range of observations, and then cross-validated on the remainder. This however does not closely approximate how a forecaster would operate in the real world, where the task is to regularly update forecasts as new data becomes available. Instead, in this study, the neural net is trained over the previous 270 days, so that the input and bias weights vary in the same way as the regression coefficients. When specified this way, the net is extremely vulnerable to over-fitting noise.

In this respect, prior studies have noted similar issues with neural networks trained over limited data sets (Lu and Wang, 2005; Papaleonidas and Iliadis, 2013).

A third finding is that the models do not predict the daily geo- magnetic data accurately, but are somewhat better at a weekly re- solution. To determine whether the models were missing signal in the data, the residuals from the forecasting models were taken. Regressing the residuals on lags finds little evidence of any structure. There is a slight negative relationship at one day, and a very mild correlation at 27 days. Otherwise, the residual from the daily data is close to random noise from a Student-t distribution with about 15 degrees of freedom. At the weekly resolution, there no serial dependence over horizons of 1–3 weeks, and only very mild dependence at four weeks. By implica- tion, the models are capturing nearly all the signal in the geomagnetic indexes, and the forecast error is close to random.

One issue that emerges in this context is whether the residual is random in the conventional sense, or may have an underlying fractal structure. In Consoloni et al. (2013) the geomagnetic indexes were found to show evidence of fractality at lower resolutions, which is consistent with the generating process, i.e., interaction between two or more stochastic series.

To evaluate this, the time series are tested using multiple scaling, a technique which can identify three parameters of interest, the order of integration, a coefficient characterizing the probability distribution, and the fractal codimension (Schertzer et al., 1997). The order of in- tegration is a measure of the memory of the series, or dependence be- tween distant time points. When the order of integration exceeds 0.5 but is less than 1, the series will trend irregularly, with repeated changes in direction. The coefficient characterizing the probability distribution is denoted α. When α = 2, the distribution is the Gaussian normal. When α = 1, the distribution is the Cauchy. When α lies in the range 1 < α < 2, the distribution has heavier tails than in the standard normal. The codimension, denoted C1, is the difference between the fractal dimension and the embedding dimension of the space (1 for a time series). As C1 approaches zero, the series becomes smoother. As C1 increases, the series is characterized by intermittent outliers.

Table 6 presents the estimated values. The coefficient of integration was estimated using the frequency domain regressions of Geweke and Porter-Hudak (1983), updated in Andrews and Guggenberger (2003). The order of integration is found to lie in the range 0.52 to 0.74, which is consistent with both nonstationarity and long memory (Baillie, 1996). The long memory traces back to cycles at multiple frequencies. The nonstationarity of the data is attributable to the up- ward trend in solar activity from the late nineteenth century until the 1960s, the Modern Maximum. All the series also exhibit thick tailed distributions, with α generally in the range of 1.4 to 1.6, and occa- sionally as low as <1.2.

The surprise is the low values for the codimension. Despite the evident volatility of the series, only irradiance shows strong evidence of an underlying fractal structure. This implies that the nonlinear varia- bility observed at high resolutions is primarily random noise. At first sight, this appears to contradict the finding that the geomagnetic in- dexes scale as multifractals at lower resolutions. On closer examination,

Table 4 Forecast accuracy, geomagnetic indexes.

Statistics are the mean absolute percent error.

Model Horizon (weeks)

1 2 3 4

Part 1: Aa Index Persistence 45.68 47.63 46.41 37.79 Regression on lags 35.73 36.02 36.64 36.77 Regression, lags and sunspots 35.71 36.52 37.35 37.42 Regression, lags and irradiance 35.92 36.23 37.09 37.22 Regression, state transition terms 35.71 36.07 36.65 36.81 Neural network 43.14 45.66 45.81 45.15 Frequency domain 35.42 35.36 35.76 34.55 Part 2: Am Index Persistence 45.54 47.46 46.77 38.61 Regression on lags 35.92 36.26 36.93 37.02 Regression, lags and sunspots 35.94 36.64 37.55 37.76 Regression, lags and irradiance 35.95 36.34 37.21 32.35 Regression, state transition terms 35.91 36.28 36.97 37.15 Neural network 43.75 45.91 45.67 45.39 Frequency domain 35.74 36.07 36.41 35.43

Table 5 Regressions for the geomagnetic indexes.

All variables are in natural logs.

RHS Variables Aa index (daily) Am index (daily) Aa Index (weekly) Am index (weekly)

Constant −2260 2.49 −1223 2.63 −1656 2.53 −1062 2.61 Irradiance 313 … 169 … 229 … 147 … Sunspots … 0.043 … 0.032 … 0.07 … 0.079 R-bar-square 0.364 0.365 0.353 0.354 0.211 0.234 0.154 0.178 Rho 0.586 0.581 0.587 0.583 0.405 0.359 0.355 0.324

All coefficients are significant at the 0.01% level or better. Rho is the coefficient of fractional differencing.

G. Reikard New Astronomy 62 (2018) 62–69

67

there is no real conflict between these results. The scaling test identifies a particular type of nonlinear variability, volatility arising from scaling symmetries over increasing horizons, typically arising from multi- plicative interactions among stochastic series. At lower resolutions, this is the case. However, at the highest resolutions, the series are to be dominated by embedding noise. If so, the large forecast errors are to be expected. At high frequencies, the signal-to-noise ratio is very low, and much of the variation is random.

6. Conclusions

The results are somewhat pessimistic. Over short horizons, on the order of days, and high resolutions, geomagnetic activity exhibits so much random variation that the forecast errors are extremely high. Over slightly longer horizons and slightly lower resolutions, there is some improvement from estimating in the frequency domain, but not a great deal. Including solar activity in the models does not yield any improvement in accuracy. The relationship between solar activity and the geomagnetic indexes, while statistically significant, is small. The state transition models also fail to generate any improvement, except at extremely short horizons. The key problem here is prediction of the state. At horizons beyond a day, the state becomes nearly impossible to predict.

The implication is that it may not be possible to achieve greater accuracy through any method. The degree of randomness in the data at high resolutions suggests that the errors cannot be reduced much fur- ther.

References

Active Cavity Radiometer Irradiance Monitor, 2014. www.acrim.com. Akaike, H., 1973. Information theory and the extension of the maximum likelihood

principle. In: Petrov, B.N., Csaki, F. (Eds.), Second International Symposium on Information Theory. Budapest, Akademiai Kiado, pp. 267–281.

Andrews, D.W.K., Guggenberger, P., 2003. A bias-reduced log-periodogram regression estimator for the long-memory parameter. Econometrica 71, 675–712 2003.

Baillie, R.T., 1996. Long memory processes and fractional integration in econometrics. J. Econometrics 73, 5–59.

Baker, D.N., McPherron, R.L., Cayton, T.E., Klebesadel, R.L., 1990. Linear prediction filter analysis of relativistic electron properties at 6.6 Re. J. Geophys. Res. 15133–15140.

Blanchard, G.T., McPherron, R.L., 1993. A bimodal representation of the response func- tion relating the solar wind electric field to the AL index. Adv. Space Res. 13, 71–74.

Blanchard, G.T., McPherron, R.L., 1995. Analysis of the linear response function relating AL to VBs for individual substorms. J. Geophys. Res. 100, 19155–19165.

Box, G.E.P., Jenkins, G.M., 1976. Time Series Analysis: Forecasting and Control. Holden- Day, San Francisco.

Buffett, B.A., 2000. Earth's core and the geodynamo. Science 288, 2007–2012. http://dx. doi.org/10.1126/science.288.5473.2007.

Buffett, B.A., 2010. Tidal dissipation and the strength of the Earth's internal magnetic field. Nature 468, 952–954. http://dx.doi.org/10.1038/nature09643.

Bunn, D.W., 2004. Modelling Prices in Competitive Electricity Markets. Wiley, New York. Cao, J.J., Cao, S.H., 2006. Study of forecasting solar irradiance using neural networks

with preprocessing sample data by wavelet analysis. Energy. 31, 3435–3445. http:// dx.doi.org/10.1016/j.energy.2006.04.001.

Chen, C., Liu, M.N., 1993. Forecasting time series with outliers. J. Forecast. 12, 13–35. Clette, F., Svalgaard, L., Vaquero, J.M., Cliver, E.W., 2014. Revisiting the sunspot number.

A 400-year perspective on the solar cycle. Space Sci. Rev. 186, 35–103. Codrescu, M.V., Fuller-Rowell, T.J., Minter, C.F., 2004. An ensemble-type Kalman filter

for neutral thermospheric composition during geomagnetic storms. Space Weather 2, S11002. http://dx.doi.org/10.1029/2004SW000088.

Consolini, G., De Marco, R., De Michelis, P., 2013. Intermittency and multifractional Brownian character of geomagnetic time series. Nonlinear Processes Geophys. 20, 455–466.

Feynman, J., Ruzmaikin, A., 2011. The sun's strange behavior: Maunder minimum or Gleissberg cycle? Sol. Phys. 272, 351–363.

Fox, R., Taqqu, M.S., 1986. Large sample properties of parameter estimates from strongly dependent stationary Gaussian time series. Ann. Statist. 14, 517–532.

Frohlich, C., 2009. Evidence of a long-term trend in total solar irradiance. Astron. Astrophys. 501 L27–U508.

Geweke, J., Porter-Hudak, S., 1983. The estimation and application of long-memory time series models. J. Time Ser. Anal. 4, 221–238.

Gholipour, A., Lucas, C., Araabi, B.N., Mirmomeni, M., Masoud Shafiee, M., 2007. Extracting the main patterns of natural time series for long-term neurofuzzy pre- diction. Neural Comput. Appl. 16, 383–393. http://dx.doi.org/10.1007/s00521-006- 0062-x.

Gleisner, H., Lundstedt, H., 1997. The response of the auroral electrojets to the solar wind modeled with neural networks. J. Geophys. Res. 102, 14269–14278.

Granger, C.W.J., 2008. Non-linear models: where do we go next – time varying parameter models. Studies in Nonlinear Dynamics and Econometrics 12. Article 1. http://www. bepress.com/snde/vol12/iss3/art1 2008.

Hamilton, J., 1989. A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica 57, 357–384.

Hamilton, J., 1990. Analysis of time series subject to changes in regime. J. Econometrics 45, 357–384.

Hansen, J., Willson, R., 1997. Forcings and chaos in interannual to decadal climate change. J. Geophys. Res. 102, 25679–25720.

Jordan, T.H., 1979. Structural geology of the Earth's interior. In: Proceedings of the National Academy of Sciences. 76. pp. 4192–4200. http://dx.doi.org/10.1073/pnas. 76.9.4192.

Kalhor, A., Araabi, B.N., Lucas, C., 2012. A new systematic design for habitually linear evolving TS fuzzy model. Expert Syst. Appl. 39, 1725–1736. http://dx.doi.org/10. 1016/j.eswa.2011.08.085.

Kalhor, A., Araabi, B.N., Lucas, C., 2011. Reducing the number of local linear models in neuro-fuzzy modeling: a split-and-merge clustering approach. Appl. Soft Comput. 11, 5582–5589. http://dx.doi.org/10.1016/j.asoc2011.05.001.

Kalman, R.E., 1960. A new approach to linear filtering and prediction problems. Trans. Am. Soc. Mech. Eng. J. Basic Eng. 83D, 35–45.

Koopmans, L.H., 1974. The Spectral Analysis of Time Series. Academic Press, New York. Krivova, N.A., Vieira, L.E.I., Solanki, S.K., 2010. Reconstruction of solar spectral irra-

diance since the Maunder minimum. J. Geophys. Res. 115 (A12). http://dx.doi.org/ 10.1029/2010JA015431.

Lovejoy, S., Schertzer, D., 2013. The Weather and Climate: Emergent Laws and Multifractal Cascades. Cambridge University Press, Cambridge.

Lundstedt, H., 1996. Solar origin of geomagnetic storms and predictions of storms with the use of neural networks. Surv. Geophys. 17, 561–573.

Lundstedt, H., 2005. Progress in space weather predictions and applications. Adv. Space. Res. 36, 2516–2523. http://dx.doi.org/10.1016/j.asr.2003.09.072.

Loskutov, A., Istomin, I.A., Kuzanyan, K.M., Kotlyarov, O.L., 2001. Testing and fore- casting the time series of the solar activity by singular spectrum analysis. Nonlinear Phenomena Complex Syst 4, 47–57.

Lu, W.Z., Wang, W.J., 2005. Potential assessment of the support vector machine method in forecasting ambient air pollutant trends. Chemosphere 59, 693–701. http://dx.doi. org/10.1016/j.chemosphere.2004.10.032.

Martin, J.D., Morton, Y.T., Zhou, Q., 2005. Neural network development for the fore- casting of upper atmosphere parameter distributions. Adv. Space Res. 36, 2480–2485. http://dx.doi.org/10.1016/j.asr.2003.12.15.

Mayaud, P.N., 1972. The Aa indices: a 100-year series characterizing the magnetic ac- tivity. J. Geophys. Res. 72, 6870–6874.

Mayaud, P.N., 1980. Derivation, meaning and use of geomagnetic indices. Geophysical Monograph 22. American Geophysical Union, Washington DC.

McPherron, R.L., 1993. Applications of artificial intelligence to space weather forecasting. In: Hruska, J., Shea, M.A., Smart, D.F., Heckman, G. (Eds.), Solar-Terrestrial Predictions – IV, Vol. 2, Magnetospheric Physics. U.S. Department of Commerce, NOAA, ERL, Boulder, Colorado, pp. 611–6263.

McPherron, R.L., 1995. Standard indices of geomagnetic activity. In: Kivelson, M.G., Russell, C.T. (Eds.), Introduction to Space Physics. Cambridge University Press, Cambridge, pp. 451–458.

Menvielle, M., Berthelier, A., 1991. The K-derived planetary indices – description and availability. Rev. Geophys. 29, 415–432.

Mirmomeni, M., Shafiee, M., Lucas, C., Araabi, B.N., 2006. Introducing a new learning method for fuzzy descriptor systems with the aid of spectral analysis to forecast solar activity. J. Atmosph. Solar Terrestrial Phys. 68, 2061–2074. http://dx.doi.org/10. 1016/j.astp2006.07.001.

Mirmomeni, M., Lucas, C., Araabi, B.N., Shafiee, M., 2007. Forecasting sunspot numbers with the aid of fuzzy descriptor models. Space Weather 5, S08004. http://dx.doi.org/ 10.1029/2006SW000289.

Mirmomeni, M., Lucas, C., Moshiri, B., Araabi, B.N., 2010. Introducing adaptive neuro- fuzzy modeling with online learning method for prediction of time-varying solar and

Table 6 Statistical properties of the data.

Time series Order of integration Alpha Fractal codimension

1) Daily resolution. Aa index 0.54 1.53 ( ± 0.06) 0.016 ( ± 0.12) Am index 0.52 1.39 ( ± 0.06) 0.008 ( ± 0.12) Irradiance 0.63 1.68 ( ± 0.05) 0.171 ( ± 0.10) Sunspots 0.71 1.17 ( ± 0.02) 0.034 ( ± 0.04)

2) Weekly resolution. Aa index 0.56 1.63 ( ± 0.04) 0.034 ( ± 0.08) Am index 0.54 1.54 ( ± 0.04) 0.028 ( ± 0.08) Irradiance 0.65 1.19 ( ± 0.07) 0.021 ( ± 0.14) Sunspots 0.74 1.41 ( ± 0.04) 0.061 ( ± 0.08)

The order of integration is the coefficient characterizing the memory of the series, or dependence between distant time points. The coefficient α characterizes the tail thickness of the probability distribution. The fractal codimension characterizes the intermittency of the series, but is also a measure of existence of scaling symmetries.

G. Reikard New Astronomy 62 (2018) 62–69

68

geomagnetic activity indices. Expert Syst. Appl. 37, 8267–8277. http://dx.doi.org/ 10.1016/j.eswa2010.05.059.

Mirmomeni, M., Lucas, C., Araabi, B.N., Moshiri, B., Bidar, M.R., 2011. Online multi-step ahead prediction of time-varying solar and geomagnetic activity indices via adaptive neurofuzzy modeling and recursive spectral analysis. Sol. Phys. 272, 189–213. http:// dx.doi.org/10.1007/s11207-011-9810-x.

National Geophysical Data Center, http://spidr.ngdc.noaa.gov/spidr, 2014. O'Brien, T.P., McPherron, R.L, 2000. A neural network model relating Delta H at a single

station to D/sub station. Adv. Space Res. 26, 209–212. http://dx.doi.org/10.1016/ S0273-1177(99)01051-0.

Orfila, A., Ballester, J.L., Oliver, R., Alvarez, A., Tintore, J., 2002. Forecasting the solar cycle with genetic algorithms. Astron. Astrophys. 386, 313–318.

Papitashvili, V.O., Clauer, C.R., Killeen, T.L., Belov, B.A., Golyshev, S.A., Levitin, A.E., 1998. Linear modeling of ionospheric electrodynamics from the IMF and solar wind data: application for space weather forecast. Adv. Space Res. 22, 113–116.

Papaleonidas, A., Iliadis, L., 2013. Neurocomputing techniques to dynamically forecast spatiotemporal air pollution data. Evolving Syst. 4, 221–233. http://dx.doi.org/10. 1007/s12530-013-9078-5.

Rangaranjan, G.K., Iyemori, T., 1997. Time variations of geomagnetic activity indices Kp and Ap: an update. Ann. Geophys. 15, 1271–1290.

Reikard, G., 2011. Forecasting space weather: can new econometric methods improve accuracy? Adv. Space Res. 47, 2073–2080.

Reikard, G., 2013. Combining frequency and time domain models to forecast space weather. Adv. Space Res. 52, 622–632.

Reikard, G., 2015. Forecasting geomagnetic activity at monthly and annual horizons: time series models. J. Atmos. Sol. Terr. Phys. 133, 111–120.

Rossi, B., Inoue, A., 2012. Out-of-sample forecast tests robust to the choice of window size. J. Business Econ. Statis. 30, 432–453.

Russell, C.T., McPherron, R.L., 1978. Semiannual variation of geomagnetic activity. J. Geophys. Res. 78, 92–108.

Scafetta, N., Willson, R.C., 2014. ACRIM total solar irradiance satellite composite vali- dation versus TSI proxy models. Astrophys. Space Sci. http://dx.doi.org/10.1007/ s10509-013-1775-9.

Schertzer, D., Lovejoy, S., Schmitt, F., Chigirinskaya, Y., Marsan, D., 1997. Multifractal

cascade dynamics and turbulent intermittency. Fractals 5, 427–471. Solanki, S.K., Krivova, N.A., 2011. Analyzing solar cycles. Science 334, 916–917. http://

dx.doi.org/10.1126/science.1212555. Solanki, S.K., Usoskin, I.G., Kromer, B., Schussler, M., Beer, J., 2004. An unusually active

sun during recent decades compared to the previous 11,000 years. Nature 431, 1084–1087.

Steinhilber, F., Beer, J., 2013. Prediction of solar activity for the next 500 years. Journal of Geophysical Research: Space Physics 118, 1861–1867. http://dx.doi.org/10.1002/ jgra.50210,2013.

Tulunay, Y., Sibeck, D.G., Senalp, E.T., Tulunay, E., 2005. Forecasting magnetopause crossing locations by using neural networks. Adv. Space Res. 36, 2378–2383. http:// dx.doi.org/10.1016/j.asr.2004.04.015.

Vandegriff, J., Wagstaff, K., Ho, G., Plauger, J., 2005. Forecasting space weather: pre- dicting interplanetary shocks using neural networks. Adv Space Res. 36, 2323–2327. http://dx.doi.org/10.1016/j.asr.2004.09.022.

Vieira, L.E.A., Solanki, S.K., Krivova, N.A., Usoskin, I.G., 2011. Evolution of solar irra- diance during the Holocene. Astron. Astrophys. 531, A6.

Vorotnikov, V.S., Smith, C.W., Hu, Q., Szabo, A., Skoug, R.M, Cohen, C.M.S., 2008. Automated shock detection and analysis algorithm for space weather application. Space Weather 6, S03002. http://dx.doi.org/10.1029/2007SW000358.

Wang, H.N., Cui, Y.M., Li, R., Zhang, L.Y., Han, H., 2008. Solar flare forecasting model supported with artificial neural network techniques. Adv. Space Res. 42, 1464–1468. http://dx.doi.org/10.1016/j.asr.2007.06.070.

Weiss, N., 2002. Dynamos in planets, stars and galaxies. Astron. Geophys. 43, 3.09–3.15. http://dx.doi.org/10.1046/j.1468-4004.2002.43309.x.

Willson, R.C., Mordvinov, A.V., 2003. Secular total solar irradiance trend during solar cycles 21–23. Geophys. Res. Lett. 30, 1199. http://dx.doi.org/10.1029/ 2002GL016038.

Wintoft, P., Lundstedt, H., 1998. Space weather modeling with intelligent hybrid systems: predicting the solar wind velocity. Adv. Space Res. 22, 59–62.

Wu, J.G., Lundstedt, H., 1997a. Geomagnetic storm predictions from solar wind data with the use of dynamic neural networks. J. Geophys. Res. 102, 14255–14268.

Wu, J.G., Lundstedt, H., 1997b. Neural network modeling of solar wind-magnetosphere interaction. J. Geophys. Res. 102, 14457–14466.

G. Reikard New Astronomy 62 (2018) 62–69

69