Is there long-range memory in solar activity on time scales shorter than the sunspot period?

Is there long-range memor y in solar a ctivity on time scales shor ter than the sunspot period? M. Ryp dal 1 and K. Ryp dal 2 1 Departmen t of Mathematics and Statistics, Univ ersity of T romsø, Norwa y . 2 Departmen t of Physics and T ec hnology , Universit y of T romsø, Norwa y . No vem b er 7, 2018 Abstract The sunsp ot num b er (SSN), the total solar irradiance (TSI), a TSI reconstruction, and the solar flare index (SFI), are analyzed for long-range p ersistence (LRP). Standard Hurst analysis yields H ≈ 0 . 9, which suggests strong LRP . Ho wev er, solar activit y time series are non-stationary due to the almost p eriodic 11 year smooth com ponent, and the analysis do es not giv e the correct H for the sto c hastic component. Better estimates are obtained by detrended fluctuations analysis (DF A), but estimates are biased and errors are large due to the short time records. These time series can b e modeled as a sto chastic pro cess of the form x ( t ) = y ( t ) + σ p y ( t ) w H ( t ), where y ( t ) is the smooth comp onent, and w H ( t ) is a stationary fractional noise with Hurst exp onent H . F rom ensem bles of n umerical solutions to the stochastic model, and application of Bay es’ theorem, we can obtain bias and error bars on H and also a test of the hypothesis that a pro cess is uncorrelated ( H = 1 / 2). The conclusions from the presen t data sets are that SSN, TSI and TSI reconstruction almost certainly are long-range p ersisten t, but with most probable v alue H ≈ 0 . 7. The SFI pro cess, how ev er, is either very w eakly persistent ( H < 0 . 6) or completely uncorrelated. Some differences b etw een stochastic prop erties of the TSI and its reconstruction indicate some error in the reconstruction sc heme. 1 In tro duction As the Sun is the main energy source driving Earth pro cesses, v ariability of Solar output has b een at the center of scien tific interest for cen turies. Since Galileo’s inv en tion of the telescope, sunsp ots ha ve been systematically coun ted. It has b een noticed that the sunsp ot num b er (SSN) fluctuates in a more or less random fashion from mon th to mon th, but that their annual av erages v ary quasi-p erio dically with a p erio d of approximately 11 y ears. Other cycles (for instance on 22 and 80 y ears) are also discernible in the sunsp ot time series, and times of almost v anishing sunsp ot activit y hav e b een detected within the historical records. One example is the Maunder minim um 1630-1720 when sunspots almost disapp eared for nearly a century . Mo dern telescop es, and since the late 1970s also space observ atories, ha v e also facilitated record- ing the v ariations of solar-flare activit y , for instance through the solar flare index (SFI) and mea- sured v ariability of total solar irradiance (TSI) and more recen tly the solar spectral irradiance (SSI). V ariations of TSI obviously represen t a forcing of Earth climate, while the effect of v ariations of 1 SSI is more subtle, since differen t parts of the solar sp ectrum affect the Earth atmosphere in dif- feren t wa ys. Other climatic influences of solar v ariability are strongly debated, suc h as the effect of c hanges in the galactic cosmic ray flux due to v ariations in the Sun’s op en magnetic flux asso ciated with the solar wind. Muc h attention has b een dev oted to study correlation b etw een solar activity and climate on cen tennial time scales through v arious reconstructions of solar activit y based on theoretical mo dels with sunsp ot data as input and by means of paleo-data on millennial scales (see review by Usoskin (2008)). This in terest is obvious since a p ositive and quan titative identification of suc h a connection w ould allow a direct empirical estimate of climate sensitivity to solar activity as reviewed b y Haigh (2007); Gr ay et al. (2010). The idea of long-range persistence in the solar records, on short as w ell as long times scales, is not new. Mandelbr ot and Wal lis (1969) applied rescaled-range ( R/S ) analysis to monthly SSN and found the c haracteristic bulge on the log( R/S ) vs. log τ curve for time scales τ around the perio d of the sunsp ot cycle, but a common slop e corresp onding to a Hurst exp onent of H ≈ 0 . 9 in the ranges 3 < τ < 30 mon ths and 30 < τ < 100 years. R uzmaikin et al. (1994) obtained similar results for the SSN, and for analysis extended to the time range 100 < τ < 3000 y ears by using a 14 C proxy for cosmic ray flux. Since such b eha vior of the R /S -curv e is reproduced b y adding a sinusoidal oscillation to a fractional Gaussian noise with Hurst exponent H , b oth groups concluded that the sto c hastic comp onent of the SSN is a strongly p ersistent fractional noise in the time range from mon ths to millennia. Ogurtsov (2004) combined SSN and SSN reconstructions with 14 C pro xies to obtain H in the range 0.9-1.0 in the time range 25 < τ < 3000 years by means of R /S analysis and first and second order detrended fluctuation analysis (DF A(1) and DF A(2)) ( Huang et al. , 1998). The statistical significance of the long-memory p ersistence h yp othesis on time scales longer than the sunsp ot p erio d has b een questioned by Oliver (1998), using the so-called scale of fluctuation approac h. Another class of approaches to the sun-climate connection has fo cused on the detection of correlations and/or common statistical signatures b etw een proxi es of solar activit y and climate signals on time scales of the sunsp ot p erio d and shorter. Since the correlation b et ween the solar cycle v ariation of TSI and climate signals lik e global temp erature anomaly (GT A) app ears to b e quite weak, some works hav e drawn the attention to a p ossible “complexity-linking” b etw een solar activit y and Earth climate which is prop osed to b e discernible by identification of a common long- range memory process in the solar and climate s ignals, represen ted b y a common Hurst exp onent H ( Grigolini et al. , 2002; Sc afetta and West , 2003, 2005, 2008; Sc afetta et al. , 2004). This view w as criticized by us in a recent Letter ( Ryp dal and Ryp dal , 2010a), and triggered a comment by Sc afetta and West (2010) and a reply; Ryp dal and Ryp dal (2010b). Our view is that trends, lik e the sunsp ot cycle in solar signals and m ulti-decadal oscillations sup erposed on a rising trend in global temp erature signals, will create spuriously high Hurst exp onen ts. F or instance ( Sc afetta and West , 2005) (their Figure 3B) obtain H ≈ 0 . 95 for solar as well as global temp erature signals. After such trends are acc oun ted for, b oth solar signals and climate signals exhibit considerably lo wer Hurst exp onen ts for the sto c hastic comp onent of the signals. In the present pap er we will fo cus on three pro xies for solar activity , SSN, TSI, and SFI, and demonstrate that the SSN and TSI are the most p ersisten t (most probable H -v alues 0.73 and 0.70, resp ectively) on time scales shorter than the sunsp ot p erio d. F urther, for the SFI the most probable v alue is 0.54, and we cannot falsify the null h yp othesis that the sto chastic component in the solar flare signals exhibit no long-range memory ( H = 0 . 5) on time scales 10-1000 days. 2 2 Estimating Hurst exp onen ts from data 2.1 The sp ecial role of the v ariogram Let x (1) , x (2) , x (3) , . . . b e a stationary pro cess (a noise). W e supp ose that there exists a contin uous- time motion X ( t ), with stationary incremen ts, suc h that x ( n ) = X ( n + 1) − X ( n ). The q ’th order structure function S q ( t ) is the q ’th statistical moment S q ( t ) ≡ E [ X ( t ) q ] . If the motion X ( t ) is statistically self-similar the structure functions are p ow er la ws S q ( t ) ∝ t ζ ( q ) , and the scaling function ζ ( q ) has a linear dep endence on q . ζ ( q ) = H q , where H is the self-similarit y exponent. F or a self-similar motion the probability densit y function (PDF) of X ( t ) /σ ( t ), where σ 2 ( t ) ≡ S 2 ( t ) ∝ t 2 H , is the same for all t . If this PDF is Gaussian the process is called a fractional Bro wnian motion (fBm). F or an imp ortan t class of pro cesses, called m ultifractal, the PDFs change shap e with v arying time scale. They t ypically ha ve higher flatness (kurtosis) for small t , often conv erging tow ards Gaussian for large t . F or such pro cesses the structure functions are still p o wer laws, but the scaling function is not linear, but has a conv ex shap e. In this paper, we shall deal with an ev en wider class of pro cesses where only the second-order structure function, sometimes called the v ariogram, is assumed to b e a pow er-la w, i.e., S 2 ( t ) ≡ E [ X ( t ) 2 ] = E [ X (1) 2 ] t 2 H . (1) W e call H the Hurst exp onen t of the pro cess x ( t ). If X ( t ) is a self-similar pro cess, then H is its self-similarit y exponent. Thus, in general w e shall not require X ( t ) to be neither Gaussian or self- similar. In fact, w e do not ev en require that the process is m ultifractal since we only assume that the second structure function is a p ow er la w. The imp ortance of the Hurst exp onent is its relation to correlations in the noise x ( t ). W e include the deriv ation of this relationship for completeness, and start b y writing 2 X ( t ) X ( s ) = X ( t ) 2 + X ( s ) 2 − [ X ( t ) − X ( s )] 2 . (2) Since X ( t ) has stationary incremen ts, and using equation (1), w e ha v e that E [( X ( t ) − X ( s )) 2 ] = E [ X ( t − s ) 2 ] = E [ X (1) 2 ] | t − s | 2 H , and E [ X ( t ) 2 ] = E [ X (1) 2 ] t 2 H , E [ X ( s ) 2 ] = E [ X (1) 2 ] s 2 H , whic h, using equation (2), yields the co-v ariance structure E [ X ( t ) X ( s )] = 1 2 E [ X (1) 2 ]  t 2 H + s 2 H − | t − s | 2 H  . (3) Cho osing X (0) = 0 w e now hav e E [ x (1) x ( n )] = E h X (1) − X (0)  X ( n + 1) − X ( n ) i = E [ X (1) X ( n + 1)] − E [ X (1) X ( n )] , 3 and by using equation (3)w e get E [ x (1) x ( n )] ∝ ( n + 1) 2 H − 2 n 2 H + ( n − 1) 2 H ∼ d 2 dn 2 n 2 H = 2 H (2 H − 1) n 2 H − 2 . This means that the correlation function of x ( t ) has algebraic deca y for all H ∈ (0 , 1) except for H = 1 / 2, for whic h x ( t ) is an uncorrelated noise . F or H > 1 / 2 the in tegral ov er the correlation function is infinite, and this prop ert y is what is denoted as long-range p ersistence. Again, w e emphasize that this prop erty only requires a p ow er-law dep endence of the v ariogram, but not that the pro cess is self-similar or ev en m ultifractal. 2.2 Detrended fluctuation analysis As men tioned in the introduction solar signals contain distinct p erio dicities, the most prominen t b eing the 11-y ear solar cycle, and these trends will distort the result of the v ariogram analysis. In the next section we will demonstrate several w a ys of detrending and ho w these tend to reduce the estimated Hurst exp onents. Here we shall just briefly describ e one systematic metho d, the detrended fluctuation analysis (DF A), whic h p erforms quite well on these data. Again we consider the discrete and stationary sto chastic pro cess { x ( k ) } and we construct the cumulativ e sum pro cess X ( i ) = P i k =0 x ( k ). Let us divide the sequence X ( i ) in to N s segmen ts of length ∆ t , where each segmen t is indexed by l = 1 , 2 , . . . , N s . Then compute an n ’th order p olynomial fit p ( n ) ∆ t,l to X ( i ) in eac h segmen t and compute the v ariance of the detrended fluctuation in the segment, V ( n ) ∆ t,l ≡ 1 ∆ t ∆ t X i =1 [ X ∆ t,l ( i ) − p ( n ) ∆ t,l ] 2 . The square root of the a v erage of this v ariance ov er all the segments, F ( n ) (∆ t ) ≡ 1 N s N s X l =1 V ( n ) ∆ t,l ! 1 / 2 is the n ’th order DF A fluctuation function. If an fBm with Hurst exponent H ≈ 0 . 5 is sup erp osed on a sinusoidal signal of comparable amplitude, ordinary v ariogram analysis will give an estimated Hurst exp onent H ≈ 1, while a third order DF A analysis (DF A(3)) will give F ( n ) (∆ t ) ∝ ∆ t H with H close to the the Hurst exp onent for the fBm. Thus, there are goo d reasons to assume that giv en sufficiently long data series DF A(3) would giv e an accurate estimate of the Hurst exp onen t for the underlying detrended sto chastic pro cess. The main problem addressed in this pap er is how to reduce the uncertainties in the estimates that arise from the limited data sets. But first we demonstrate the effect of detrending on the solar data without dealing with uncertain ties. 2.3 The effect of trends In Figure 1a we hav e plotted the SFI, SSN and a TSI reconstruction ov er the last four sunsp ot cycles, and the instrumen tal TSI PMOD comp osite for the last tw o cycles. All data are given as 4 daily av erages. The thick smo oth curves are moving av erages weigh ted by a Gaussian windo w with 1 y ear standard deviation. In order to extract the stochastic comp onen t the smoothed signal should b e subtracted, but the signals are still strongly non-stationary b ecause the amplitudes of the rapid fluctuations are muc h higher at sunsp ot maxima than at minima. These amplitudes turn out to ha ve a distinct statistical dep endence on the lo cal moving av erage y ( t ). In fact, the v ariance of the fluctuations around y ( t ) is roughly prop ortional to y ( t ), i.e. V ar[ x ( t ) | y ( t ) = y ] ∝ y , as sho wn in Figure 2. This means that the mean- and amplitude-detrended signals z ( t ) ≡ ( x ( t ) − y ( t )) / p y ( t ) , (4) whic h is plotted in Figure 1b, are approximately stationary . The sunsp ot minima are still discernible in the detrended signals since flares and sunsp ots are absent for long time interv als around these minima. This implies that if statistics for the active sun is sough t, it ma y be reasonable to eliminate p erio ds of v ery low activity around sunsp ot minima from the analysis. In Figure 3a we hav e plotted v ariograms of the four raw signals x ( t ). This yields Hurst exponents in the range 0 . 88 < H < 0 . 97 and corresp onds to the results obtained in Figure 3B of Sc afetta and West (2005). In Figure 3b w e show results of DF A(3) analysis, which yields Hurst exp onents in the range 0 . 55 < H < 0 . 67. Results similar to the latter is obtained by computing v ariograms of the detrended signal in Figure 1b. 3 A sto c hastic mo del As mentioned in section 2.2 the limited length of the observed data records makes it difficult to compute error bars in the estimates of Hurst exp onents directly from the data. Essentially the en tire data set is used to estimate one v alue, and then w e do not hav e more data to assess the statistical spread in this estimate. What w e need to know to compute suc h error bars is the PDF of estimated v alues ˆ H in an imaginary ensemble of realizations of data sets of the same length as the observ ed record. Such an ensemble can b e generated syn thetically from a mo del that is assumed to hav e the the same statistical prop erties as the observ ed data, including an hypothesized v alue of H that can b e v aried. F rom the numerically generated ensembles one can construct a conditional probabilit y density p ( ˆ H | H ) for obtaining an estimated v alue ˆ H for the Hurst exp onent, given that the “true” exp onent is H . Then, by means of Ba y es’ theorem we can obtain the conditional PDF p ( H | ˆ H ), whic h gives us the probabilit y of ha ving a “true” Hurst exp onen t H provided w e hav e estimated a v alue ˆ H from the observ ational data. The width of p ( H | ˆ H ) giv es us the error bar of our estimate. The details of this pro cedure will b e giv en in the next section, but first we shall describ e the mo del used to generate synthetic data, which follo ws naturally from equation (4). W e ha v e already observed that z ( t ) defined from this expression is an approximately stationary sto c hastic pro cess and that v ariograms estimated from the data are appro ximate p ow er laws, from which a Hurst exp onent can b e estimated. PDFs of the signals on different time scales turn out to hav e appro ximately the same shape on the differen t time scales, but they are not alw a ys Gaussian. W e shall therefore model z ( t ) as a self-similar (in general non-Gaussian) pro cess with self-similarit y exp onen t H . By denoting this fractional noise pro cess b y w H ( t ), we can write x ( t ) = y ( t ) + σ p y ( t ) w H ( t ) , (5) 5 where y ( t ) is a 11-year-perio dic oscillation with amplitude A , and w H ( t ) is a fractional non-Gaussian noise with unit v ariance. The relation betw een parameters σ and A are easily estimated from the sunsp ot data by comparing the amplitude of the smo othed signal in Figure 1a with the amplitude of the detrended noise in Figure 1b. The fractional noise w H ( t ) is generated by taking w H ( t ) = Z H ( t + 1) − Z H ( t ) where Z H ( t ) is a pro cess defined b y the in tegral Z H ( t ) = c H Z K H ( t, t 0 ) d Z ( t 0 ) . (6) Here c H is a normalization constant and K H ( t, t 0 ) = ( t − t 0 ) H − 1 / 2 + − ( − t 0 ) H − 1 / 2 + . The pro cess Z ( t ) is chosen to b e a L ´ evy pro cess (a pro cess with stationary and indep endent in- cremen ts) with zero mean and unit v ariance. If Z ( t ) is gaussian, implying that it is a Bro wnian motion, then Z H ( t ) is a fractional Brownian motion with self-similarit y exp onen t H . In this case w H ( t ) is a fractional Gaussian noise with Hurst exp onent H . Using equation (6) together with the indep endence of incremen ts in Z ( t ) w e obtain the relation E [ Z H ( t ) 2 ] = c 2 H Z K H ( t, t 0 ) 2 dt 0 , and using the scaling relation K H ( at, at 0 ) = a H − 1 / 2 K H ( t, t 0 ) it follo ws that E [ Z H ( t ) 2 ] ∝ t 2 H . F or any given distribution 1 on Z ( t ) (and hence on w H ( t )) we can simulate the pro cess Z H ( t ) using the metho d of Sto ev and T aqqu (2004). F or the Sunsp ot Num b er the estimated PDF of the noise w H ( t ) is shown in figure 6(b). W e see that it is reasonable to choose a Gaussian distribution for w H ( t ). F or the other time series the Gaussian appro ximation is not v alid (See figures 4(b) and 10(b)), and here w e fit the data to stable distributions using a maximum likelihoo d estima- tor. These stable distributions are then truncated at v alues corresp onding roughly to the largest observ ed datap oint in the time series. The truncation ensures that the moments of w H ( t ) exist and preven ts the presence of unrealistically large even ts that can influence the analysis of esti- mated Hurst exp onen ts. T o demonstrate that our estimated distributions are reasonable we hav e compared the estimated PDFs for the fluctuations in the solar data with PDFs estimated from realizations of the sto chastic mo dels. These are presen ted in figures 4(b), 6b and 10(b). 4 A Ba y esian estimate 4.1 Uniform sample space W e assume that the observed signal is a realization of the sto chastic equation (5) where w H ( t ) is a fractional noise with the measured distribution and with Hurst exp onen t H ∈ [0 , 1]. W e construct a uniform (all v alues of H ∈ I ≡ [0 , 1] are equally probable) sample space S of realizations of these pro cesses. F or each realization with a giv en H we measure (estimate) a v alue ˆ H . In general ˆ H 6 = H . W e can then, for instance, compute the conditional probability densit y p ( ˆ H | H ) from an ensemble of n umerical solutions to the sto c hastic mo del where for eac h realization H is chosen randomly in 1 In order for Z ( t ) to b e a well-defined L ´ evy pro cess these distributions must b e infinitely divisible. How ev er, since w e do not require our mo del to hold on time scales shorter than a day , this has no practical imp ortance. 6 the in terv al I . The conditional PDF p ( H | ˆ H ) can also b e computed directly from the ensemble or from Bay es’ theorem: p ( H | ˆ H ) = p ( ˆ H | H ) p ( H ) p ( ˆ H ) . (7) Here p ( H ) is b y construction of S uniform on the unit interv al I and hence p ( H ) = 1. How ever, p ( ˆ H ) is not necessarily uniform and m ust b e computed from the ensem ble. 4.2 A non-uniform sample space The uniform sample space is based on a mo del of realit y prior to any observation of ˆ H that any v alue of H ∈ I is equally probable, i.e. w e hav e no prior preference to what H should b e. This ma y not b e the most plausible starting p oin t. Another is to ask the question; is the sto c hastic comp onen t in the signal long-range correlated or not? In other w ords, is H = 1 / 2 or is H 6 = 1 / 2? In the uniform sample space w e ha ve of course that p ( H = 1 / 2) = 0, so this sample space is not a goo d starting p oin t if w e wan t to assign a non-zero prior probabilit y to the h yp othesis that H = 1 / 2. What we can do is to construct a non-uniform sample space S ( p 0 ) consisting of an ensem ble of n umerical realizations of the fractional pro cess where a fraction p 0 of the realizations ha ve H = 1 / 2 and the remaining fraction 1 − p 0 is c hosen with H drawn randomly from the set I H 6 =1 / 2 = [0 , 1 / 2) S (1 / 2 , 1]. Let us formulate the following hypothesis: Hyp othesis H : The observ ed signal is describ ed by the sto c hastic equation with an uncorrelated ( H = 1 / 2) sto chastic noise term. W e then define an observ ation: Observ ation H o : DF A(3) applied to the observed signal yields an estimated Hurst exp onen t ˆ H = H o . F rom the constructed nonuniform ensemble we can compute the conditional probability that the h yp otesis H is true provided the prior probabilit y of this b eing true is p 0 ≡ p ( H ) and we ha v e made the observ ation by DF A of the physical time series that the estimated Hurst exponent is ˆ H = H o . First w e must create an ensem ble S with N elements. N p 0 elemen ts are realizations with H = 1 / 2. W e call this sub ensemble S H =1 / 2 . The remaining N (1 − p 0 ) elemen ts are realizations with H drawn randomly from I H 6 =1 / 2 . W e call this sub ensem ble S H 6 =1 / 2 . By doing this w e ha v e implicitly assumed zero probability to ha ve H outside the unit interv al, and hence that p ( H ) + p ( ¯ H ) = p 0 + (1 − p 0 ) = 1. Here, the sym b ol ¯ H denotes the null h yp othesis, H ∈ I H 6 =1 / 2 . The conditioned sub ensem ble S ( H o ) is the set of all elements for whic h the estimated ˆ H is in an  -neigh b orho o d of the observed v alue H o . Here  should b e a small n umber, but not smaller than the measurement error. Let us consider the fraction p ( H o |H ) of elements in S H =1 / 2 whic h belong to S ( H o ) and the fraction p ( H o | ¯ H ) of S H 6 =1 / 2 whic h belong to S ( H o ). These fractions are the conditional probabilities that ˆ H ∈ ( H o − , H o +  ) provided the hypotheses H or the h yp othesis ¯ H are true, resp ectively . It is now conv enient to in tro duce the Bay es factor: R ≡ p ( H o |H ) p ( H o | ¯ H ) . (8) 7 Note that if  is small, b oth numerator and denominator in this ratio are prop ortional to  , and hence R is indep endent of  . The Bay es factor is large ( R  1) if the probability of making the observ ation H = H o is m uc h larger if the hypothesis H is true than if ¯ H is true. Intuitiv ely then, making the observ ation H o should supp ort the hypothesis and make p ( H| H o ) close to unity . On the other hand, w e hav e that R  1 if the probability of H = H o is m uch larger if ¯ H is true than if H is true. In that case one should in tuitiv ely expect p ( H| H o ) to b e v anishingly small. This can b e quantified b y means of Ba yes’ theorem, whic h yields, p ( H| H o ) p ( ¯ H| H o ) = R p ( H ) p ( ¯ H ) = R p 0 1 − p 0 . Since p ( ¯ H| H o ) = 1 − p ( H | H o ) we can insert this in the equation abov e and solv e for p ( H| H o ): p ( H| H o ) =  1 + 1 R 1 − p 0 p 0  − 1 . (9) W e observ e that in addition to the Bay es factor also the prior probabilt y p 0 en ters the expression, but the result is not extremely sensitiv e to the c hoice of p 0 pro vided it is not very close to 0 or 1. The prior probability seeks to quantify our prior knowledge ab out the truth of the hypothesis H , and if that kno wledge is p o or, the most reasonable choice may b e p 0 = 1 / 2, which mak es (1 − p 0 ) /p 0 = 1 and p ( H| H o ) = R 1 + R . In general, the conditional probabilities entering the Bay es factor may b e computed from empirical data sets or from theoretical mo dels. In our case, we are facing v ery limited data sets, and w e therefore compute them from an ensemble of numerically generated realizations of the sto c hastic equation (5). The estimated Bay es factors for the v arious solar pro xies are presen ted in T able 1. 5 Results Since we are making Bay esian estimates, the results dep end on our prior probablities. Using the uniform ensem ble, this probability distribution is uniform on the unit interv al. Using the non- uniform sample space the prior probability parameter is p 0 . 5.1 The data The data records analyzed in the following are selected mainly for illustration of the p ow er and limitations of the analysis metho ds, although they also represen t interesting pro xies that reflect differen t asp ects of solar activit y . All records analyzed hav e differen t lengths, the longest is the sunsp ot num b er, and the shortest is the TSI comp osite. It will b ecome clear that the error bars in the estimates of the Hurst exp onents dep end on the lengths of the records, which is why the smallest error bars are found in the sunsp ot num b er and the largest in the TSI comp osite. It is conceiv able that other existing data sets could b e used to reduce these error bars, and certainly the near future will provide suc h data sets. The solar flare index (SFI) data ha ve b een do wnloaded from the NOAA National Geoph ysical Data Center (NNDGC) in Boulder, Colorado. The file contains daily data from the b eginning of 1966 to the end of 2008. There are no missing data points in this record. There is certainly a wealth 8 of other solar-flare data av ailable, but to our knowledge this is the longest unbrok en time-record of solar flare activit y with daily resolution. The daily sunspot n umber has been downloaded from the Solar Influences Data Analysis Cen ter (SIDC) in Brussels, Belgium. The file con tains data from 1818, but hav e missing data p oints up to 1848. After this year the data record is complete, so in our analysis we hav e used the unbrok en record from 1848 to 2011. The TSI reconstruction is based on Krivova et al. (2007) and the data file is retrieved from Max-Planc k Institute of Solar System Research (MPS) in Lindau am Hartz, Germany . The data file contains daily data from 1611, but with many gaps of missing data p oin ts b efore 1963. F or our analysis we therefore use the unbrok en record from 1963 to the end of 2004. The TSI comp osite data record is based on F r¨ olich et al. (2000) and is do wnloaded from Ph ysik alisches-Meteorologisc hes Observ atorium Dav os/W orld Radiation Center (PMOD/WR C), Switzerland, and is version d41 62 1007. The daily data set go es from Jan uary 1976 to July 2010, but con tains 7% missing data p oin ts. The missing data dominates in the early part of the record, and for this reason we hav e discarded all data b efore 1983. F or the data after 1983 there are only a small fraction of the data p oints missing. These missing data p oints hav e b een reconstructed from the sto chastic mo del equation (5), using a non-Gaussian white noise for w h ( t ). The probabilit y densit y function (PDF) for this noise term is estimated from the data, and the drift term y ( t ) is c hosen to b e the smo othed trend curve shown in figure 2. Since this smo othed trend curv e is based on a weigh ted mo ving av erage y ( t ) is prop erly defined only up to tw o years after the start and tw o y ears b efore the end the data record, so the time series analyzed in the follo wing go es from 1985 to 2008. 5.2 Uniform sample space Figure 4a sho ws the DF A(3) fluctuation function for the observ ed SFI time series (low er curv e) in a log-log plot. The estimated Hurst exp onen t is ˆ H ≈ 0 . 67. F or the SFI the DF A(3) estimate is particularly sensitiv e to the solar minim um p erio ds because in these in terv als there are long in terv als with no solar flares and zero SFI, which gives rise to spurious persistence of the signal. If interv als where the index is zero is excluded from the analysis on all time scales, w e obtain the upp er curv e in the figure, whic h yields ˆ H ≈ 0 . 55. W e b elive this figure is more representativ e for the solar-flare statistics in the active sun. Figure 4b shows the PDF of the detrended and amplitude-adjusted SFI (sho wn in Figure 1c) together with the PDF of the syn thetically generated SFI-signal. The dashed curv e is a fit of a stable distribution to the estimated PDF. The low er panel Figure 4c shows a realization of the syn thetically generated un-detrended signal. This signal should b e compared to the upp er trace in Figure 1a. In Figure 5a (crosses) w e show the mean estimated ˆ H computed from ensem bles of synthetic realizations of the SFI for 11 different v alues of H in the unit interv al. ˆ H is estimated from DF A(3) applied to the synthetic x ( t ) (i.e. the signal with the p erio dic trend) and the deviation from the straight dashed line indicates the systematic bias of the DF A(3) metho d applied to suc h a signal. When the pro cess is p ersisten t ( H > 0 . 5) the metho d p erforms quite w ell, with a slight o verestimation of H when the pro cess is strongly p ersistent ( H → 1). The circles (partly hidden b y the crosses) are computed from DF A(3) applied to the fractional noise directly . The fact that the tw o curv es fall on top of eac h other shows that the DF A(3) metho d remov es the influence of the p erio dic trend. There is also a statistical spread in the estimated ˆ H indicated by the error bars, whic h has a Gaussian conditional PDF p ( ˆ H | H ). This statistical spread is very imp ortant in 9 our analysis, since it giv es us the opportunity to use Bay es’ theorem to compute the error bars on the observed Hurst exponent. The spread in estimated ˆ H dep ends on the length of the synthetic data record, which we hav e c hosen equal to the observ ational record. As explained in the previous section the conditional PDF p ( H | ˆ H ) can be computed from the same ensem ble. In Figure 5b the conditional cumulativ e distribution P ( H | ˆ H ) ≡ R H −∞ p ( H 0 | ˆ H ) dH 0 is computed from the conditional PDF p ( ˆ H | H ) b y means of the uniform sample space and Bay es’ theorem for ˆ H = 0 . 55. F rom this distribution it is easy to compute the conditional mean E ( H | H o ) (which is the b est estimate for H given the observ ation H o ) and the 95% confidence interv al for this estimate. These are giv en in T able 6, and do es not rule out the p ossibility that the SFI sto chastic pro cess is uncorrelated ( H = 0 . 5). Note that E ( H | H o = 0 . 55) = 0 . 54, i.e. that the b est estimate for H and the observ ation H o are slightly different. Below, we shall see that this difference is more substantial in the sunsp ot n umber. Figures 6 and 7 sho w the same for the sunsp ot n umber. Here the PDFs are nearly Gaussian, and the 95% confidence in terv al is narro wer than for the SFI b ecause the observ ed sunspot record is m uch longer. The b est estimate for H is 0.73, while H o = 0 . 78, but the narrow confidence interv al excludes the possibility that H = 0 . 5. The results for the TSI are shown in Figures 8 and 9 and for the TSI reconstruction in Figures 10 and 11. The estimated Hurst exp onent is somewhat smaller ( H = 0 . 62) in the reconstruction compared to the observ ational data ( H = 0 . 70). The PDF also appears more skew in the former. These differences suggests some systematic errors in the reconstruction pro cedure that has effect on these timescales shorter than the sunsp ot p erio d. The reason for the larger confidence interv als for the TSI compared to the reconstruction, is that the latter time series is longer. The mean estimated Hurst exp onents and their 95% confidence in terv als are summarized in T able 1. Except for the SFI this analysis indicates that there is a significant p ersistence in all solar activit y pro xies analyzed, but their Hurst exp onents are considerably low er than obtained by most previous authors. The largest H is found in the sunsp ot num b er, with the most probable v alue H = 0 . 73, and with 95% confidence H < 0 . 76. 5.3 Non-uniform sample space F rom T able 1 it is tempting to conclude that H > 1 / 2 also for the SFI, since the H = 1 / 2 v alue is barely within the 95% confidence interv al. Ho w ever, this result is based on an assumption on a uniform prior probabilit y , i.e. we as sume all v alues of H ∈ I is equally probable prior to the observ ation. If, how ever, prior to observ ation we assume that there is a finite probability that the SFI-pro cess is uncorrelated ( p 0 6 = 0) on the relev ant time scales, section 4.2 shows that the computation is differen t and the result will dep end on p 0 . The conditional probabilit y p ( H| H o ) is giv en by equation(9), where the Bay es factor R must b e computed from the numerically generated ensem ble of realizations generated for given v alues of p 0 . The curves of p ( H| H o ) as a function of p 0 are shown for TSI and SFI in Figure 12. It appears that the conditional probabilit y that H = 1 / 2 is negligible for the TSI , unless we assume a prior probability for this very close to unit y (in that case our prior prejudice ov ershadows an y computation). How ever, for the SFI we observ e that even an agnostic prior view ( p 0 = 0 . 5) will result in a higher probability after observ ation p ( H | H o ) ≈ 0 . 8 and that p ( H| H o ) > 0 . 8 if p 0 > 0 . 5. 10 6 Discussion and conclusions The emphasis in this paper is on analysis metho d for discerning long-range memory in short data records with near-p erio dic trends in lo cal mean and v ariance. The data records analyzed hav e b een c hosen mainly for illustration of these metho ds. Still, from these data w e can conclude that Hurst exp onen ts for all pro xies analyzed are significan tly smaller than obtained recently by Sc afetta and West (2005) and in the classical papers b y Ruzmaikin et al. (1994) and Mandelbr ot and Wal lis (1969). W e can also conclude that TSI and SSN time series on the one hand are considerably more p ersisten t than the SFI time series, and at present it cannot be excluded that the solar flare activit y is an uncorrelated pro cess on time scales longer than a few mon ths. W e also hav e indications that the particular TSI reconstruction analyzed here (see T able 1) exhibits lo wer p ersistence and a sk ewer PDF of the fluctuations than the observ ed TSI. More accurate and longer reconstruction time series are emerging, and as the TSI time series are also growing longer, the metho ds devised here should b e able to pinp oin t differences in the stochastic properties of the observed and reconstructed signals that might hav e to b e addressed b y the groups working with reconstructions. As measurements of solar sp ectral irradiance (SSI) o ver sev eral solar cycles b ecome av ailable, time series of irradiance in the UV-range will b e interesting to analyze. P ersistence in TSI and SSI will b e in teresting to compare to corresp onding prop erties in solar wind parameters, geomagnetic indices and terrestrial climate parameters in order to assess the existence of p ossible sun-climate coupling ( Sc afetta and West , 2003; Ryp dal and Ryp dal , 2010a), and to identify the coupling mech- anisms. Sto c hastic prop erties of the proxy data at solar maxima could b e compared with those at solar minima, and if a systematic difference exists, suc h measures of reconstructed data near the Maunder and Dalton minima should b e compared to those aw a y from these minima. The ob jective of this exercise would b e to establish whether sto chastic measures of pro xy data could b e used as reliable pro xies for solar activity . References F r¨ olic h, (2010), Observ ations of Irradiance V ariations, Sp ac e Scienc e R ev. , 94 , 15. Gra y , L. J., et al. (2010), Solar influences on climate, R ev. Ge ophys , 48 , R G4001. Grigolini, P ., et al. (2002), Diffusion en tropy and w aiting time statistics of hard-x-ra y solar flares, Phys. R ev. E , 65 , 046203. Haigh J. D., (2007),The Sun and the Earths Climate, Living R ev. Solar Phys. , 4 , 2. [Online Article]: h ttp://www.livingreviews.org/lrsp-2007-2 Huang, N.E, Z. Seng, S. R. Long, M. C. W u, H. H. Shih, Q. Zheng, N. Y en, C. C. Th ung, H. H. Liu (1998), The empirical mo de decomp osition and the Hilb ert sp ectrum for nonlinear and non-stationary time series analysis, Proc. R. So c. Lond. A 454 , 903. Kriv ov a, N. A., L. Balmaceda, and S. K. Solanki (2007), Reconstruction of solar total ir- radiance since 1700 from the surface magnetic flux, Astr onomy & Astr ophysics , 467 , 335, h ttp://www.mps.mpg.de/pro jects/sun-climate/data.html 11 Mandelbrot, B. B., and J. R. W allis (1969), Global dependence in geophysical records, Water R esor. R es. , 5 , 321. Ogurtso v, M. G. (2004), New evidence for long-term persistence in the Sun’s activit y , Solar Physics , 220 , 93. Oliv er, R., and J. L. Ballester (1998), Is there memory in solar activity?, Physic al R eview , 58 , 5650. Ruzmaikin, A., J. F eynman, and P . Robinson (1994), Long-term p ersistence of solar activity , Solar Phys , 148 , 395. Ryp dal, M., and K. Ryp dal, (2010), T esting Hypotheses ab out Sun-Climate Complexit y Linking, Phys. R ev. L ett. , 104 , 128501. Ryp dal, M., and K. Rypdal, (2010), Ryp dal and Ryp dal reply , Phys. R ev. L ett. , 105 , 219802. Scafetta, N., and B. J. W est (2003), Solar Flare In termittency and the Earths T emperature Anoma- lies, Phys. R ev. L ett. , 90 , 248701. Scafetta, N., et al. (2004), Solar turbulence in earths global and regional temp erature anomalies, Phys. R ev. E , 69 , 026303. Scafetta, N. and B. J. W est (2005), Multiscaling Comparativ e Analysis of Time Series and Geo- ph ysical Phenomena, Complexity , 10 , 51. Scafetta, N. and B. J. W est (2008), Is climate sensitiv e to solar v ariability , Physics T oday , p. 50, Marc h 2008. Scafetta, N., and B. J. W est (2010), Comment on “T esting Hyp otheses ab out Sun-Climate Com- plexit y Linking, Phys. R ev. L ett. , 105 , 218801. author = Stilian Stoev and Murad S. T aqqu, title = Simulation metho ds for linear fractional stable motion and F ARIMA using the fast F ourier transform, journal = F ractals, y ear = 2004, v olume = 12, pages = 95–121 Sto ev S., and M. S.T aqqu (2004), Simulation metho ds for linear fractional stable motion and F ARIMA using the fast F ourier transform, F r actals , 12 , 95. Usoskin, I. G., S. K. Solanki, and G. A. Ko v altsov (2007), Grand minima and maxima of solar activit y: new observ ational constraints, Astr onomy & Astr ophysics , 471 , 301. Usoskin, I. G., (200), A History of Solar Activit y ov er Millennia, Living R ev. Solar Phys. , 5 , 3. [Online Article]: http://www.livingreviews.org/lrsp-2008-3 12 Figure 1: (a): Time series for Sunsp ot Number (SN), Solar Flare Index (SFI), T otal Solar Irradiance (TSI) reconstruction, and the TSI PMOD comp osite. The smo oth curves are obtained b y conv olving the time series with a gaussian kernel with a one-year standard deviation. The data are normalized suc h that the smo oth comp onen ts hav e equal maximal amplitude, i.e. the difference b etw een the maximal and minimal v alue is the same for the four smoothed signals. (b): In this figure we hav e subtracted the smo othed signals from the original times series in order to de-trend the time series and then divided the detrended signals by the square ro ot of the smoothed signals in (a). Prior to this amplitude adjustmen t, the origins for eac h of the signals are shifted to the minimal v alues of the smo oth signals in (a). 13 Figure 2: F or eac h of the signals x ( t ) in Figure 1a the plots sho w the v ariance V ar[ x ( t ) | y ( t ) = y ] conditioned up on the v alue of the corresp onding smo othed signal y ( t ). When the v ariance’s dep endence on the smo oth signal of the fluctuations are well appro ximated by a linear function, we can mo del the standard deviation’s dep endence on the smo oth signal b y the square-ro ot function. TSI reconstruction TSI (PMOD) Sunsp ot Number Solar Flare Index H o 0 . 61 0 . 71 0 . 78 0 . 55 R H p ( H | H o ) dH 0 . 62 0 . 70 0 . 73 0 . 54 95% confidence 0 . 57 < H < 0 . 69 0 . 62 < H < 0 . 81 0 . 68 < H < 0 . 76 0 . 49 < H < 0 . 61 R 7 . 0 × 10 − 3 1 . 6 × 10 − 2 3 . 8 T able 1: 1. row: Hurst exp onen t H o as estimated b y DF A(3) from the observ ed time series. 2. row: The most probable Hurst exp onents R H p ( H | H o ) dH computed from the uniform sample space. 3. ro w: The 95% confidence interv als computed from the uniform sample space. 4. row: The Bay es factor computed from the non-uniform sample space. The computaion for the SSN is not shown since Bay es factor in this case is obviously extremely small. 14 Figure 3: V ariograms computed from the raw (undetrended) data. (b): DF A(3) v ariogram (square of fluctuation function) computed from raw data. The v alues for H giv en in the figure is half the slop e of the dotted lines. These lines are dra wn to indicate the typical slop e of the corresponding v ariograms. 15 Figure 4: (a): DF A(3) fluctuation function squared computed from the SFI ra w data (circles), and DF A(3) fluctuation function computed from the SFI ra w data, but b y omitting all in terv als on all time scales that contain data p oints where the SFI is zero. (b): PDF of fluctuations the detrended SFI data series (circles), and PDF of the corresp onding data generated by the sto chastic mo del. (c): A realization of a model-generated SFI time series. 16 Figure 5: (a): Ensem ble mean E [ ˆ H ] of estimated Hurst exponent for SFI computed b y the DF A(3) metho d. The crosses are computed from an ensem ble of realizations generated from the stochastic mo del where the fractional noise w H ( t ) used in the mo del has Hurst exp onent H and the PDF giv en in Figure 4b. The circles (partly hidden b y the crosses) are computed from DF A(3) applied to the fractional noise directly . The Gaussian statistical spread in the estimated ˆ H is given by the error bars. (b): The conditional cum ulativ e distribution P ( H | ˆ H ) ≡ R H −∞ p ( H 0 | ˆ H ) dH 0 computed from the conditional PDF p ( ˆ H | H ) b y means of the uniform sample space and Ba yes’ theorem for ˆ H = 0 . 55. 17 Figure 6: (a): DF A(3) fluctuation function squared computed from the SNN ra w data (circles). (b): PDF of fluctuations the detrended SSN data series (circles), and PDF of the corresp onding data generated b y the sto c hastic mo del. (c): A realization of a model-generated SSN time series. Figure 7: Same as Figure 5, but for SSN time series. The reason for the narro wer confidence interv al in panel (b) compared to those in Figure 5 is that the SSN time series analyzed is considerably longer than the SFI time series. 18 Figure 8: DF A(3) fluctuation function squared computed from theTSI raw data (circles). (b): PDF of fluctuations the detrended TSI data series (circles), and PDF of the corresp onding data generated by the stochastic mo del. (c): A realization of a model-generated TSI time series. 19 Figure 9: Same as Figure 5 and Figure 7, but for the TSI time series. The reason for the wider confidence interv al in panel (b) compared to those in Figure 5 and 7 is that the TSI time series analyzed is considerably shorter than the SFI and SSN time series. 20 Figure 10: DF A(3) fluctuation function squared computed from the TSI reconstruction ra w data (circles). (b): PDF of fluctuations the detrended TSI reconstruction data series (circles), and PDF of the corresponding data generated b y the sto c hastic mo del. (c): A realization of a mo del-generated TSI reconstruction time series. 21 Figure 11: Same as Figure 9, but for the TSI reconstruction time series. The reason for the slightly narro wer confidence in terv al in panel (b) compared that in Figure 9 is that the TSI reconstruction time series analyzed is somewhat longer than the TSI time series. 22 Figure 12: The conditional probability p ( H| H o ) versus prior probability p 0 for (a):TSI, and (b): SFI, computed from the non-uniform sample space. 23