Extreme Value Theory Applied to the Millennial Sunspot Number Series

The Sun is a variable star, and studies of variations of solar activity may shed light on the magnetic activity of cool stars. Solar magnetic activity is observed for the last century as synoptic images of the Sun, and for the last 400 years in the form of sunspot counts and drawings (see, e.g., review papers by Vaquero 2007; Clette et al. 2014; Hathaway 2015; Usoskin 2017). Solar activity depicts a great deal of variability during the last centuries, from the very quiet Maunder minimum in the second half of seventeenth century (Eddy 1976; Usoskin et al. 2015) to the modern high-activity episode in the second half of the twentieth century (Solanki et al. 2004). Indirect proxies, such as cosmogenic radionuclides in terrestrial archives, can help in reconstructions of solar activity for several millennia in the past (Beer et al. 2012; Usoskin 2017). However, there is still an open question as to whether the period of the last 400 years (since 1610 AD) covered by direct solar observations is representative for the entire range of solar variability. In other words, are the observed Maunder minimum and the Modern grand maximum "typical" for solar activity variations?

In some branches of science, extreme values of significant variables have a special value and meaning. In those cases, one needs to use a statistical theory devoted specifically to an analysis of rare values corresponding to extreme situations. One of the options is the use of the extreme value theory (EVT). In recent years, the EVT has been applied in terrestrial (Acero et al. 2014; Wi et al. 2016) and solar climatology (Asensio Ramos 2007; Acero et al. 2017). Extreme values of solar activity have been studied directly using the sunspot number series recoded during the last four centuries (Asensio Ramos 2007; Acero et al. 2017) and indirectly using geomagnetic indices compiled during the last recent decades (Siscoe 1976). However, studies using longer timescales have not been performed until now.

Here, we analyze reconstructed multi-millennial records of solar activity, based on cosmogenic isotopes, by modern statistical methods, viz. EVT, to evaluate the properties of the extremes of solar activity on the long timescale.

Two sets of cosmogenic radionuclide data were used as tracers of solar activity (Beer et al. 2012; Usoskin 2017): ¹⁴C in tree trunks and ¹⁰Be in polar ice. The decadal sunspot numbers reconstructed from both radionuclides were considered in this study, as published by Usoskin et al. (2016), and denoted as SN 14C and SN 10Be, respectively. The study period is considered as from 6755 BC through 1645 AD for SN 10Be and from 6755 BC to 1895 AD for SN 14C (Figure 1).

Zoom In Zoom Out Reset image size

EVT is a branch of statistics dealing with the distribution of excesses, trying to study and quantify the behavior of a process at unusually large or small levels. It seeks to assess, from a given ordered sample of a given variable, the probability of events that are more extreme than any previously observed (Coles 2001), being one of the statistical disciplines most commonly used in the last few decades in many fields such as finance, hydrology, and life and earth sciences. EVT is used to study extreme values in a timeseries and aims to predict the occurrence probability of rare events. Among different approaches used in the EVT, we consider the peaks-over-threshold (POT) one. This technique is based on the definition of a high enough threshold and a fit of the exceedances over the threshold to the prescribed statistic.

The POT technique considers all sample values that exceed a predefined upper threshold u. The probability distribution of the exceedances over the threshold can be modeled using the generalized Pareto distribution (GPD). The POT approach needs independent observation to avoid too short-range dependences in the timeseries, when one data point is linked to neighboring ones, and extreme values may cluster together. The cosmogenic radionuclide timeseries studied in this work show independent observations separated by 10 years. Therefore, it is not necessary to apply a declustering procedure usual in this technique. Both timeseries were subjected to a GPD analysis. In the asymptotic limit for sufficiently large thresholds, with the observed radionuclide timeseries SN(t), the distribution of independent exceedances X_u(t) = SN(t)−u with SN(t) > u follows a GPD given by

$$\begin{eqnarray}&&P(X\lt x)=1-{\left(1+\displaystyle \frac{\xi \cdot x}{\sigma }\right)}^{-1/\xi }.\end{eqnarray} \tag{ 1 }$$

with x > 0 and 1+ξx/σ > 0, where σ is the scale parameter, and ξ is the shape parameter ($\xi \ne 0$). Negative shape parameter values indicate that the distribution has an upper bound, while values positive or zero indicate that the distribution has no upper limit (Coles 2001). The scale parameter σ gives information about the variability of the distribution.

To apply the above approach to the EVT, it is necessary to choose a correct threshold of extremes from which the exceedances must be evaluated. This threshold u will be not too high in order to the number of exceedances to be large enough to minimize the uncertainty of the GPD parameters, but not too low because it would violate the asymptotic basis of the model (Coles 2001). Two methods are available for the threshold selection: the mean residual life plot and the assessment of the stability of the parameter estimates. First, the mean residual life plot was considered, which involves plotting the "mean excess" (the mean value of X_u(t)) against u, for a range of values of u. Such a plot is expected to be linear above the threshold at which the GPD model becomes valid. But as mentioned by Coles (2001), the interpretation of the mean residual life plot is not always simple in practice. Figure 2 shows the mean residual life plot for both radionuclides (black line) and its confidence interval (CI) (dashed lines) based on the approximate normality of sample means (Coles 2001). Note that there is an approximately linear relationship for the intervals [54, 59] and [55, 60] of the threshold u for SN 14C and SN 10Be, respectively, just when the CIs became larger. For lower values of the threshold u, there is another linear zone. However, these lower values of u provide an excessive number of exceedances that do not allow the use of the EVT, as aforementioned. Second, the parameter stability plot involves plotting the parameter estimates from the GPD model against u, for a range of values of u. These parameters are the shape parameter ξ and the reparametrized scale parameter σ* = σ_u − ξu (see Section 4.3.4 of Coles 2001, for details). The parameter estimates should be stable (i.e., near constant) above the threshold at which the GPD model becomes valid. The two bottom panels in Figure 2 show the parameter plots and their CIs for SN 14C (left) and SN 10Be (right). This method also confirmed that u in the interval [54, 59] for SN 14C and [55, 60] for SN 10Be timeseries were optimal thresholds, keeping the shape and scale parameters near constant. Then, from both methods, the thresholds selected were [54, 59] for SN 14C and [55, 60] for SN 10Be which correspond approximately with the values in the interval ranged from the ninetieth to the ninety-fifth percentile for each timeseries.

Zoom In Zoom Out Reset image size

Therefore, the independent exceedances for both radionuclides timeseries were subjected to a GPD analysis. The parameters of the GP distribution in Equation (1) were estimated by the maximum likelihood method using the in2extRemes statistical R software package for extreme values (Gilleland & Katz 2013). Once the parameters were estimated, the CI for each parameter was evaluated by a bootstrap procedure using 500 replicates (Gilleland & Katz 2013).

Besides, an important property to predict the probability of the occurrence of future extreme events is the return level (RL), which is commonly used to convey information about the likelihood of rare events. The N-year RL is the level expected to be exceeded once every N years on average, and it was estimated for different values of N using both approaches. More details about RL estimations using POT can be found in Coles (2001). The aforementioned procedure was also used to estimate the RLs and their CIs with the bootstrap procedure.

Table 1 shows the results for different thresholds analyzed and the corresponding number of exceedances for each threshold, ranging u from 54.0 to 58.6 for SN 14C and from 56.0 to 60.4 for SN 10Be. Tables 2 and 3 show the results for the GPD parameters when fitting the data to that distribution for SN 14C and SN 10Be, respectively, and for the different thresholds chosen. These tables list the estimates for the scale and shape parameters and the 95% CIs obtained by bootstrapping. The shape parameter for the different thresholds appears mostly negative, even within the 95% CI, implying an upper bound of the extreme values for the SN 10Be and SN 14C radionuclide timeseries. This result is more relevant for SN 10Be because for all the thresholds considered, except for u equal to the ninety-fourth percentile, the shape parameter is negative within the 95% CI. For SN 14C, the shape parameter is almost always negative, but the 95% CI sometimes slightly extends to positive values. These results imply, for both radionuclides, that the highest extreme values in relation to solar activity have been reached in the past and are not expected to be exceeded in the future. In order to assess the accuracy of the threshold excess model fitted to each radionuclide timeseries, different diagnostic plots were used. Figure 3 shows diagnostic plots for the GP fit to the threshold exceedances from the SN 10Be radionuclide timeseries. The plots confirm the validity of the fitted model: top panel in Figure 3 shows a quantile–quantile (QQ) plot of empirical data quantiles against GP fit quantiles leading to similar distributions with the points lying on the line y = x (solid line). Middle panel in Figure 3 shows a QQ-plot of randomly generated data from the fitted GP against the empirical data quantiles with 95% confidence bands with the points being nearly linear too. Besides, the corresponding density estimate seems consistent between the empirical density of the observed maxima SN with the modelled GP fit density (bottom panel in Figure 3). Figure 4 also shows the diagnostic plot for the GP fit but for the SN 14C timeseries also validating the fitted model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned above, a usual procedure for interpreting extreme values uses the RL. As the original timeseries recorded decadal observations, then the N -decadal RL was estimated for N  = 100, and 1000 decades, corresponding to RLs of 1000 and 10,000 years, respectively, and the 95% CI was estimated using the bootstrap technique. Tables 4 and 5 list estimates of the RLs for the two radionuclides timeseries and for the two values of N chosen. The result is not trivial. One would expect that the 1000 year RL lies inside the exceedances range considered for both radionuclides: [54.02, 78.84] for SN 14C, and [55.97, 76.77] for SN 10Be for the lowest value of the threshold considered, because the observed timeseries is longer than 1000 years. However, the results show that even the 10,000 year RL corresponding to values of N greater than the observed period lie actually also inside the mentioned intervals. For SN 14C and considering the different thresholds, the 10,000 year RL varies from 73.9 and 74.5, values lower than the highest-observed value 78.8. Similarly for SN 10Be, the 10,000 year RL ranges between 74.2 and 76.4, which is also lower than the highest-observed value 76.8. The RL plots are shown in Figure 5 for both radionuclides. An extrapolation for the 1000 decade RLs given by the solid line indicates a value lower than the highest-observed one, although the maximum for both radionuclides are inside the 95% CI uncertainty of the RL. Therefore, an interesting result is that the RL is reaching a plateau, and it will be unlikely that values greater than the already observed ones will be reached in the future. These results must be considered with some caution, however, because RLs were estimated for a longer period than the observed one.

Zoom In Zoom Out Reset image size

In this work, we applied the POT technique, in the context of the Extreme Value Theory, to SN 14C and SN 10Be timeseries for the last nine millennia. The shape and the scale parameter for the generalized Pareto distribution, which models the probability distribution of the exceedances, were estimated. Different thresholds were considered and validated to be accepted as good thresholds in the sense of the EVT. As aforementioned, the information provided by the shape parameter is relevant for an understanding of the behavior of the timeseries in the sense of extreme values. In this study, the shape parameter was found mostly negative for all the thresholds considered and for both radionuclides, thus revealing the existence of an upper bound for the extremes of the timeseries. Besides, to interpret the extreme values, the N-year RLs were estimated for two values: N = 1000 and 10,000 years. As expected, the 1000 year RL lies inside the exceedances range considered for both radionuclides, but surprisingly the 10,000 year RL also lies in the exceedances range despite of the fact that the observed period is shorter than 10,000 years, being probable a higher value. These results suggest that the highest extreme values of the timeseries for both radionuclides in relation to solar activity have been reached in the past and are not expected to be exceeded in the nearest future. Moreover, the RLs are reaching a plateau, and it will be unlikely that the sunspot numbers will attain values greater than the already observed ones in the future. This generally agrees with a similar analysis carried out using the sunspot number based in telescopic data for the last four centuries (Acero et al. 2017).

Thus, the main characteristics of the solar activity have already been observed in the telescopic period (the last four centuries), which covers the full range of solar variability from a grand minimum to a grand maximum. In any case, this approach and the results are limited by the decadal resolution of the isotope data used in this work.

Support from the Junta de Extremadura (Research Group Grant GR15137) and from the Ministerio de Economía y Competitividad of the Spanish Government (AYA2014-57556-P) is gratefully acknowledged. I.U. acknowledges support by the Center of Excellence ReSoLVE (Academy of Finland Project No. 272157).