- Email: [email protected]
Solar activity forecast: Spectral analysis and neurofuzzy prediction
ARTICLE IN PRESS
Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603 www.elsevier.com/locate/jastp
Solar activity forecast: Spectral analysis and neurofuzzy prediction Ali Gholipoura,b,, Caro Lucasa,b, Babak N. Araabia,b, Masoud Shafieec a
Electrical and Computer Engineering Department, Control and Intelligent Processing Center of Excellence, University of Tehran, Tehran, Iran b School of Cognitive Sciences, Institute for studies in theoretical Physics and Mathematics (IPM), Tehran, Iran c Faculty of Electrical Engineering, Amir Kabir University of Technology, Tehran, Iran Received 22 October 2003; received in revised form 24 November 2004; accepted 7 December 2004
Abstract Active research in the last two decades indicates that the physical precursor and solar dynamo techniques are preferred as practical tools for long-term prediction of solar activity. But why should we omit more than 23 cycles of solar activity history, and just use empirical methods or simple autoregressive methods on the basis of observations for the latest eight cycles? In this article, a method based on spectral analysis and neurofuzzy modeling is proposed that is capable of issuing very accurate long-term prediction of sunspot number time series. A locally linear neurofuzzy model is optimized for each of the principal components obtained from singular spectrum analysis, and the multi-step predicted values are recombined to make the sunspot number time series. The proposed method is used for solar cycles 22 and 23 and the results are remarkably good in comparison to the predictions made by solar dynamo and precursor methods. An early prediction of the maximum smoothed international sunspot number for cycle 24 is 145 in 2011–2012. r 2005 Elsevier Ltd. All rights reserved. Keywords: Solar activity; Sunspot numbers; Forecasting; Singular spectrum analysis; Neurofuzzy modeling; Long-term prediction
1. Introduction Most of the space weather phenomena are influenced by variations in solar activity. During the years of solar maximum there are more solar flares causing significant increase in solar cosmic ray intensity. The high-energy particles disturb communication systems and affect the lifetime of satellites. Coronal mass ejections and solar Corresponding author. Electrical and Computer Engineer-
ing Department, Control and Intelligent Processing Center of Excellence, University of Tehran, Tehran, Iran. Fax: (+9821) 872 5029. E-mail addresses: [email protected], ali_gholipour@ yahoo.com (A. Gholipour), [email protected] (C. Lucas).
flares are the origin of shocks in solar wind and cause geomagnetic disturbances in the earth’s magnetosphere. The high rate of geomagnetic storms and sub-storms results in atmosphere heating and drag of Low Earth Orbit (LEO) satellites. Long-term solar activity forecasting is especially useful to space mission centres because the orbital trajectory parameters of satellites are greatly affected by the changing solar activity. There are also many arguments in recent years on correlations between solar activity and terrestrial climate (Laut, 2003; Oh et al., 2003; Tsiropoula, 2003). The level of solar activity is usually expressed by Zurich or International sunspot number ðRZ Þ or the 2800 MHz (10.7 cm) radio flux ðF10:7Þ index. The sunspot numbers are commonly used in numerical analysis because they are
1364-6826/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2004.12.001
ARTICLE IN PRESS 596
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
available for more than 23 cycles, while the F 10:7 index is preferred in space weather studies because of its higher signal-to-noise ratio and direct relationship to the solar EUV-UV emission. The correlation coefficient between the two indices is 0.987 and a simple least-squares estimation yields the following formula (Sofia et al., 1998): F 10:7 ¼ 58:7 þ 0:918RZ . Much effort has been dedicated to find physical evidence for the solar activity, to discover relationships between the level of activity and the effects of solar events on earth and environment, which is useful in reconstructing the sunspot number time series; and, interestingly, observations like the records of the 11 Be concentration in polar ice show that the period of high solar activity during the last 60 years is unique throughout the past 1150 years (Usokin et al., 2003). So, although the solar activity presents some clear periodicities, its prediction is quiet difficult; and despite the many statistical and experimental arguments of finding meaningful structures in the observations, it more and more shows long-term and large-scale changes. However, with the current accurate records of up to 300 years, seeking for longterm predictions of about 10 years is reasonable and seems feasible. Various numerical prediction techniques have been used for the sunspot number time series, e.g. Fourier analyses, curve fitting, artificial intelligence, and neural networks (Gholipour et al., 2003). These methods, although very accurate in short-term predictions, are not reliable in long term. The geomagnetic precursor techniques are based on statistical relationships between geomagnetic activity indices in the declining phase of one cycle (at the end of the cycle) and the activity level of the next solar maximum (Brown, 1992; Thompson, 1993; Joselyn, 1997). The first studies have used the geomagnetic aa index for prediction, while it is argued that the predictions on the basis of Ap index are more accurate (Sofia et al., 1998). All the precursor methods are investigated via the statistical tests for just eight solar cycles including eight numbers for solar maximums, and with such few data samples the best fitting is linear autoregressive estimation, and the reliability of the model is doubtful. Schatten et al. proposed a more acceptable method on the physical basis of solar dynamo theory and introduced the solar dynamo amplitude (SODA) index to predict the solar activity several years in advance (Schatten et al., 1978; Schatten and Sofia, 1987). This method has been used several times to forecast the last two cycles (Schatten and Pesnell, 1993; Schatten et al., 1996; Schatten, 2002), and seems to be more reliable than the common precursor methods, but has not been validated yet. The current forecasting approaches are developed in the direction of ignoring the long history of sunspot number time series
and trying to find simple empirical models on the basis of a few observations of the peak of the latest cycles and the geomagnetic indices. The timescale of these indices is much smaller than the variations in solar cycle, and their long-term trends cannot reveal the nonlinearities or complex periodicities of solar activity. In this article, a decomposition method based on singular spectrum analysis is used to make an intuitive nonlinear black box modeling technique applicable to long-term prediction of sunspot number time series. The singular spectrum analysis (SSA) (Vautard and Ghil, 1989; Vautard et al., 1992) performs a data adaptive filtering in the lag coordinate space of data and yields the principal components of the time series, which have narrow band frequency spectra and obvious temporal patterns. The principal components can be classified into linear or nonlinear trends, periodic or quasi-periodic patterns, and colored noise components, most of which are long-term predictable. Thus, after adaptive noise cancellation, one can design appropriate models for the components, and recombine the predicted components to obtain the original time series. The proposed method is successfully tested for the 22nd and 23rd solar cycles and the results are compared to the predictions by precursor and solar dynamo techniques. The mathematical description of the SSA algorithm is presented in Section 2; the locally linear neurofuzzy modeling technique along with the appropriate learning algorithm is considered in Section 3, and the results and their comparison are presented in Section 4. The last section contains some concluding remarks.
2. Spectral analysis SSA is defined as a new tool to extract information from short and noisy chaotic time series (Vautard et al., 1992). It relies on the Karhunen–Loeve decomposition of an estimate of covariance matrix based on M lagged copies of the time series. Thus, as the first step, the embedding procedure is applied to construct a sequence fX~ ðtÞg of M-dimensiona vectors from time series fX ðtÞ : t ¼ 1; . . . ; Ng: X~ ðtÞ ¼ ðX ðtÞ; X ðt þ 1Þ; . . . ; X ðt þ M 1ÞÞ, t ¼ 1; . . . ; N 0 ; N 0 ¼ N M þ 1.
ð1Þ
The N 0 M trajectory matrix (D) of the time series has the M-dimensional vectors as its columns, and is obviously a Hankel matrix (the elements on the diagonals i+j ¼ const are equal). In the second step, the M M covariance matrix C X is calculated and its eigenelements fðlk ; rk Þ : k ¼ 1; . . . ; Mg are determined by singular value decomposition (SVD). Each eigenvalue, lk ; estimates the partial variance in the direction of rk ; and the sum of all eigenvalues equals the total
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
597
Fig. 1. The first 12 principal components of sunspot number time series obtained from SSA: there are linear and nonlinear trends and multi-periodic components, which are the sources of quasi-periodic variations and a sense of chaotic behavior.
variance of the original time series. In the third step, the time series is projected onto each eigenvector and yields the corresponding principal components: Ak ðtÞ ¼
M X
X ðt þ j 1Þrk ðjÞ.
(2)
j¼1
Each of the principal components, being a nonlinear or linear trend or a periodic or quasi-periodic pattern, has narrow band frequency spectra and well-defined characteristics to be estimated. As the fourth step, the time series is reconstructed by combining the associated principal components: RK ðtÞ ¼
Ut 1 XX Ak ðt j þ 1Þrk ðjÞ. M t k2K j¼Lt
(3)
The normalization factor ðM t Þ; and the lower ðLt Þ and upper ðU t Þ bounds of the reconstruction procedure differ for the center and edges of the time series, and are defined by the following formula: ðM t ; Lt ; U t Þ 8 1 > > ; 1; t ; > > > t > > > > < 1 ; 1; M ; ¼ M > > > > > > 1 > > > : N t þ 1 ; t N þ M; M ;
For the sunspot number time series, M is chosen as 35 to cover at least three solar cycles. According to the singular spectrum logarithmic plot, the first twelve components are chosen as principal components, and the others are omitted to enhance signal-to-noise ratio. Fig. 1 shows the principal components; the first component is a long-term trend with a possible period of more than 100 years. The second and third components are the main periodic patterns with approximate period of eleven years and are the sources of quasi-periodic variations, and the other components can be either modeled by linear or nonlinear trends or periodic patterns. An efficient neurofuzzy modeling framework is introduced in the next section, which automatically constructs the best linear or nonlinear model for each of the components. After applying the neurofuzzy modeling and prediction, the extrapolated components are recombined via the fourth stage of SSA to form the long-term prediction of sunspot number time series.
3. Neurofuzzy modeling 1ptpM 1; 0
MptpN ;
ð4Þ
N 0 þ 1ptpN:
The fundamental approach with the locally linear neuro fuzzy (LLNF) model is dividing the input space into small linear subspaces with fuzzy validity functions (Nelles, 1999, 2001). Any produced linear part with its validity function can be described as a fuzzy neuron. Thus, the total model is a neurofuzzy network with one
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
598
hidden layer, and a linear neuron in the output layer, which simply calculates the weighted sum of the outputs of locally linear neurons: y^ i ¼ oi0 þ oi1 u1 þ oi2 u2 þ þ oip up , y^ ¼
M X
(5)
y^ i fi ðuÞ,
(6)
i¼1
where u ¼ ½u1 u2 up T is the model input, M is the number of locally linear neurons, and oij denotes the linear estimation parameters of the ith neuron. The validity functions are chosen as normalized Gaussians: m ðuÞ , fi ðuÞ ¼ PMi j¼1 mj ðuÞ ðup cip Þ2 1 ðu1 ci1 Þ2 þ þ mi ðuÞ ¼ exp 2 s2i1 s2ip 1 ðu1 ci1 Þ2 ¼ exp 2 s2i1 ! 1 ðup cip Þ2 . exp 2 s2ip
(7) !!
ð8Þ
The M:p parameters of the nonlinear hidden layer are the parameters of Gaussian validity functions: center ðcij Þ and standard deviation ðsij Þ: Optimization or learning methods are used to adjust the two sets of parameters, the rule consequent parameters of the locally linear models (oij ’s) and the rule premise parameters of validity functions (cij ’s and sij ’s). Global optimization of linear consequent parameters is simply obtained by the least-squares technique (Nelles, 2001). An incremental tree-based learning algorithm, e.g. locally linear model tree (LoLiMoT), is appropriate for tuning the rule premise parameters, i.e. determining the validation hypercube for each locally linear model. In each iteration the worst performing locally linear neuron is determined to be divided. All the possible divisions in the p-dimensional input space are checked and the best is performed. The splitting ratio can be simply adjusted as 12; which means that the locally linear neuron is divided into two equal halves. The fuzzy validity functions for the new structure are updated; their centers are the centers of the new hypercubes, and the standard deviations are usually set as 0.7. The algorithm is as follows: 1. The initial model: start with a single locally linear neuron, which is a globally optimal linear leastsquares estimation over the whole input space with F1 ðuÞ ¼ 1 and M ¼ 1: 2. Find the worst neuron: Calculate a local loss function, e.g. MSE for each of the i ¼ 1; . . . ; M locally linear neurons, and find the worst performing neuron.
3. Check all divisions: The worst neuron is considered for further refinement. The validation hypercube of this neuron is divided into two halves with an axis orthogonal split. Divisions in all dimensions are tried, and for each of the p divisions the following steps are carried out: (a) Construction of the multi-dimensional validity functions for both generated hypercubes. (b) Local estimation of the rule consequent parameters for both newly generated neurons. (c) Calculation of the total loss function or error index for the current overall model. 4. Validate the best division: The best of the p alternatives in Step 3 is selected. If it results in reduction of loss functions or error indices on training and validation data sets, the related validity functions and neurons are updated, the number of neurons is incremented M ¼ M þ 1; and the algorithm continues from Step 2, otherwise the learning algorithm is terminated.
This automatic learning algorithm provides the best linear or nonlinear model with maximum generalization, and performs well in prediction applications. Just one parameter should be defined before running the algorithm: the embedding dimension (number of regressors from time series as the inputs to the neurofuzzy model). The well-known Takens’ theorem gives a sufficient condition as a lower bound for the embedding dimension of chaotic systems (Takens, 1981). There are also many other approaches, e.g. Mees et al. (1987), Aleksic (1991), Kennel et al. (1992), Cao (1997) and Ataei et al. (2003). But, because most of the principal components of sunspot number time series show linear behavior, the Quenouille’s partial autocorrelation function test is optimal for most of them, and is a good estimation for the nonlinear ones. Note that because the data samples are very few, most of the advanced methods cannot provide a good estimation of the embedding dimension. The incremental learning algorithm is applied to the principal components, and one locally linear neurofuzzy model is automatically constructed for each of the principal components. Fig. 2 shows the block diagram of this algorithm. Optimizing the neurofuzzy models for each of the principal components is obtained using the proposed learning algorithm on separate training and validation sets. Like a cross validation procedure, the validation set is randomly chosen as a 10% part of the total data set. The trained locally linear models are used as predictor models for the corresponding components, and the long-term prediction of sunspot number time series is simply reconstructed from the predicted principal components via the fourth stage of SSA algorithm. The total model consists of relatively large number of parameters, but these parameters are
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
599
Fig. 2. The proposed method consists of four steps: (1) decomposition of the time series to nonlinear and periodic principal components via SSA stages 1 to 3, (2) locally linear modeling of the principal components, (3) multi-step prediction of the reconstructed principal components and (4) reconstructing the prediction of time series by combining the predicted principal components via SSA stage 4.
Table 1 Final configuration of the neurofuzzy models in predicting the solar cycle 22 PC# Regressors Neurons NMSE
1 20 1 0.005
2
3
4
5
6
5 3 0.048
4 5 0.024
5 1 0.019
6 1 0.020
5 3 0.005
7 15 1 0.028
8 12 2 0.023
9 5 1 0.019
10 8 5 0.005
11
12
14 1 0.005
15 1 0.005
Number of regressors, number of neurons, and mean square error (MSE) of the validation sets are presented for each of the principal components.
optimized independently for the principal components, and over-parameterization is avoided by considering the error indices on separate validation sets. Therefore, every box named LLNF in Fig. 2 is a complete trainable predictor model and is validated statistically through the fourth step of its learning procedure.
4. Solar activity forecasting The proposed method is used in long-term prediction of the last two solar cycles (cycles 22 and 23), where the results can be compared to actual values and the predictions made by precursor and solar dynamo techniques. The 22nd solar cycle, with an exceptionally large maximum, resulted in fault of many of the numerical and precursor methods, and was the starting point for the success of the solar dynamo technique; three years before the peak, Schatten et al. published a prediction of 210 25 for F 10:7 (160 25 for RZ ) (Schatten and Sofia, 1987), and the smoothed international sunspot number for the peak was actually 158.5. By the suggested method of this article, two predictions are made for the 22nd cycle; one is 158.0, 9 years before the maximum, and the other is 163.5, 5 years before the maximum. The predictor model parameters
are presented in Table 1, where the number of regressors is determined by Quenouille’s partial autocorrelation function test, and the number of neurons is optimized by the automatic tree-based learning algorithm. The mean-square error index of the validation sets for each of the principal components is also shown in this table. The predictions are shown in Figs. 3–5, and a comparison of the results is presented in Table 2. For the 23rd solar cycle, the precursor and solar dynamo techniques were used many times. The actual value of smoothed sunspot number for the maximum in 2000 was considerably smaller than the two latest peaks, and the early prediction by solar dynamo technique was absolutely wrong: it was 170 25 for RZ ; and the 23rd cycle was argued to be comparable to the 21st and 22nd cycles (Schatten and Pesnell, 1993). Even the late prediction by a combination of precursor methods in 1997 (160 30 for RZ ) did not include the actual value in one sigma interval (Joselyn, 1997). In 1998 the following precursor formulas were used to predict the peak of the cycle (Sofia et al., 1998): RZ ¼ 8:45ð1:34Þaa þ 13:53ð2:12Þ, RZ ¼ 18:73ð4:65ÞAp 25:70ð2:37Þ.
ARTICLE IN PRESS 600
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
Fig. 3. Long-term prediction of solar cycle 22 by the proposed method (20 solar cycles are used for training); (A) dotted: the observed values, solid line from 1920 to 1981: the reconstructed and one step ahead predicted values, and thick line from 1981 to 1995: multi-step predicted values.
Fig. 4. Long-term prediction of solar cycle 22 by the proposed method (only 8 solar cycles are used for training); (A) dotted: the observed values, solid line from 1920 to 1985: the reconstructed and one step ahead predicted values, and thick line from 1985 to 1994: multi-step predicted values.
Fig. 5. Long-term prediction of solar cycle 22 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 1985: the reconstructed and one step ahead predicted values, and thick line from 1985 to 1994: multi-step predicted values. Table 2 Several methods in long-term prediction of solar cycle 22 via the smoothed peak of sunspot number time series
The proposed method The proposed method Precursor Solar dynamo
Time of prediction
Predicted value for maximum
Reference
1981 1985 1987 1987
158 163.5 170 25 160 25
— — Sofia et al. (1998) Schatten and Sofia (1987)
The observed value is 158.5 in 1990.
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
601
Table 3 Several methods in long-term prediction of solar cycle 23 via the smoothed peak of sunspot number time series
The proposed method Solar dynamo Precursor Solar dynamo The proposed method Precursor (aa index) Precursor (Ap index) Dynamo-based
Time of prediction
Predicted value of maximum
Reference
1991 1993 1996 1996 1996 1998 1998 1998
141 170 25 160 30 138 30 132 148 25 130 40 143 30
— Schatten and Pesnell (1993) Joselyn (1997) Schatten et al. (1996) — Sofia et al. (1998) Sofia et al. (1998) Sofia et al. (1998)
the observed value is 120 in 2000.
Fig. 6. Long-term prediction of solar cycle 23 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 1995: the reconstructed and one step ahead predicted values, and thick line from 1996 to 2003: multi-step predicted values.
Table 4 Error indices of the predictions made by the proposed method; root-mean-square error (RMSE), and maximum absolute error of the long-term predicted values Predicted solar cycle
Prediction horizon
RMSE
Maximum absolute error
22 22 23 23
5 9 5 9
6.6250 7.4724 7.1851 7.8059
17.5 19.8 22.3 21.5
depicted in Fig. 7. Note that an early prediction of this cycle by the solar dynamo technique is 175 50 in F10:7; which is very close to the peak of the 23rd cycle (Schatten, 2002, 2003).
5. Discussion and conclusions years years years years
The next predictions by the solar dynamo method were better in accuracy (Schatten et al., 1996; Sofia et al., 1998). The history of these predictions, and two predictions made by the suggested method of this paper is presented in Table 3. Fig. 6 shows the multi-step prediction of cycle 23 by the suggested method Table 4. The proposed method is used to predict the next cycle (cycle 24); the peak of the next solar maximum is expected to be near 145 for RZ : This means that it will be about 20% larger than the 23rd cycle and 10% smaller than the 22nd cycle. The shape of the next cycle is
The current practical methods in long-term prediction of solar activity, the empirical solar dynamo and precursor techniques, are investigated via the observations of the latest eight solar cycles because the measurements of geomagnetic activity and solar magnetic field before that time were not available. The best statistical methods for these few observations yield linear autoregressive models, and the physical solar dynamo method is not validated yet. This paper goes back to the numerical techniques to extract the most available information from the history of sunspot number time series. By applying a nonlinear spectral analysis tool, the singular spectrum analysis, the principal patterns of the sunspot number time series are extracted. An optimal linear or locally linear neurofuzzy model is constructed for each of the components by an incremental learning algorithm and
ARTICLE IN PRESS 602
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
Fig. 7. Long-term prediction of solar cycle 24 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 2001: the reconstructed and one step ahead predicted values, and thick line from 2001 to 2012: multi-step predicted values.
is used in multi-step prediction. The original time series is achieved by recombining the extrapolated components. The long-term predictions by this method are superior in comparison to predictions made by the precursor and solar dynamo techniques for the last two solar cycles, and the peak of the next solar maximum is predicted to be near 145 in 2011–2012.
Acknowledgements The authors would like to thank the anonymous reviewers for their helpful comments; especially they wish to appreciate the first reviewer for a complementary discussion on the importance of the subject.
References Aleksic, Z., 1991. Estimating the embedding dimension. Physica D 52, 362–368. Ataei, M., Khaki-Sedigh, A., Lohmann, B., Lucas, C., 2003. Determination of embedding dimension using multiple time series based on singular value decomposition. In: Troch, I., Breitenecker, F. (Eds.), Proceedings of the Fourth International Symposium on Mathematical Modeling. Vienna, Austria, pp. 190–196. Brown, G.M., 1992. The peak of solar cycle 22: predictions in retrospect. Annales Geophysicae 10, 453–470. Cao, L., 1997. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D 110, 43–50. Gholipour, A., Abbaspour, A., Araabi, B.N., Lucas, C., 2003. Enhancements in the prediction of solar activity by locally linear model tree. Proceedings of MIC2003: 22nd International Conference on Modeling, Identification and Control, Innsbruck, Austria, pp. 158–161. Joselyn, J.A., 1997. Panel achieves consensus prediction of solar cycle. EOS Transactions American Geophysical Union 78, 211–212. Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining embedding dimension for phase space reconstruction
using a geometrical construction. Physical Review A 45 (6), 3403–3411. Laut, P., 2003. Solar activity and terrestrial climate: an analysis of some purported correlations. Journal of Atmospheric and Solar-Terrestrial Physics 65 (7), 801–812. Mees, A.I., Rapp, P.E., Jenning, L.S., 1987. Singular value decomposition and embedding dimension. Physical Review A 36 (1), 340–346. Nelles, O., 1999. Nonlinear system identification with local linear neuro-fuzzy models. Ph.D. Thesis, TU Darmstadt, Shaker Verlag, Aachen. Nelles, O., 2001. Nonlinear system identification. Springer, Berlin. Oh, H., Ammann, C.M., Naveau, P., Nychka, D., OttoBliesner, B.L., 2003. Multi-resolution time series analysis applied to solar irradiance and climate reconstructions. Journal of Atmospheric and Solar-Terrestrial Physics 65 (2), 191–201. Schatten, K., 2002. Solar activity prediction: timing predictors and cycle 24. Journal of Geophysical Research 107 (A11), 1377 doi:10.1029/2002JA009404. Schatten, K., 2003. Correction to Solar activity prediction: timing predictors and cycle 24 by Kenneth Schatten. Journal of Geophysical Research 108 (A3), 1100 doi:10.1029/ 2003JA009843. Schatten, K.H., Pesnell, W.D., 1993. An early solar dynamo prediction: cycle 23–cycle 22. Geophysical Research Letters 20, 2257–2278. Schatten, K.H., Sofia, S., 1987. Forecast of an exceptionally large even numbered solar cycle. Geophysical Research Letters 14, 632. Schatten, K.H., Scherrer, P.H., Svalgaard, L., Wilcox, J.M., 1978. Using Dynamo theory to predict the sunspot number during solar cycle 21. Geophysical Research letters 5, 411. Schatten, K.H., Myers, D.J., Sofia, S., 1996. Solar activity forecast for solar cycle 23. Geophysical Research Letters 23 (6), 605–608. Sofia, S., Fox, P., Schatten, K.H., 1998. Forecast update for activity cycle 23 from a Dynamo-based method. Geophysical Research Letters 25 (22), 4149–4152. Takens, F., 1981. Detecting strange attractors in turbulence. In: Rand, D.A., Young, L.S. (Eds.), Lecture Notes in Mathematics, Springer, Berlin, pp. 366–381.
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603 Thompson, R.J., 1993. A technique for predicting the amplitude of solar cycle. Solar Physics 148, 383. Tsiropoula, G., 2003. Signatures of solar activity variability in meteorological parameters. Journal of Atmospheric and Solar-Terrestrial Physics 65 (4), 469–482. Usokin, I.G., Solanki, S.K., Schussler, M., Mursula, K., Alanko, K., 2003. Millennium-scale sunspot number reconstruction: evidence for an unusually active Sun since
603
the 1940s. Physical Review Letters, doi: 10.1103/PhysRevLett.91.211101. Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynamics with applications to paleoclimatic time series. Physica D 35, 395–424. Vautard, R., Yiou, P., Ghil, M., 1992. Singular spectrum analysis: a toolkit for short noisy chaotic signals. Physica D 58, 95–126.
Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603 www.elsevier.com/locate/jastp
Solar activity forecast: Spectral analysis and neurofuzzy prediction Ali Gholipoura,b,, Caro Lucasa,b, Babak N. Araabia,b, Masoud Shafieec a
Electrical and Computer Engineering Department, Control and Intelligent Processing Center of Excellence, University of Tehran, Tehran, Iran b School of Cognitive Sciences, Institute for studies in theoretical Physics and Mathematics (IPM), Tehran, Iran c Faculty of Electrical Engineering, Amir Kabir University of Technology, Tehran, Iran Received 22 October 2003; received in revised form 24 November 2004; accepted 7 December 2004
Abstract Active research in the last two decades indicates that the physical precursor and solar dynamo techniques are preferred as practical tools for long-term prediction of solar activity. But why should we omit more than 23 cycles of solar activity history, and just use empirical methods or simple autoregressive methods on the basis of observations for the latest eight cycles? In this article, a method based on spectral analysis and neurofuzzy modeling is proposed that is capable of issuing very accurate long-term prediction of sunspot number time series. A locally linear neurofuzzy model is optimized for each of the principal components obtained from singular spectrum analysis, and the multi-step predicted values are recombined to make the sunspot number time series. The proposed method is used for solar cycles 22 and 23 and the results are remarkably good in comparison to the predictions made by solar dynamo and precursor methods. An early prediction of the maximum smoothed international sunspot number for cycle 24 is 145 in 2011–2012. r 2005 Elsevier Ltd. All rights reserved. Keywords: Solar activity; Sunspot numbers; Forecasting; Singular spectrum analysis; Neurofuzzy modeling; Long-term prediction
1. Introduction Most of the space weather phenomena are influenced by variations in solar activity. During the years of solar maximum there are more solar flares causing significant increase in solar cosmic ray intensity. The high-energy particles disturb communication systems and affect the lifetime of satellites. Coronal mass ejections and solar Corresponding author. Electrical and Computer Engineer-
ing Department, Control and Intelligent Processing Center of Excellence, University of Tehran, Tehran, Iran. Fax: (+9821) 872 5029. E-mail addresses: [email protected], ali_gholipour@ yahoo.com (A. Gholipour), [email protected] (C. Lucas).
flares are the origin of shocks in solar wind and cause geomagnetic disturbances in the earth’s magnetosphere. The high rate of geomagnetic storms and sub-storms results in atmosphere heating and drag of Low Earth Orbit (LEO) satellites. Long-term solar activity forecasting is especially useful to space mission centres because the orbital trajectory parameters of satellites are greatly affected by the changing solar activity. There are also many arguments in recent years on correlations between solar activity and terrestrial climate (Laut, 2003; Oh et al., 2003; Tsiropoula, 2003). The level of solar activity is usually expressed by Zurich or International sunspot number ðRZ Þ or the 2800 MHz (10.7 cm) radio flux ðF10:7Þ index. The sunspot numbers are commonly used in numerical analysis because they are
1364-6826/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2004.12.001
ARTICLE IN PRESS 596
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
available for more than 23 cycles, while the F 10:7 index is preferred in space weather studies because of its higher signal-to-noise ratio and direct relationship to the solar EUV-UV emission. The correlation coefficient between the two indices is 0.987 and a simple least-squares estimation yields the following formula (Sofia et al., 1998): F 10:7 ¼ 58:7 þ 0:918RZ . Much effort has been dedicated to find physical evidence for the solar activity, to discover relationships between the level of activity and the effects of solar events on earth and environment, which is useful in reconstructing the sunspot number time series; and, interestingly, observations like the records of the 11 Be concentration in polar ice show that the period of high solar activity during the last 60 years is unique throughout the past 1150 years (Usokin et al., 2003). So, although the solar activity presents some clear periodicities, its prediction is quiet difficult; and despite the many statistical and experimental arguments of finding meaningful structures in the observations, it more and more shows long-term and large-scale changes. However, with the current accurate records of up to 300 years, seeking for longterm predictions of about 10 years is reasonable and seems feasible. Various numerical prediction techniques have been used for the sunspot number time series, e.g. Fourier analyses, curve fitting, artificial intelligence, and neural networks (Gholipour et al., 2003). These methods, although very accurate in short-term predictions, are not reliable in long term. The geomagnetic precursor techniques are based on statistical relationships between geomagnetic activity indices in the declining phase of one cycle (at the end of the cycle) and the activity level of the next solar maximum (Brown, 1992; Thompson, 1993; Joselyn, 1997). The first studies have used the geomagnetic aa index for prediction, while it is argued that the predictions on the basis of Ap index are more accurate (Sofia et al., 1998). All the precursor methods are investigated via the statistical tests for just eight solar cycles including eight numbers for solar maximums, and with such few data samples the best fitting is linear autoregressive estimation, and the reliability of the model is doubtful. Schatten et al. proposed a more acceptable method on the physical basis of solar dynamo theory and introduced the solar dynamo amplitude (SODA) index to predict the solar activity several years in advance (Schatten et al., 1978; Schatten and Sofia, 1987). This method has been used several times to forecast the last two cycles (Schatten and Pesnell, 1993; Schatten et al., 1996; Schatten, 2002), and seems to be more reliable than the common precursor methods, but has not been validated yet. The current forecasting approaches are developed in the direction of ignoring the long history of sunspot number time series
and trying to find simple empirical models on the basis of a few observations of the peak of the latest cycles and the geomagnetic indices. The timescale of these indices is much smaller than the variations in solar cycle, and their long-term trends cannot reveal the nonlinearities or complex periodicities of solar activity. In this article, a decomposition method based on singular spectrum analysis is used to make an intuitive nonlinear black box modeling technique applicable to long-term prediction of sunspot number time series. The singular spectrum analysis (SSA) (Vautard and Ghil, 1989; Vautard et al., 1992) performs a data adaptive filtering in the lag coordinate space of data and yields the principal components of the time series, which have narrow band frequency spectra and obvious temporal patterns. The principal components can be classified into linear or nonlinear trends, periodic or quasi-periodic patterns, and colored noise components, most of which are long-term predictable. Thus, after adaptive noise cancellation, one can design appropriate models for the components, and recombine the predicted components to obtain the original time series. The proposed method is successfully tested for the 22nd and 23rd solar cycles and the results are compared to the predictions by precursor and solar dynamo techniques. The mathematical description of the SSA algorithm is presented in Section 2; the locally linear neurofuzzy modeling technique along with the appropriate learning algorithm is considered in Section 3, and the results and their comparison are presented in Section 4. The last section contains some concluding remarks.
2. Spectral analysis SSA is defined as a new tool to extract information from short and noisy chaotic time series (Vautard et al., 1992). It relies on the Karhunen–Loeve decomposition of an estimate of covariance matrix based on M lagged copies of the time series. Thus, as the first step, the embedding procedure is applied to construct a sequence fX~ ðtÞg of M-dimensiona vectors from time series fX ðtÞ : t ¼ 1; . . . ; Ng: X~ ðtÞ ¼ ðX ðtÞ; X ðt þ 1Þ; . . . ; X ðt þ M 1ÞÞ, t ¼ 1; . . . ; N 0 ; N 0 ¼ N M þ 1.
ð1Þ
The N 0 M trajectory matrix (D) of the time series has the M-dimensional vectors as its columns, and is obviously a Hankel matrix (the elements on the diagonals i+j ¼ const are equal). In the second step, the M M covariance matrix C X is calculated and its eigenelements fðlk ; rk Þ : k ¼ 1; . . . ; Mg are determined by singular value decomposition (SVD). Each eigenvalue, lk ; estimates the partial variance in the direction of rk ; and the sum of all eigenvalues equals the total
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
597
Fig. 1. The first 12 principal components of sunspot number time series obtained from SSA: there are linear and nonlinear trends and multi-periodic components, which are the sources of quasi-periodic variations and a sense of chaotic behavior.
variance of the original time series. In the third step, the time series is projected onto each eigenvector and yields the corresponding principal components: Ak ðtÞ ¼
M X
X ðt þ j 1Þrk ðjÞ.
(2)
j¼1
Each of the principal components, being a nonlinear or linear trend or a periodic or quasi-periodic pattern, has narrow band frequency spectra and well-defined characteristics to be estimated. As the fourth step, the time series is reconstructed by combining the associated principal components: RK ðtÞ ¼
Ut 1 XX Ak ðt j þ 1Þrk ðjÞ. M t k2K j¼Lt
(3)
The normalization factor ðM t Þ; and the lower ðLt Þ and upper ðU t Þ bounds of the reconstruction procedure differ for the center and edges of the time series, and are defined by the following formula: ðM t ; Lt ; U t Þ 8 1 > > ; 1; t ; > > > t > > > > < 1 ; 1; M ; ¼ M > > > > > > 1 > > > : N t þ 1 ; t N þ M; M ;
For the sunspot number time series, M is chosen as 35 to cover at least three solar cycles. According to the singular spectrum logarithmic plot, the first twelve components are chosen as principal components, and the others are omitted to enhance signal-to-noise ratio. Fig. 1 shows the principal components; the first component is a long-term trend with a possible period of more than 100 years. The second and third components are the main periodic patterns with approximate period of eleven years and are the sources of quasi-periodic variations, and the other components can be either modeled by linear or nonlinear trends or periodic patterns. An efficient neurofuzzy modeling framework is introduced in the next section, which automatically constructs the best linear or nonlinear model for each of the components. After applying the neurofuzzy modeling and prediction, the extrapolated components are recombined via the fourth stage of SSA to form the long-term prediction of sunspot number time series.
3. Neurofuzzy modeling 1ptpM 1; 0
MptpN ;
ð4Þ
N 0 þ 1ptpN:
The fundamental approach with the locally linear neuro fuzzy (LLNF) model is dividing the input space into small linear subspaces with fuzzy validity functions (Nelles, 1999, 2001). Any produced linear part with its validity function can be described as a fuzzy neuron. Thus, the total model is a neurofuzzy network with one
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
598
hidden layer, and a linear neuron in the output layer, which simply calculates the weighted sum of the outputs of locally linear neurons: y^ i ¼ oi0 þ oi1 u1 þ oi2 u2 þ þ oip up , y^ ¼
M X
(5)
y^ i fi ðuÞ,
(6)
i¼1
where u ¼ ½u1 u2 up T is the model input, M is the number of locally linear neurons, and oij denotes the linear estimation parameters of the ith neuron. The validity functions are chosen as normalized Gaussians: m ðuÞ , fi ðuÞ ¼ PMi j¼1 mj ðuÞ ðup cip Þ2 1 ðu1 ci1 Þ2 þ þ mi ðuÞ ¼ exp 2 s2i1 s2ip 1 ðu1 ci1 Þ2 ¼ exp 2 s2i1 ! 1 ðup cip Þ2 . exp 2 s2ip
(7) !!
ð8Þ
The M:p parameters of the nonlinear hidden layer are the parameters of Gaussian validity functions: center ðcij Þ and standard deviation ðsij Þ: Optimization or learning methods are used to adjust the two sets of parameters, the rule consequent parameters of the locally linear models (oij ’s) and the rule premise parameters of validity functions (cij ’s and sij ’s). Global optimization of linear consequent parameters is simply obtained by the least-squares technique (Nelles, 2001). An incremental tree-based learning algorithm, e.g. locally linear model tree (LoLiMoT), is appropriate for tuning the rule premise parameters, i.e. determining the validation hypercube for each locally linear model. In each iteration the worst performing locally linear neuron is determined to be divided. All the possible divisions in the p-dimensional input space are checked and the best is performed. The splitting ratio can be simply adjusted as 12; which means that the locally linear neuron is divided into two equal halves. The fuzzy validity functions for the new structure are updated; their centers are the centers of the new hypercubes, and the standard deviations are usually set as 0.7. The algorithm is as follows: 1. The initial model: start with a single locally linear neuron, which is a globally optimal linear leastsquares estimation over the whole input space with F1 ðuÞ ¼ 1 and M ¼ 1: 2. Find the worst neuron: Calculate a local loss function, e.g. MSE for each of the i ¼ 1; . . . ; M locally linear neurons, and find the worst performing neuron.
3. Check all divisions: The worst neuron is considered for further refinement. The validation hypercube of this neuron is divided into two halves with an axis orthogonal split. Divisions in all dimensions are tried, and for each of the p divisions the following steps are carried out: (a) Construction of the multi-dimensional validity functions for both generated hypercubes. (b) Local estimation of the rule consequent parameters for both newly generated neurons. (c) Calculation of the total loss function or error index for the current overall model. 4. Validate the best division: The best of the p alternatives in Step 3 is selected. If it results in reduction of loss functions or error indices on training and validation data sets, the related validity functions and neurons are updated, the number of neurons is incremented M ¼ M þ 1; and the algorithm continues from Step 2, otherwise the learning algorithm is terminated.
This automatic learning algorithm provides the best linear or nonlinear model with maximum generalization, and performs well in prediction applications. Just one parameter should be defined before running the algorithm: the embedding dimension (number of regressors from time series as the inputs to the neurofuzzy model). The well-known Takens’ theorem gives a sufficient condition as a lower bound for the embedding dimension of chaotic systems (Takens, 1981). There are also many other approaches, e.g. Mees et al. (1987), Aleksic (1991), Kennel et al. (1992), Cao (1997) and Ataei et al. (2003). But, because most of the principal components of sunspot number time series show linear behavior, the Quenouille’s partial autocorrelation function test is optimal for most of them, and is a good estimation for the nonlinear ones. Note that because the data samples are very few, most of the advanced methods cannot provide a good estimation of the embedding dimension. The incremental learning algorithm is applied to the principal components, and one locally linear neurofuzzy model is automatically constructed for each of the principal components. Fig. 2 shows the block diagram of this algorithm. Optimizing the neurofuzzy models for each of the principal components is obtained using the proposed learning algorithm on separate training and validation sets. Like a cross validation procedure, the validation set is randomly chosen as a 10% part of the total data set. The trained locally linear models are used as predictor models for the corresponding components, and the long-term prediction of sunspot number time series is simply reconstructed from the predicted principal components via the fourth stage of SSA algorithm. The total model consists of relatively large number of parameters, but these parameters are
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
599
Fig. 2. The proposed method consists of four steps: (1) decomposition of the time series to nonlinear and periodic principal components via SSA stages 1 to 3, (2) locally linear modeling of the principal components, (3) multi-step prediction of the reconstructed principal components and (4) reconstructing the prediction of time series by combining the predicted principal components via SSA stage 4.
Table 1 Final configuration of the neurofuzzy models in predicting the solar cycle 22 PC# Regressors Neurons NMSE
1 20 1 0.005
2
3
4
5
6
5 3 0.048
4 5 0.024
5 1 0.019
6 1 0.020
5 3 0.005
7 15 1 0.028
8 12 2 0.023
9 5 1 0.019
10 8 5 0.005
11
12
14 1 0.005
15 1 0.005
Number of regressors, number of neurons, and mean square error (MSE) of the validation sets are presented for each of the principal components.
optimized independently for the principal components, and over-parameterization is avoided by considering the error indices on separate validation sets. Therefore, every box named LLNF in Fig. 2 is a complete trainable predictor model and is validated statistically through the fourth step of its learning procedure.
4. Solar activity forecasting The proposed method is used in long-term prediction of the last two solar cycles (cycles 22 and 23), where the results can be compared to actual values and the predictions made by precursor and solar dynamo techniques. The 22nd solar cycle, with an exceptionally large maximum, resulted in fault of many of the numerical and precursor methods, and was the starting point for the success of the solar dynamo technique; three years before the peak, Schatten et al. published a prediction of 210 25 for F 10:7 (160 25 for RZ ) (Schatten and Sofia, 1987), and the smoothed international sunspot number for the peak was actually 158.5. By the suggested method of this article, two predictions are made for the 22nd cycle; one is 158.0, 9 years before the maximum, and the other is 163.5, 5 years before the maximum. The predictor model parameters
are presented in Table 1, where the number of regressors is determined by Quenouille’s partial autocorrelation function test, and the number of neurons is optimized by the automatic tree-based learning algorithm. The mean-square error index of the validation sets for each of the principal components is also shown in this table. The predictions are shown in Figs. 3–5, and a comparison of the results is presented in Table 2. For the 23rd solar cycle, the precursor and solar dynamo techniques were used many times. The actual value of smoothed sunspot number for the maximum in 2000 was considerably smaller than the two latest peaks, and the early prediction by solar dynamo technique was absolutely wrong: it was 170 25 for RZ ; and the 23rd cycle was argued to be comparable to the 21st and 22nd cycles (Schatten and Pesnell, 1993). Even the late prediction by a combination of precursor methods in 1997 (160 30 for RZ ) did not include the actual value in one sigma interval (Joselyn, 1997). In 1998 the following precursor formulas were used to predict the peak of the cycle (Sofia et al., 1998): RZ ¼ 8:45ð1:34Þaa þ 13:53ð2:12Þ, RZ ¼ 18:73ð4:65ÞAp 25:70ð2:37Þ.
ARTICLE IN PRESS 600
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
Fig. 3. Long-term prediction of solar cycle 22 by the proposed method (20 solar cycles are used for training); (A) dotted: the observed values, solid line from 1920 to 1981: the reconstructed and one step ahead predicted values, and thick line from 1981 to 1995: multi-step predicted values.
Fig. 4. Long-term prediction of solar cycle 22 by the proposed method (only 8 solar cycles are used for training); (A) dotted: the observed values, solid line from 1920 to 1985: the reconstructed and one step ahead predicted values, and thick line from 1985 to 1994: multi-step predicted values.
Fig. 5. Long-term prediction of solar cycle 22 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 1985: the reconstructed and one step ahead predicted values, and thick line from 1985 to 1994: multi-step predicted values. Table 2 Several methods in long-term prediction of solar cycle 22 via the smoothed peak of sunspot number time series
The proposed method The proposed method Precursor Solar dynamo
Time of prediction
Predicted value for maximum
Reference
1981 1985 1987 1987
158 163.5 170 25 160 25
— — Sofia et al. (1998) Schatten and Sofia (1987)
The observed value is 158.5 in 1990.
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
601
Table 3 Several methods in long-term prediction of solar cycle 23 via the smoothed peak of sunspot number time series
The proposed method Solar dynamo Precursor Solar dynamo The proposed method Precursor (aa index) Precursor (Ap index) Dynamo-based
Time of prediction
Predicted value of maximum
Reference
1991 1993 1996 1996 1996 1998 1998 1998
141 170 25 160 30 138 30 132 148 25 130 40 143 30
— Schatten and Pesnell (1993) Joselyn (1997) Schatten et al. (1996) — Sofia et al. (1998) Sofia et al. (1998) Sofia et al. (1998)
the observed value is 120 in 2000.
Fig. 6. Long-term prediction of solar cycle 23 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 1995: the reconstructed and one step ahead predicted values, and thick line from 1996 to 2003: multi-step predicted values.
Table 4 Error indices of the predictions made by the proposed method; root-mean-square error (RMSE), and maximum absolute error of the long-term predicted values Predicted solar cycle
Prediction horizon
RMSE
Maximum absolute error
22 22 23 23
5 9 5 9
6.6250 7.4724 7.1851 7.8059
17.5 19.8 22.3 21.5
depicted in Fig. 7. Note that an early prediction of this cycle by the solar dynamo technique is 175 50 in F10:7; which is very close to the peak of the 23rd cycle (Schatten, 2002, 2003).
5. Discussion and conclusions years years years years
The next predictions by the solar dynamo method were better in accuracy (Schatten et al., 1996; Sofia et al., 1998). The history of these predictions, and two predictions made by the suggested method of this paper is presented in Table 3. Fig. 6 shows the multi-step prediction of cycle 23 by the suggested method Table 4. The proposed method is used to predict the next cycle (cycle 24); the peak of the next solar maximum is expected to be near 145 for RZ : This means that it will be about 20% larger than the 23rd cycle and 10% smaller than the 22nd cycle. The shape of the next cycle is
The current practical methods in long-term prediction of solar activity, the empirical solar dynamo and precursor techniques, are investigated via the observations of the latest eight solar cycles because the measurements of geomagnetic activity and solar magnetic field before that time were not available. The best statistical methods for these few observations yield linear autoregressive models, and the physical solar dynamo method is not validated yet. This paper goes back to the numerical techniques to extract the most available information from the history of sunspot number time series. By applying a nonlinear spectral analysis tool, the singular spectrum analysis, the principal patterns of the sunspot number time series are extracted. An optimal linear or locally linear neurofuzzy model is constructed for each of the components by an incremental learning algorithm and
ARTICLE IN PRESS 602
A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603
Fig. 7. Long-term prediction of solar cycle 24 by the proposed method; (A) dotted: the observed values, solid line from 1920 to 2001: the reconstructed and one step ahead predicted values, and thick line from 2001 to 2012: multi-step predicted values.
is used in multi-step prediction. The original time series is achieved by recombining the extrapolated components. The long-term predictions by this method are superior in comparison to predictions made by the precursor and solar dynamo techniques for the last two solar cycles, and the peak of the next solar maximum is predicted to be near 145 in 2011–2012.
Acknowledgements The authors would like to thank the anonymous reviewers for their helpful comments; especially they wish to appreciate the first reviewer for a complementary discussion on the importance of the subject.
References Aleksic, Z., 1991. Estimating the embedding dimension. Physica D 52, 362–368. Ataei, M., Khaki-Sedigh, A., Lohmann, B., Lucas, C., 2003. Determination of embedding dimension using multiple time series based on singular value decomposition. In: Troch, I., Breitenecker, F. (Eds.), Proceedings of the Fourth International Symposium on Mathematical Modeling. Vienna, Austria, pp. 190–196. Brown, G.M., 1992. The peak of solar cycle 22: predictions in retrospect. Annales Geophysicae 10, 453–470. Cao, L., 1997. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D 110, 43–50. Gholipour, A., Abbaspour, A., Araabi, B.N., Lucas, C., 2003. Enhancements in the prediction of solar activity by locally linear model tree. Proceedings of MIC2003: 22nd International Conference on Modeling, Identification and Control, Innsbruck, Austria, pp. 158–161. Joselyn, J.A., 1997. Panel achieves consensus prediction of solar cycle. EOS Transactions American Geophysical Union 78, 211–212. Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining embedding dimension for phase space reconstruction
using a geometrical construction. Physical Review A 45 (6), 3403–3411. Laut, P., 2003. Solar activity and terrestrial climate: an analysis of some purported correlations. Journal of Atmospheric and Solar-Terrestrial Physics 65 (7), 801–812. Mees, A.I., Rapp, P.E., Jenning, L.S., 1987. Singular value decomposition and embedding dimension. Physical Review A 36 (1), 340–346. Nelles, O., 1999. Nonlinear system identification with local linear neuro-fuzzy models. Ph.D. Thesis, TU Darmstadt, Shaker Verlag, Aachen. Nelles, O., 2001. Nonlinear system identification. Springer, Berlin. Oh, H., Ammann, C.M., Naveau, P., Nychka, D., OttoBliesner, B.L., 2003. Multi-resolution time series analysis applied to solar irradiance and climate reconstructions. Journal of Atmospheric and Solar-Terrestrial Physics 65 (2), 191–201. Schatten, K., 2002. Solar activity prediction: timing predictors and cycle 24. Journal of Geophysical Research 107 (A11), 1377 doi:10.1029/2002JA009404. Schatten, K., 2003. Correction to Solar activity prediction: timing predictors and cycle 24 by Kenneth Schatten. Journal of Geophysical Research 108 (A3), 1100 doi:10.1029/ 2003JA009843. Schatten, K.H., Pesnell, W.D., 1993. An early solar dynamo prediction: cycle 23–cycle 22. Geophysical Research Letters 20, 2257–2278. Schatten, K.H., Sofia, S., 1987. Forecast of an exceptionally large even numbered solar cycle. Geophysical Research Letters 14, 632. Schatten, K.H., Scherrer, P.H., Svalgaard, L., Wilcox, J.M., 1978. Using Dynamo theory to predict the sunspot number during solar cycle 21. Geophysical Research letters 5, 411. Schatten, K.H., Myers, D.J., Sofia, S., 1996. Solar activity forecast for solar cycle 23. Geophysical Research Letters 23 (6), 605–608. Sofia, S., Fox, P., Schatten, K.H., 1998. Forecast update for activity cycle 23 from a Dynamo-based method. Geophysical Research Letters 25 (22), 4149–4152. Takens, F., 1981. Detecting strange attractors in turbulence. In: Rand, D.A., Young, L.S. (Eds.), Lecture Notes in Mathematics, Springer, Berlin, pp. 366–381.
ARTICLE IN PRESS A. Gholipour et al. / Journal of Atmospheric and Solar-Terrestrial Physics 67 (2005) 595–603 Thompson, R.J., 1993. A technique for predicting the amplitude of solar cycle. Solar Physics 148, 383. Tsiropoula, G., 2003. Signatures of solar activity variability in meteorological parameters. Journal of Atmospheric and Solar-Terrestrial Physics 65 (4), 469–482. Usokin, I.G., Solanki, S.K., Schussler, M., Mursula, K., Alanko, K., 2003. Millennium-scale sunspot number reconstruction: evidence for an unusually active Sun since
603
the 1940s. Physical Review Letters, doi: 10.1103/PhysRevLett.91.211101. Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynamics with applications to paleoclimatic time series. Physica D 35, 395–424. Vautard, R., Yiou, P., Ghil, M., 1992. Singular spectrum analysis: a toolkit for short noisy chaotic signals. Physica D 58, 95–126.