Modeling and Periodicity Analysis of Sunspot Time Series 1700-2015

This paper deals with the modeling of sunspot time series 1700-2015. Different models are fitted from four different classes of mathematical and stochastic models in order to describe this series. A special attention is made for the periodicity analysis of this series through these models as well as their main properties. Citation: Al-Khayat BYT, Hamdi MS (2018) Modeling and Periodicity Analysis of Sunspot Time Series 1700-2015. J Appl Computat Math 7: 385. doi: 10.4172/2168-9679.1000385


Introduction
Solar activity, as indexed by sunspots, is observations describing the average number of sunspots observed annually [1]. This series is highly asymmetric and has a cycle of length varying from 7 to 14 years [2]. Several models are fitted to sunspot series in mathematical, statistical and solar physical litrecheres in order to capture their main features as well as prediction [3]. This series consists of 317 observations describing the number of sunspots observed annually between 1700 and 2015 [4]. The data taken from the Royal Observatory of Belgium, SILSO, it is highly asymmetric as shown in Figure 1.
The main feature of this series is a cycle of activity varying in duration between 9 and 14 years, with an average period of approximately 11 years [5]. The series exhibits another feature, namely different gradients of "ascensions" and "dissensions", i.e., in each cycle the rise to the maximum has steeper gradient than the fall to the next minimum [6][7][8].
The era of linear modeling of sunspot time series began by Yules' autoregressive model in 1927. Box and Jenkins [1] fitted a second order linear autoregressive model, AR (2), to the period 1770-1869. Tong and Lim [8] fitted a non-linear model to the period 1700-1920. Subba and Gabr [3] have fitted a bilinear model to the period 1700-1920. Thanoon (1988) fitted a subset model of the model of Tong and Lim [8]. The eventual forecasting function of that model produces an asymmetric limit cycle of period 22 years. The limit cycle consists of two sub cycles whose "rise" and "fall" times are (5,6) and (4,7) years, respectively. Thanoon [5] who fitted a threshold model with piece wise non-linear dynamic to the same period 1700-1920 does another modification. The eventual forecasting function of that model has periodic nature but not systematic. It has a chaotic behavior, that behavior might be an indicator of chaotic behavior of the solar system. An adaptive spline threshold model fitted to the same period by Lewis and Stevens [2]. This model produces an asymmetric limit cycle of period 137 years. Li (2001) focused on predicting the magnitudes of the next cycle maximum as well as the next cycle minimum and the time from the initial minimum of a cycle to its maximum.
In this paper, we fit five different models to the complete sunspot series 1700-2015, one linear and four nonlinear models. These models are a delay model, a Henon map, a linear autoregressive model, and two non-linear autoregressive models. Orders of linear and non-linear autoregressive models are fixed at 11 and the first 20 observations are considered as initial values in all fitted models. The periodic property of sunspot series is studied through these models.

Model I
Delay difference equations are a type of difference equations in which the present of the unknown function at a certain time is given in terms of the values of the function at previous times. Many delay models are fitted to the data using the well-known least-squares method. Ignoring the noise term, the following delay model is then identified: x(n)=0.0126 × (n-1) [ 1-0.0037 × (n-8)] × (n-9). (1) A comparison between the real data and the obtained data from the fitted model is shown in the next ( Figure 2).
Using the forward Euler scheme; the previous model can be viewed as a discrete analogue of the following one-dimensional delay differential; where μ,β,α are some constants. The study of the solution of this delay equation is outside the scope of this paper.

Model II
The following dynamical system model, known as Henon map, is fitted to the series: where y(n)=0.8849x(n-1).
A comparison between the real data and the obtained data from the fitted model is shown in the next ( Figure 3).

Model III
An autoregressive (AR) model is a form of a stochastic difference equation; it is a representation of a time series such that the output variable depends linearly on its own previous values and on a stochastic term known as the "noise term".
Ignoring the noise term, this model can be written as:  (Table 1).
A comparison between the observed and the fitted series by the AR(11) model is shown in the following (Figure 4).
Threshold autoregressive models are non-linear time series models. These models are an extension of autoregressive models in order to allow for higher degree of flexibility in model parameters through a regime switching behavior. These models consist of more than one autoregressive parts (regimes), each for a different regime, Tong [6], Tong [7].
The graph of this power spectral density is shown in the following ( Figure 5).
It is clear that there is a dominant peak at the frequency f=0.59 radians, which corresponds to the frequency ω=0.59/2=0.0939 cycles per year. The corresponding period is P=1/0.0939=10.6496 years per cycle.
We may investigate the periodicity of sunspot series through the AR(11) model. It is known that an AR model could generate the pseudo-periodic behavior in the process when it has complex roots, Priestley (1981, P. 131). The characteristic equation of the AR(11) model has one real root, 0.9389, and the following ten complex roots: In order to investigate the periodicity of this series we consider the decomposition of the autoregressive operator in the following table. Roots of the characteristic equation are given in column (1). For real root r, the first order factor 1-r -1 is given in column 2, is the backward shift operator. For a pair of complex roots, the factor in column 2 is the associated second order factor, i.e. the factor associated with the roots a ± bi which is: In column (3) the frequency, in cycles per year, associated with each of these factors is given. For real roots this frequency is zero if Orders of difference equations are fixed at 11, two threshold autoregressive models to sunspot series 1700-2015 using least square method. The first model consist of two regimes and denoted by TAR(2; 11, 11), while the second consist of three regimes and denoted by TAR (3; 11, 11,11).

Model V
Ignoring the noise terms, the deterministic part of the fitted TAR(3; 11, 11, 11) model is written as:
The following Table 3 gives the residuals variances, the values of the Bayesian information criterion, BIC, as well as the long term behavior of each model.
According to the BIC criteria, the best model is the AR(11) which has the minimum value.

Investigation of Periodicity of Sunspot Series
We may find the spectral density function of the sunspot series based on the fitted AR(11) model, see; e.g.. Priestly (1981).    Table 3: Estimated autoregressive coefficients of the TAR(3; 11, 11, 11) model.
In column (4) the period, in years per cycles, p=ω -1 , is given. The last column of Table 4 indicates the location of each of the roots with respect to the non-stationary region, which are given by 2 2 a b + The harmonic nature of the fitted AR(11) model to the data is due to the fact that the characteristic equation of this model associated with the autoregressive operator has a pair of complex roots near the unit circle. From the last (Table 5), we note that the pair of complex roots associated with the frequency ω=09738 cycle per year is very close to the unit circle while no other roots are equally as close. Hence, the approximate 11 year cycle of sunspot series is detected by the AR(11) model.

Discussion
The cyclical nature as well as the eleven years period of the sunspot time series has been captured by the three time series models: AR(11), TAR(2; 11, 11) and TAR(3; 11, 11, 11).     We may also investigate the periodicity of sunspot series through the effectual forecasting function (eff). This function is obtained by generation the deterministic part of the model (by ignoring the noise terms) using last values of the series as initial conditions.
The effectual forecasting function of the TAR(2; 11, 11) model gives as limit cycle of period 11 years with ascend periods of 4 years and descend periods of 7 years, as shown in next ( Figure 6).
On the other hand, the effectual forecasting function of the TAR(3; 11, 11, 11) model has chaotic behavior, as shown in next (Figure 7).
The following Table 6 shows a comparison between the prediction performances for the next 11 years based on the fitted time series models.