Automatic algorithm for monotone trend removal Călin Vamoş “T. Popoviciu” Institute of Numerical Analysis, Romanian Academy, P.O. Box 68, 400110 Cluj-Napoca, Romania Received 6 October 2005; revised manuscript received 29 November 2006; published 16 March 2007 The available numerical algorithms for trend removal require a direct subjective intervention in choosing critical parameters. In this paper an algorithm is presented, which needs no initial subjective assumptions. Monotone trends are approximated by piecewise linear curves obtained by dividing into subintervals the signal values interval, not the time interval. The slope of each linear segment of the estimated trend is proportional to the average one-step displacement of the signal values included into the corresponding subinterval. The evalu- ation of the trend removal is performed on statistical ensembles of artificial time series with the random component given by realizations of autoregressive of order one stochastic processes or by fractional Brownian motions. The accuracy of the algorithm is compared with that of two well-tested methods: polynomial fitting and a nonparametric method based on moving average. For stationary noise the results of the algorithm are slightly better, but for nonstationary noise the preliminary results indicate that the polynomial fitting has the best accuracy. As a verification on a real time series, the time periods with monotone variation of global average temperature over the last 1800 years are established. The removal of a nonmonotone trend is also briefly discussed. DOI: 10.1103/PhysRevE.75.036705 PACS numbers: 05.40.Ca, 02.60.Ed I. INTRODUCTION In many practical situations processing experimental ob- servations of superposed phenomena having different time scales is needed. Statistical methods can be used to separate these components from the global signal. For example, trend removal is performed by different methods from the simpler ones least-squares fit of a parametric family of functions or smoothing by moving averageto more complex as discrete differencing 1or wavelet analysis 2. More sophisticated methods should be used only if their results are significantly better than the results of the simpler ones. A common draw- back of the available methods is a direct subjective interven- tion in choosing critical characteristics such as functional form of the trend, number of averagings or differentiations of the signal, length of the averaging interval, type of the basis functions, etc. In this paper a trend removal algorithm, which needs no initial assumptions, is proposed. It is based on a statistical method, which approximates the trend by a piecewise linear curve obtained by dividing into subintervals the signal values interval, not the time interval. The slope of each linear seg- ment of the estimated trend is proportional to the average one-step displacement of the signal values included into the corresponding subinterval, therefore the method is referred to as average conditional displacement ACD. Here only the case of the monotone trend is considered. The extension of the ACD method to a nonmonotone trend is discussed in the last section and will be presented in a separate future paper. In practice there are situations when the trend monotony is specially required, for example, ascertainment of the time periods with monotone variation of the global atmospheric temperature, problem discussed in Sec. V. But even when the trend monotony is not significant, the initial separation of a monotone component from the global nonmonotone trend may be useful. Even if initially there is no indication that a persistent phenomenon has a contribution to the generation of the time series, however, rendering evident a significant monotone component is a reason to test the existence of such a phenomenon. Obviously, the main interest in real applica- tions is to estimate the whole trend, including the turning points. As a rule, the trend removal methods do not take into account whether the estimated trend is or is not monotone. The exceptions are the methods in which the trend is explic- itly looked for as a monotone function, usually a linear, ex- ponential, or logarithmic function. The weakness of this ap- proach consists in the limited number of available monotone functional forms. If an enrichment of the functional forms is attempted, for example, using polynomials of order greater than one, then the monotony property is lost. From this point of view, the advantage of the ACD method is that it can describe a much richer set of monotone trends as piecewise linear functions. Furthermore, if the ACD method does not succeed to remove a monotone trend, this is useful informa- tion to characterize the time series properties because it shows that no monotone component could be associated to that series. Besides the request that the estimated trend should be monotone, another essential feature of the ACD method is its capacity to support the development of an automatic algo- rithm. The automatic choice of the values of the ACD pa- rameters is performed by means of a numerical estimation of the standard deviation of the signal random component. The quality of the estimated standard deviation and the accuracy of the automatic ACD algorithm are tested on artificial sig- nals with stationary noises of type AR1autoregressive of order one. Even without a thorough analysis, the ACD method is also applied on signals with nonstationary noise. Such sig- nals occur in detrended fluctuation analysis DFA3, ex- tensively used for the detection of long-range correlations in time series with deterministic trend and stationary 1 / f noise. A detailed analysis of DFA and a review of its applications are presented in 4. The initial signal is summed up, result- PHYSICAL REVIEW E 75, 036705 2007 1539-3755/2007/753/03670516©2007 The American Physical Society 036705-1
ing in a secondary one with nonstationary noise. An essential step in DFA is the repeated removal of polynomial trends by means of the least-squares method from the secondary signal and DFA could be significantly improved if the trends were removed by means of an automatic algorithm like the exten- sion of the ACD method for the nonmonotone trend. Since ACD is a method which does not develop any of the trend removal methods used at present, there is no obvious term of comparison for its performance evaluation. Therefore I have chosen two reference methods based on secondary criteria. The first one is the polynomial fitting with different degrees used in DFA. The second one is “the jump process” presented in 5which has some similarities with the method used in the ACD algorithm to dampen the fluctuations. Thus I have chosen for comparison one representative from each of the two major types of trend removal methods parametric and nonparametricextensively reviewed in 5. The perfor- mance of the ACD algorithm is slightly better than that of the reference methods, but its main quality is the possibility to apply it automatically, without any subjective intervention. In the next section I describe the ACD algorithm. Then, in Sec. III I present the types of artificial time series used to analyze the ACD algorithm and an estimation of the magni- tude order of the standard deviation of a stationary noise superposed on a monotone trend. Section IV contains the analysis of the ACD parameters and the description of the automatic algorithm associated with the ACD method. The performance comparison of the ACD algorithm with the polynomial fit and the jump process is made in Sec. V, where I also analyze the variation of the global temperature anoma- lies over the last 1800 years. In Sec. VI nonmonotone trends and short time series are analyzed. In the last section I sum- marize the results and briefly discuss the extension of the ACD method to a nonmonotone trend. II. ACD METHOD Consider a time series x n , with 1 n N, generated by the stochastic process X n = f n + Z n , 1 where Z n is a discrete stochastic process with Z n =0 and f n are the values of the trend f tat the moments t n = n -1t, t being the sampling interval. The time takes values within the finite interval t 0,1. We assume that Z n does not depend on the trend values f n and f thas a slow and monotone variation and a vanishing temporal average f ¯ = N -1 n=1 N f n =0. If the stochastic process Z n is stationary, then a math- ematical justification of the ACD method can be given. We denote by p z zthe probability distribution of Z n and by p x x , nthat of X n , which depends explicitly on the time in- dex n because X n is a nonstationary process. According to Eq. 1, the two distributions are identical with the exception of a translation by f n , i.e., p x x , n= p z x - f n . The infinitesi- mal interval about a given real number is denoted by I = - /2, + /2. The average number of values of the time series x n lying within I is given by N , with N = n=1 N p z - f n . 2 We define the average one-step displacement with the ini- tial value in the neighborhood of g= 1 N n=1 N-1 p z - f n X n X n I , 3 where X n = X n+1 - X n . Using the conditional probability den- sity for successive values of the stationary stochastic process Z n denoted by p z zz, we can write X n X n I = - + x - p z x - f n+1 - f n dx , where relation p x x , n +1 , n= p z x - f n+1 - f n is used. From the simple change of variables z = x - f n+1 and from the definition of the conditional probability and the consistency condition for the joint probability density p z z, z p z - f n = - + p z z, - f n dz , it follows that Eq. 3becomes g= 1 N n=1 N-1 f n+1 - f n p z - f n + . 4 The first term is the average of the one-step variation of the trend within the neighborhood of . A similar relation is ob- tained if the final value is included in I . The second term in Eq. 4 = 1 N n=1 N-1 - + z + f n - p z z, - f n dz 5 measures the difference between gand the trend average slope. Let us investigate when this term vanishes. If the trend variation at one time step is much smaller than that of the noise  f t n t Z for all n, then the sum can be approxi- mated by an integral. If, in addition, the trend is linear f t = at + b and we make the change of variable = - f t, then we obtain = a N t -f 1 -f N d - + z - p z z, dz . Consider that the noise Z n has the symmetry property p z z = p z -zand p z z, z= p z -z,-z. It follows that for = f 1 + f N /2 the term vanishes and its value increases when approaches the extreme values of the time series. If the noise is moderately asymmetric, then vanishes for a different but close to f 1 + f N /2. Under the conditions of the previous paragraph, since f n+1 - f n = at, from Eqs. 4and 2it follows that g = at, i.e., we find the exact slope of the linear trend. In the general case of nonlinear trends, g/ t is only an approxi- mation of the trend slope f t. The trend Ftestimated by the ACD method is calculated from the requirement that its CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-2
derivative expressed with respect to the function values not to the argumentshould be proportional to the average one- step displacement g= FF -1 t . 6 This relation holds only if Ftis invertible, i.e., if gpre- serves the same sign over all its domain of definition. Our goal is to find the numerical quantity corresponding to the theoretical one defined in Eq. 3using the values of a time series x n obtained as a realization of the stochastic process X n . The infinitesimal intervals I must be replaced with finite ones. We divide the domain of the time series values into disjoint intervals I s = s , s+1 , s =1,2,..., S, so that any value x n is contained into an interval I s . Denote by N s the number of values x n lying within I s corresponding to the theoretical quantity N defined in Eq. 2. The intervals I s can be chosen in many ways. The simplest solution is to use homogeneous intervals, i.e., the values of N s should dif- fer from each other by a unit at the most. Then the total number of intervals S is the single parameter describing the distribution of the series values. For the final form of the ACD algorithm we use a distribution into nonhomogeneous intervals described in Appendix B. The one-step variation of the time series is x n = x n+1 - x n . For a given s, we compute the sample average of x n under the condition that the initial or final values should be in- cluded into the interval I s gˆ s = 1 2N s x n I s x n + x n+1 I s x n , 7 which corresponds to the definition in Eq. 3. If all the val- ues gˆ s have the same sign, then we can use them to deter- mine a numerical approximation of the trend by a piecewise linear curve denoted F ˆ tand corresponding to the theoreti- cal one in Eq. 6. Since the initial f 1 and the final f N values of the trend are unknown, we use instead the extreme values of the time series x n . The domain of definition t 0, Tof the function F ˆ tis different from that of the real trend f t. Therefore the nu- merical values of the estimated trend are obtained either by translation F ˆ n = F ˆ + nt, or by scaling F ˆ n = F ˆ nT / N. The optimal estimated trend is given by the requirement that the difference x n - F ˆ n should have a minimum standard deviation. If the fluctuations described by causes the quantities gˆ s to have different signs, then the estimated trend F ˆ n cannot be determined and the fluctuations are smoothed by means of a moving average. Therefore the numerical algorithm of trend removal using the ACD method consists of a succession of trend extractions and moving average smoothings. Denote by x n i the time series obtained after a succession of i extrac- tions and smoothings. Initially x n 0 = x n . First, at each step i, we try to remove an estimated monotone trend F ˆ n i x n i+1 = x n i - F ˆ n i . 8 If gˆ s does not have the same sign F ˆ n i =0 for all n, then x n i+1 is computed by moving averaging. After i steps, the estimated trend is the sum of the removed components F ˆ n = i F ˆ n i 9 and the estimated random component is zˆ n = x n - F ˆ n . 10 The time series processing is interrupted when the standard deviation of the residual x n i is times smaller than that of the initial series is a given positive real number, or if the standard deviation of x n i+1 is larger than that of the previ- ous step x n i , or if by adding the component F ˆ n i+1 the estimated trend F ˆ n is not monotone. For the interior values two-sided moving averages are performed over an interval of length 2K +1. If n K n N - K, then the average is taken over the first n + K the last N - n + K +1values. This asymmetric average forces the val- ues near the time series boundaries to follow the variations of the interior values. The initial averages are performed on small intervals such that the trend should be deformed as little as possible. For the first smoothing we use K =1. If the quantities gˆ s do not acquire the same sign, then K is gradu- ally increased by a unit for each new smoothing up to a maximum value K f and then the next smoothings are com- puted keeping the same value for K. In this way a compro- mise is made between the computing efficiency and the re- quirement that the trend should not be distorted when the noise is small. A numerical problem of the ACD algorithm described above is the possibility that some values gˆ s i could be very close to zero and then the segments of the estimated trend in the corresponding intervals I s would be almost parallel to the time axis. In such a case the estimated trend would be arti- ficially deformed and its length would be much longer than the length of the initial series T 1. In order to eliminate this possibility we impose the additional condition that the absolute value of the slope gˆ s i should have an inferior bound. Denote by n s min n s max the time step when the signal takes for the first time the last timea value in the interval I s . The difference n s = n s min - n s max is the number of time steps during which the time series takes all the values in I s . Then we choose the minimum value of the estimated trend slope as gˆ s i s+1 - s / n s . Since the ACD algorithm is rather complex, it is affected by many error sources. First, there is the theoretical error given by Eq. 5which is larger when the trend nonlinearity is greater, the time step is larger, and the coordinate is closer to the extreme values of the time series. Then, the numerical implementation induces other errors such as the replacement of the infinitesimal intervals I with the finite ones I s , the use of the extreme values of the time series instead of those of the real trend, and the trend distortion by moving average. In spite of the error diversity, the evaluation AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-3
presented in Secs. V and VI shows that the final error of the ACD algorithm is not larger than the errors of the usual trend estimation methods. This result could be explained by the fact that the errors of the successive estimations of the trend components F ˆ n i compensate each other. III. ARTIFICIAL TIME SERIES The ACD algorithm performance is analyzed using statis- tical ensembles of time series with as various as possible characteristics. We generate the time series according to Eq. 1superposing realizations of a given stochastic process Z n on a monotone trend. Since in Sec. V the ACD algo- rithm is compared with the polynomial fit, we do not analyze the algorithm for a polynomial trend, but for f t= t / a - t, t 0,1, with a 1. If the trend were of the same functional form as the function used in the least-squares-fit method, then the ACD method would be disadvantaged. We choose this functional form because its slope f t= a / a - t 2 has a nonhomogeneous distribution with its extreme values near the boundaries of the definition domain: the minimum at t =0 and the maximum at t =1. It is more difficult to remove such a trend than one with the maximum slope in the interior of the definition domain because numerical algorithms are more difficult to implement near the boundaries of the time series. In our case the ratio between the maximum and the minimum slope is f 1/ f 0= 1-1/ a -2 . We give to the parameter a values within the interval a 1.1,2.0, so that the ratio of the extreme values of the slope varies between 121 and 4. The resolution of the time series is varied with one order of magnitude choosing the number of values N within the interval N 100,1000. The majority of the numerical tests are performed on sig- nals with a stationary random component z n so that they obey the theoretical considerations in Sec. II. We numeri- cally generate realizations of a stationary stochastic process Z n of type AR1autoregressive of order onedefined by Z n = Z n-1 + G n , where 0 1 and G n is a Gaussian noise with null mean and standard deviation G . The properties of the AR1stochastic process are well-known 6. Its prob- ability distribution is Gaussian with null mean and standard deviation Z 2 = G 2 1- 2 -1 . When the value of the parameter increases, the successive values of the noise become more correlated. In the numerical tests takes values within the interval 0,0.99, white noise being obtained for =0. By means of the noise standard deviation Z we control the ratio between the variation due to the noise and that due to the trend. In order to generate signals of both types, the val- ues of Z are chosen within the interval Z 0.1,10. Fig- ure 1acontains two signals with AR1noise, one domi- nated by noise, the other by trend. Taking into account the intention that the extension of the ACD method for a nonmonotone trend should be used in DFA, I have chosen for preliminary tests a 1/ f noise, i.e., Z n Z n+h  h , where 0 1. Such stationary time series with unit standard deviation are generated applying the nu- merical method in 7for a spatially correlated noise on a one-dimensional time lattice. An important step in DFA is the polynomial trend removal from a series obtained by sum- ming up the initial one y n = k=1 n x k , which is a nonstationary fractional Brownian motion. Figure 1bshows two signals y n obtained by summing up a time series x n with a 1/ f noise with = 0.4. Although the theory presented in Sec. II does not hold for nonstationary noise, we remove the trend replacing in Eqs. 7–10x n and z n with y n and k=1 n z k , respectively. In time series processing a critical parameter is the ratio of the overall trend variation to the amplitude of noise fluc- tuations quantitatively expressed as 0 = f / Z , 11 where f is the trend standard deviation and Z is the noise standard deviation. Figure 2 shows the dependence of 0 on the noise standard deviation Z for the extreme values of the trend parameter a. The global variation of 0 ranges over three magnitude orders providing a sufficient diversity of the time series characteristics on which the ACD algorithm is tested. The signals in Fig. 1acorrespond to the two possible 0 200 400 600 800 1000 −4 −2 0 2 4 6 8 10 12 n x n (a) 0 200 400 600 800 1000 0 500 1000 1500 n x n (b) FIG. 1. Examples of signals used to analyze the performance of the ACD method. The noise is superposed on the trend f t= t / a - t, t 0,1represented with a thin line. The values of the trend parameter are a =2 thick lineand a = 1.1 point markers. aStationary AR1 noise with = 0.9 and Z =1. bNonstationary noise given by a realization of a fractional Brownian motion obtained by summing up a 1 / f noise with = 0.4. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-4
cases: 0 1 for a = 1.1 when the signal is dominated by trend and 0 1 for a = 2.0 when the signal is dominated by noise. Thus a time series with AR1noise is characterized by four parameters: the series length N related to the reso- lution of the series, the trend parameter a, the parameter describing the correlation between the successive terms of the noise, and the noise standard deviation Z . In building an automatic algorithm for time series pro- cessing it is useful to have a simple criterion to estimate the order of magnitude of the noise standard deviation Z using only the values of the original series x n . The random com- ponent z n is actually one of the unknown quantities which has to be determined by numerical processing. Denote by m x n = x m+n - x n the variation of order m N of the time series x n . According to the justification given in Appendix A, a workable estimation of the magnitude order of Z is given by Z est = 1 2 ˆ m 0 x n , 12 where ˆ ·is the sample standard deviation and m 0 is the smallest integer number m N /2 for which m 0 +1 x n m 0 x n , 13 where ·is the usual quadratic norm. I have tested the accu- racy of the estimation 12on 1000 time series with AR1 noise, each of them with parameters N, a, , and Z ran- domly chosen within their maximum range of variation. In Fig. 3 a remarkable correlation between the actual value of Z and the estimated one from Eq. 12is noticeable. Only for Z 1 there are some estimated values almost one order of magnitude greater than the real value. A percentage of 2.8% time series for which there is no value m 0 satisfying Eq. 13is missing from the figure. Such situations occur when the noise standard deviation is very small see Appen- dix Aand then we consider Z est =0. In the next section we describe a method to obtain a nonvanishing estimation even in these cases. In Appendix A we also prove that estimation 12is meaningful even if the trend is not monotone, such that it can be applied without additional restrictions on an arbitrary signal. IV. ACD PARAMETERS In order to apply the ACD method there are two things to be done: 1establish the values of two parameters the maximum value of the ratio between the initial signal stan- dard deviation and the final residual standard deviation and the maximum value K f of the semilength of the averaging interval; and 2distribute the signal values into disjoint intervals I s . In this section we show how these tasks can be automatically accomplished. First we analyze the influence of the parameter on the ACD performance. Since no criterion to choose the values of the other parameter N f has been established, its values are chosen randomly according to a homogeneous probability distribution over its variation range. The parameter N f con- trols the speed of the damping of the signal fluctuations by the moving average. Its value is calculated from the value of the semilength K f of the averaging interval N f =2K f +1 which takes values in the interval K f / N 0.01,0.1. For the shortest signals N = 100the minimum value of K f is K f =1. In this stage of the analysis we use homogeneous intervals I s , such that they are completely described by a single pa- rameter, i.e., the total number of intervals S. The minimum value of this parameter is S min =2, and the maximum one S max = N / N min , where ·is the integer part function and N min = 14 is the minimum number of the signal values neces- sary to recognize a Gaussian white noise, which is deter- mined in Appendix A. In the end we shall introduce a non- homogeneous distribution of the data values. The optimum value of the parameter is found by means of the estimation 12of the noise standard deviation. It is known that the sample mean of a random variable with a Gaussian distribution has the standard deviation N r times smaller than that of the random variable, where N r is the number of realizations. Considering that the reduction by 10 −1 10 0 10 1 10 −2 10 −1 10 0 10 1 10 2 σ Z η 0 a=1.1 a=2.0 FIG. 2. Average ratio of the trend standard deviation to the noise standard deviation in terms of standard deviation of the noise for the maximum value a = 1.1and the minimum one a = 2.0of the trend parameter for time series with AR1noise. 10 −1 10 0 10 1 10 −2 10 −1 10 0 10 1 10 2 σ Z σ Z est FIG. 3. Correlation between the standard deviation of an AR1 noise superposed on a monotone trend and the estimated standard deviation of the noise. AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-5
N r times of the standard deviation is a measure of the maxi- mum effectiveness which a statistical method can have, we assume that in the case of the ACD method the standard deviation of the final residual can be at best Z / N as well. Then we assume that the optimum value of is opt = Nˆ x Z est , 14 where ˆ x is the sample standard deviation of the original signal. As noticed in the previous section, there are situations when Eq. 12does not allow the determination of Z est and then opt cannot be calculated. Since such situations occur if Z is small with respect to the variation due to the trend, then first the trend is removed from the time series x n and the estimation 12is applied to the residual obtained. The trend removal is performed by means of the ACD algorithm with S = S max homogeneous intervals and K f = 0.1N. If the trend removal is not possible, then we apply Eq. 12to the differ- ence between the initial signal and its moving average performed with K f = 0.01N. In this way, we always obtain a value Z est by means of which opt in Eq. 14can be calculated. Figure 4 shows the dependence of the accuracy of the ACD method on the three ACD parameters. We evaluate the accuracy of the trend removal by means of the index = F ˆ n - f n /f n . 15 The smaller , the more alike the original and the estimated trends are. If the method ACD does not succeed to remove the trend F ˆ n =0, then = 1. Therefore 1 means that the estimated trend differs more from the real one than a vanish- ing trend and the estimated trend removal would cause a distortion of the information contained by the signal. In such cases the trend removal is not recommended 8and in Sec. V we describe a method to recognize such situations. The results in Fig. 4 are obtained for statistical ensembles of time series with lengths N = 100 and 1000 and with trend param- eters a = 1.1 and 2.0. The number of time series in a statistical ensemble has a value between 300 and 3000, such that for each type of time series the standard error for the average evaluation index is smaller than 5%. Figure 4ashows the average evaluation index with respect to the parameter for random values of the param- eters K f and S. The time series have the random component either stationary of type AR1with Z =1 and = 0.9 con- tinuous linesor a nonstationary fractional Brownian motion with = 1.4 dashed lines. For opt the average index reaches a stationary value in all cases, this kind of behavior showing that the value given by Eq. 14assures the retrieval of all available information from the original signal. As a rule, for opt , decreases for increasing , showing the 10 −2 10 −1 10 0 10 1 0.2 0.5 1 1.5 2 ρ /ρ opt 〈 θ 〉 (a) 10 −2 10 −1 10 0 10 1 0.02 0.05 0.1 0.5 1 ρ /ρ opt 〈θ〉 (b) 0.02 0.04 0.06 0.08 0.1 0.2 0.5 1 1.5 2 K f /N 〈 θ 〉 (c) 0 0.2 0.4 0.6 0.8 1 −0.5 0 0.5 1 1.5 r S 〈 ln(θ /θ ∗ ) 〉 (d) FIG. 4. Dependence of the average evaluation index of the ACD algorithm on the ACD parameters. The time series are characterized by different values of their length N and of the parameter a in the trend formula f t= t / a - t, t 0,1. The couple N , atakes the values 100,2.0asterisk marker , 100,1.1cross marker , 1000,2.0circle marker , and 1000,1.1square marker . In panels a, c, and dthe random component of the signals is a stationary AR1noise with Z =1 and = 0.9 continuous linesand a nonstationary fractional Brownian motion with = 1.4 dashed lines. In panel bthe random component is a stationary AR1noise with Z =0.1 and = 0.9 continuous lines, and =0 dashed lines. aand bThe variable parameter is the maximum value of the ratio between the standard deviation of the initial signal and that of the final residual, while the other two ACD parameters have random values. A log-log scale is used. cThe variable parameter is the maximum value K f of the semilength of the averaging interval while = opt and S takes random values. A log-linear scale is used. dThe evaluation index for S homogeneous intervals divided to the evaluation index * obtained using nonhomogeneous intervals in terms of the ratio r S = S - S min / S max - S min while = opt and K f = 0.1N. A standard error for smaller than 5% is obtained with statistical ensembles containing 500 a, 3000 b, and 300 cand d time series. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-6
importance of a right choice of opt , because if has too small a value, then the trend is not completely removed. The only situation in which, for opt , increases for increas- ing is for time series with N =100 and a = 2.0 asterisk marker. For such time series, for greater the trend removal leads to worse results, i.e., decreasing the standard deviation of the final residual, the estimated trend becomes more and more different from the real one. This behavior is due to the fact that when the signal is dominated by noise and its reso- lution is low, the estimated trend follows the noise fluctua- tions, not the trend. As seen in Fig. 4asignals with the same characteristics as the ones previously discussed but with N = 1000 and marked with open circles, increasing the signal resolution can render accessible the information on the trend. These results show that the optimum value of the pa- rameter is correctly given by Eq. 14. Equation 14implies a special procedure to calculate Z est when Z has small values. The validity of this procedure is verified in Fig. 4bfor time series with AR1noise with the minimum standard deviation Z = 0.1. In all cases the average evaluation index reaches a stationary value for = opt after significant decreases. This behavior is in accordance with the fact that when the noise is small, the trend can be easily removed. Another observation is that the results de- pend on the parameter , i.e., on the correlation between successive values of noise. For Z 1 the results are similar with those previously discussed for signals dominated by noise and we do not present them here. With the optimum value of the first parameter given by Eq. 14, we go on to analyze the parameter K f controlling efficiency of the moving averaging in the ACD method while the data values are distributed into a random number S of homogeneous intervals. The larger K f is, the longer the av- eraging intervals and the more strongly damped the noise fluctuations are, but at the same time the more distorted the trend is. Figure 4ccontains the results for signals with AR1noise with Z =1 and = 0.9 continuous linesand with nonstationary fractional Brownian motion with = 1.4 dashed lines. Over the entire variation range of K f the av- erage evaluation index has approximately the same value. Also for other values of Z and the behavior of is similar, therefore we could choose for K f an arbitrary value. However, the number of averagings decreases with almost two orders of magnitude when K f varies and to save com- puting time we use its maximum value K f = 0.1N. Finally we analyze the optimal distribution of the data values into intervals I s , using for and K f the optimum val- ues determined above. When we distribute the time series values into disjoint intervals I s , we have to take into account two opposite requests. If the noise fluctuations are much smaller than the trend variation, then it is recommended to use a large number of intervals I s in order to describe the trend shape as accurately as possible. On the contrary, if the noise fluctuations are much larger, then it is better to use a smaller number of intervals I s , each of them containing more signal values so that the noise fluctuations may be smoothed as much as possible in the average slope given by Eq. 7. Therefore we distribute the signal values into a number of disjoint intervals inversely proportional with the estimated standard deviation of the noise Z est given by Eq. 12, S est = x max - x min / Z est , 16 where x max and x min are, respectively, the maximum and the minimum value of the time series x n . If the value obtained from Eq. 16is smaller than S min larger than S max = N / N min , then we impose S est = S min S est = S max . Since, as a rule, N is not exactly divisible by S est , the distribution of the series values into S est homogeneous intervals is performed so that the number of values N s for different s may differ at most with a unit. The distribution of the signal values into S est homoge- neous intervals must be corrected in order to take into ac- count some numerical characteristics of the ACD algorithm. In Appendix B we describe an algorithm of splitting and merging the homogeneous intervals, in the end resulting S * nonhomogeneous intervals. The ACD method with nonho- mogeneous intervals cannot be evaluated as easy as that for homogeneous ones. The values of the series x n can be dis- tributed into a given number of intervals in many ways. In an exhaustive analysis, the algorithm described in Appendix B should be compared with all the possible distributions of the signal values for all the allowed values of S. However, we are satisfied with an automatic algorithm having an accuracy comparable with that of the simplest nonautomatic form of the ACD method. Therefore we use as a reference in our comparisons the ACD method with homogeneous intervals of signal values. Figure 4dpresents the average of the quantity ln/ * , where is the evaluation index in Eq. 15 obtained for different numbers S of homogeneous intervals and * is the evaluation index for S * nonhomogeneous inter- vals obtained according to Appendix B. Since the maximum value S max = N / N min of the parameter S depends on the sig- nal length N, in Fig. 4dthe ratio r S = S - S min / S max - S min is used as an independent variable. The result of a trend removal is better when the evaluation index is smaller and then the positive values of the quantity ln/ * indicate that the ACD algorithm with S * nonhomogeneous intervals gives a more correct estimated trend than with ho- mogeneous intervals. In Fig. 4done notices that in most of the presented cases the accuracy of the two forms of the ACD method is approxi- mately equal and in the rest of the cases the accuracy ob- tained with nonhomogeneous intervals is better. Only in a single case signals with 1 / f noise of length N = 1000 and trend parameter a = 1.1there are values of S for which the accuracy obtained with homogeneous intervals is signifi- cantly better. For all the analyzed signals with AR1noise results not presented herethe ACD algorithm with nonho- mogeneous intervals have comparable or better accuracy. From these results we can conclude that although the distri- bution of the signal values into S * nonhomogeneous intervals is not the best solution, however, its accuracy is good enough to be used in the automatic form of the ACD method. Since the verification of the optimum values of the param- eters and K f has been performed using the distribution of the signal values into homogeneous intervals, it is necessary to repeat the tests using nonhomogeneous intervals. The re- AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-7
sults are almost identical to those in Fig. 4 and we do not present them here. Thus we have obtained an automatic form of the ACD method for monotone trend removal theoretically and nu- merically substantiated for time series with stationary noise which can be useful as well if the noise is nonstationary. Figure 5 shows the average evaluation index for the final automatic form of the ACD algorithm with nonhomogeneous intervals applied to the artificial time series discussed in Sec. III, with different lengths N. In all the cases, when N in- creases, the accuracy of the trend removal is improving. For signals with large noise and small slope Z =10 and a = 2.0 for all the values of N,  10 showing that the estimated trend follows more the noise in the signal than the real trend. Increasing the trend slope a = 1.1or decreasing the noise Z =1, the evaluation index reaches the usefulness limit  1. For the rest of the time series, even for those with 1/ f noise, we obtain 1, indicating that they allow an efficient trend removal by means of the ACD algorithm. On a standard 3 GHz Pentium 4 PC processing a time series with N = 1000 required 0.4 CPU seconds for AR1 noise and 0.8 CPU seconds for 1/ f noise. The maximum number of partial estimated trend components was 4, and it was not significantly influenced by the series length or noise type. The number of the moving averagings increases di- rectly proportional to the series length. For N =1000, on av- erage, 80 smoothings were needed and the maximum number of smoothings was 140. V. EVALUATION OF THE ACD ALGORITHM The evaluation of the ACD method is performed by com- parison with two trend removal methods used at present, one parametrical, the other nonparametrical. The trend is re- moved from time series of the same type as those described in Sec. III. The comparison of the accuracy of different methods is realized by means of the quantity ln/ ACD , where is the evaluation index in Eq. 15of the method to be compared with and ACD is the evaluation index of the automatic ACD algorithm. We use the logarithm of the quan- tity / ACD because this way the cases when ACD are better visualized. The result of a trend removal is better when the evaluation index is smaller and then the positive values of the quantity ln/ ACD indicate that the ACD method gives a more accurate estimation of the trend than the refer- ence method. Figure 6 shows the results of the comparison with the polynomial fit of degrees from 1 to 10. First we notice that the polynomial fitting does not dispose of a parameter equivalent to of the ACD algorithm, i.e., the polynomial trend is estimated by a single step and there is no possibility to control the magnitude of the final residual. In addition let us remind that if the signal is dominated by noise, then, as the degree of the polynomial increases, the number of coef- ficients to be determined by least-squares fit increases too and they are more strongly affected by the random compo- nent of the signal. Therefore the accuracy of the polynomial fit for noise dominated signals decreases with the polynomial degree. Reversely, if the signal is dominated by trend varia- tion, a greater polynomial degree implies an improvement of the trend extraction results. Since ACD does not depend on the polynomial degree, the curves plotted in Fig. 6 describe the dependence of the quan- tity ln poly on the polynomial degree, but translated with a constant -ln ACD , different for each type of time series. Hence comparing Figs. 6band 6cone can see that, as expected, the polynomial fitting gives better results if the polynomial has a larger smallerdegree when the signal is dominated by trend noise. This contradictory behavior makes very difficult an automatic choice of the optimum value of the polynomial degree for all the possible types of time series. From the first three panels of Fig. 6, it is difficult to draw a definite conclusion regarding the relation between the ACD accuracy and that of the polynomial fit. For Z =1 continu- ous lines in Fig. 6a, the ACD algorithm has better perfor- mance for short time series N = 100, when the fitted poly- nomial follows more the noise fluctuations than the trend, especially for higher polynomial degrees. If the resolution is increased N = 1000, then for degrees higher than 3, the two methods have comparable accuracy. For the signals domi- nated by trend in Fig. 6b Z = 0.1the accuracies are com- parable for low resolution. At high resolution the results of the two methods are very sensitive to the polynomial degree. For degrees smaller than 5 the ACD accuracy is significantly better and for degrees higher than 5 significantly worse. For the signals dominated by noise in Fig. 6c Z =10, the ACD is better in almost all cases. Although the accuracy of the ACD algorithm is not significantly better than that of the polynomial fitting, this comparison indicates that the auto- matic form of the ACD algorithm has been obtained preserv- ing performance comparable with that of existing algorithms. For signals with 1 / f noise dashed lines in Fig. 6a, except the polynomials with degree 1 or 2, the ACD algo- 200 400 600 800 1000 10 −1 10 0 10 1 N 〈 θ 〉 FIG. 5. Dependence of the average evaluation index of the ACD algorithm on the time series length. The time series differ from each other by the parameter a in trend formula f t= t / a - t, t 0,1 continuous lines correspond to a = 2.0 and dashed lines to a = 1.1 and by the standard deviation of the AR1noise with = 0.9 marker corresponds to Z = 10, to Z =1, and to Z = 0.1. The signals with 1 / f noise are marked with . A standard error smaller than 5% is obtained with statistical ensembles containing 300 time series. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-8
rithm has worse results. An explanation of this behavior could be that many times it is difficult to differentiate a non- stationary noise containing a so-called “stochastic trend” from a variation due to a deterministic trend. For such sig- nals the ability of the ACD algorithm to adjust itself to the signal shape is a drawback since the estimated trend de- scribes not only the deterministic trend but also the nonsta- tionary noise. Finally, Fig. 6dshows the fraction of the time series for which the estimated polynomial trend is not monotone. Ob- viously, the polynomials of order 1 preserve the monotony, but for higher orders more than one-half of the estimated trends are not monotone. It is interesting to remark that the polynomials of odd degree allow a better identification of the trend monotony. The second comparison method is the “jump process” in 5and it consists in an iterative weighted averaging x n i+1 = x n i + Rx n-1 i -2x n i + x n+1 i , 17 where ratio R 0 R 0.5and iteration parameter M i Mare user-specified constants. At boundaries a sym- metric extension of data is employed. As in the case of the polynomial fitting, the averaging in Eq. 17has opposite effects on the accuracy of the estimated trend. When the effect of the averaging is greater larger values of the param- eters R and M, the noise fluctuations are more strongly smoothed, but at the same time the trend shape is more distorted. Figure 7 contains the results of the comparison of the ACD algorithm with the jump process for different values of the parameter R and for M = 100 fixed. For signals with AR1noise, excepting a few cases with Z = 0.1 in Fig. 7b, the ACD method has better accuracy. For signals with 1/ f dashed lines in Fig. 7a, although not so bad as for the polynomial fit, the results of the ACD algorithm continue to be worse than those of the comparison method. As discussed above, comparing Figs. 7band 7cone notices that, as expected, the jump process has better results if the averaging parameter R is smaller largerwhen the signal is dominated by trend noise. From Fig. 7d, which shows the fraction of time series for which the estimated trend is not monotone, it follows that when R increases, more and more estimated trends preserve the monotony of the real trend. The fraction of the nonmonotone estimated trends depends also on the standard deviation of the noise and on the signal length and it becomes smaller if we increase the number of iterations M. From the comparisons discussed above, we conclude that, on average, the ACD algorithm removes monotone trends with slightly better accuracy than the other two methods if the noise is stationary. In the case of nonstationary noise the accuracy is worse. Nevertheless, the ACD algorithm pre- serves the advantage of its automatic form. In addition, the ACD algorithm removes only monotone trend, not non- monotone trend as the existing algorithms do. To illustrate these special characteristics we estimate by means of the three algorithms the monotone trend of a real time series from paleoclimatology. Figure 8ashows the global mean annual temperature anomalies during the period A.D. 200– 2 4 6 8 10 −0.5 0 0.5 1 1.5 q 〈 ln(θ poly /θ ACD ) 〉 (a) 2 4 6 8 10 −2 0 2 4 q 〈 ln(θ poly /θ ACD ) 〉 (b) 2 4 6 8 10 −0.5 0 0.5 1 1.5 q 〈 ln(θ poly /θ ACD ) 〉 (c) 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 q nonmonotone fraction (d) FIG. 6. Comparison of the ACD accuracy with that of the polynomial fit in terms of the polynomial degree q. The time series are characterized by different values of their length N and of the parameter a in the trend formula f t= t / a - t, t 0,1. The couple N , a takes the values 100,2.0asterisk marker , 100,1.1cross marker , 1000,2.0circle marker , and 1000,1.1square marker . aThe random component is a stationary AR1noise with Z =1 and = 0.9 continuous linesand a nonstationary fractional Brownian motion with = 1.4 dashed lines. bThe random component is a stationary AR1noise with Z =0.1 and = 0.9 continuous linesand =0 dashed lines. cThe random component is a stationary AR1noise with Z =10 and = 0.9 continuous linesand =0 dashed lines. d. The random component is as in panel abut the represented quantity is the fraction of the time series for which the polynomial estimated trend is not monotone. A standard error for poly smaller than 5% is obtained with statistical ensembles containing 200 time series. AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-9
1995 with respect to the Northern Hemisphere mean annual temperature over 1856–1980 discussed in 9and freely ac- cessible from 10. The temperature anomalies series contains N = 1796 val- ues, many of them repeating themselves. In fact there are only 48 distinct values and the distribution of the time series values into disjoint intervals used in the ACD algorithm de- mands that all values should be distinct. To satisfy this re- quest without distorting the initial signal, we superpose on the original values numbers randomly generated with a ho- mogeneous probability distribution on a range 1000 times smaller than the minimum difference between two distinct signal values. The ACD algorithm has estimated the trend in a single removal performed after 129 averagings. The series values are distributed into S * = 31 nonhomogeneous intervals obtained by splitting S est = 9 homogeneous intervals. The es- timated standard deviation of the noise is Z est = 0.056 in com- parison with the standard deviation of the initial series ˆ x = 0.063. In Fig. 8bthe ACD estimated trend is compared with those estimated by the other two methods. As it follows from Fig. 6d, most of the estimated trends by nonlinear polynomial fit are nonmonotone. The linear trend, the only 0 0.1 0.2 0.3 0.4 −1 0 1 2 R 〈 ln(θ jump /θ ACD ) 〉 (a) 0 0.1 0.2 0.3 0.4 −1 0 1 2 R 〈 ln(θ jump /θ ACD ) 〉 (b) 0 0.1 0.2 0.3 0.4 −1 0 1 2 3 4 R 〈 ln(θ jump /θ ACD ) 〉 (c) 0 0.1 0.2 0.3 0.4 0 0.2 0.4 0.6 0.8 1 R nonmonotone fraction (d) FIG. 7. Comparison of the ACD accuracy with that of the jump process in terms of the averaging parameter R for a fixed number of iterations M =100. The time series are characterized by different values of their length N and of the parameter a in the trend formula f t = t / a - t, t 0,1. The couple N , atakes the values 100,2.0asterisk marker , 100,1.1cross marker , 1000,2.0circle marker , and 1000,1.1square marker . aThe random component is a stationary AR1noise with Z =1 and = 0.9 continuous linesand a nonstationary fractional Brownian motion with = 1.4 dashed lines. bThe random component is a stationary AR1noise with Z =0.1 and = 0.9 continuous linesand =0 dashed lines. cThe random component is a stationary AR1noise with Z =10 and = 0.9 continuous linesand =0 dashed lines. d. The random component is as in panel abut the represented quantity is the fraction of the time series for which the estimated trend is not monotone. A standard error for jump smaller than 5% is obtained with statistical ensembles containing 200 time series. 500 1000 1500 2000 −0.4 −0.2 0 0.2 t (years) temperature anomalies (°C) (a) 500 1000 1500 2000 −0.24 −0.23 −0.22 −0.21 −0.2 t (years) temperature anomalies (°C) (b) ACD linear jump FIG. 8. Global mean annual temperature anomalies during the period A.D. 200–1995 with respect to the Northern Hemisphere mean annual temperature over 1856–1980 9and the trend estimated by means of ACD algorithm continuous line. In panel bthe trends estimated by means of the linear fit dashed lineand jump process dotted lineare also presented. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-10
polynomial trend which is assuredly monotone, exaggerates the global variation of the time series and does not provide any information on the slope variation over the analyzed interval. In the case of the jump process, the time series must be averaged until all the nonmonotone variations of the esti- mated trend are eliminated. For R = 0.4, M = 329 637 averag- ings were needed such that the computing time became pro- hibitive. Besides that, even if the estimated trend has the same shape as the ACD trend, due to the large number of averagings, the magnitude of the global temperature varia- tion is underestimated. As shown in 8, it is possible that, when the noise is large, the trend removal should not be recommended. The series in Fig. 8ais such a case. In order to verify that the estimated trend F ˆ n is significant, we build surrogate series xˆ n = F ˆ n + ˆ n . Since the serial dependence of the noise has a significant influence on the estimated trend quality, we choose the surrogate series ˆ n as realizations of an AR1 stochastic process. We note that the white noise is obtained for = 0, so it is implicitly included among the possible sur- rogate series. From the standard deviation and the autocorre- lation at one time step of the estimated noise in Eq. 10we compute the value of the parameter used in surrogates generation. The fraction of the estimated trends for the sur- rogate series xˆ n having the same sign of the global variation as the estimated trend F ˆ n for the original series is a mea- sure of the probability that the estimated trend should corre- spond to a real one in the original series. For example, from the estimated trends for 100 surrogate series of the time se- ries in Fig. 8a, only 58 are decreasing showing that the association of a trend to the initial climatological series is hazardous. Climatologists are interested in the time periods with a monotone temperature variation associated with geophysical processes of global scale, as, for example, the global warm- ing in the last century. I have applied the surrogate series method described above on the global temperature anomalies over time intervals measured in centuries, i.e., intervals t 1 , t 2 with t 1 = 200,300, ... ,1900, t 2 = 300,400, ... ,1900,1995 and t 2 t 1 . The results are presented in Fig. 9. There are only three time periods to which we can associate a monotone trend with a confidence level of 99%. Two of them are re- lated to the global warming over the last 2 centuries and the third one corresponds to the 10th century. For 19 time peri- ods the ACD algorithm did not succeed in associating a monotone trend to the temperature variation. The rest of the estimated trends cannot be considered “significant,” a con- clusion which coincides with that in 9. VI. NONMONOTONE TREND So far we have analyzed only time series with a monotone trend. In this section we consider time series without trend or with nonmonotone trend. In these cases, as well, the ACD evaluated trend is monotone and it may or may not corre- spond to a monotone component of the actual trend. There- fore an evaluation of the statistical significance of the esti- mated trend is needed, either using the surrogate series method presented in the previous section, or applying a sta- tistical test for monotone trend detection. In the following we compare the results of the ACD algorithm with the Mann- Kendall test 11. The Mann-Kendall test is based on the quantity S = P - M, where P is the number of the pairs x n x m with n m and M is the number of pairs x n x m with n m. If x n are independent observations, then for N 10 the random vari- able Z MK = S -1/ S if S 0 0 if S =0 S +1/ S if S 0 with S = NN -12N +5/18 1/2 follows a standard normal distribution. The null hypothesis that there is no monotone trend is rejected when the computed Z MK value is greater in absolute value than the critical value z /2 , where is the chosen significance level. This test was extended to a serial dependent time series 12, but here we restrain ourselves to the simple original form of the test. The Mann-Kendall test is used especially in hydrology, other environmental sciences, and econometrics for short time series, even containing only several dozens of values. In order to obtain a complete analysis of the ACD algorithm performance, we process time series of lengths down to the inferior limit of validity of the Mann-Kendall test, i.e., N =10. To make the ACD algorithm applicable to such time series, we modify it as follows: for N 2N min = 28, instead of nonhomogeneous intervals we use S = 2 homogeneous inter- vals to distribute the signal values. According to Appendix A, in such cases the estimated noise standard deviation has a 200 400 600 800 1000 1200 1400 1600 1800 2000 400 600 800 1000 1200 1400 1600 1800 1997 t 1 (years) t 2 (years) − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − + − − − − − − − + + − − − − − − − + + + − − − − − − − + + + − − − − − − − + + + − − − − − − − − + + + + + − − − − − − − − + + + + + + − − − − − − + + + + + + + − − − − + + + + + + + − + + − − + − − − + − + − − + + − − + − + + − − − − + FIG. 9. Time periods t 1 , t 2 over which the global average tem- perature has an increasing signor decreasing signtrend. The three bigger signs correspond to significant trends with a con- fidence level of 99%. The ACD method did not succeed to estimate a monotone trend for the time periods to which no sign is attached. AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-11
larger error, however, the ACD algorithm can be applied without any further modifications. The so-called type I error of a statistical test occurs when the null hypothesis is rejected although it is true. Figure 10 shows the probability of type I error occurrence for the Mann-Kendall test with 5% significance level applied to AR1noise without trend, in terms of the correlation param- eter , for different lengths of the time series. The figure also shows the relative frequency of the monotone trends esti- mated by means of the ACD algorithm for the same time series. One notices that the ACD algorithm estimates incor- rectly a monotone trend much more frequently than it could occur as a random event. The error of the ACD algorithm diminishes when the time series length increases, however, it remains much greater than the probability of type I error occurrence for the Mann-Kendall test. All the trend estimation methods present errors due to the lack of discrimination between monotone and nonmonotone trends. As shown in Fig. 6dfor nonlinear polynomial fit and Fig. 7dfor the jump process, these methods estimate many times a nonmonotone trend when the actual one is monotone. Being designed to remove the monotone component of the trend, the ACD algorithm has an opposite behavior. The type II error of a statistical test occurs when the null hypothesis is accepted although it is false. Our null hypoth- esis is false if the signal contains a monotone trend. To ana- lyze this kind of situation we introduce a trend composed of a monotone linear trend and a sinusoid f t= pt + a sin 2t, t 0,1, where p and a are positive real parameters. This formula allows us to make a continuous transition from a monotone to a nonmonotone trend. When a = a 0 =1/ 2 0.159, the function f thas an inflection point at t = 0.5 with the tangent parallel to the Ox axis. Therefore if a a 0 the trend is monotone, and if a a 0 it is nonmonotone. The power of a statistical test is defined as one minus the probability of type II error. Figure 11 shows the power of the Mann-Kendall test for time series with Z =1 and the mono- tone trend described in the previous paragraph with a = a 0 and different values of the parameter p. As expected, the closest shape to the ideal test power power vanishing for p = 0 and rapidly increasing up to 1is obtained for white noise =0and longer series N =50. For correlated noise the results are worse because the probability of type I error is much larger than the significance level. Like in Fig. 10, the relative frequency of the monotone trends estimated by the ACD algorithm for the same time series is much greater than the Mann-Kendall test power, reinforcing the previous con- clusion that the ACD algorithm estimates the monotone com- ponent of the trend. In Fig. 12 we give an example of a monotone component estimated by the ACD algorithm for a nonmonotone trend. We have used a signal with N = 1000 values given by the nonmonotone trend f t= t +2a 0 sin 2t, t 0,1, without superposed noise. One notices that the estimated trend is significantly different from the “real” monotone component of the trend, which is a straight line. This example shows that for a nonmonotone trend one cannot define unequivocally a monotone component and the ACD monotone estimated trend is only one of the infinitely existing possibilities. To evaluate the ACD algorithm results for nonmonotone trend, the same method as that in Sec. V has been used. The time series have N =10 and, respectively, N = 50 values, in order to provide information on the ACD performance for very short time series too. The trend is given by f t= t + a sin 2t with variable a, such that the trend is continu- ously modified from a linear one a =0to the nonmonotone one represented in Fig. 12 a =2a 0 . Figure 13ashows the dependence of the average evaluation index of the ACD al- 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 N=10 N=100 N=1000 φ α FIG. 10. The probability of the type I error for the Mann- Kendall test with 5% significance level applied to AR1noise with different serial dependence dashed line. The continuous line rep- resents the relative frequency of the monotone estimated trends by the ACD algorithm for the same time series. For the results of the ACD algorithm a standard error smaller than 4% is obtained with statistical ensembles containing 1000 time series. 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 p power FIG. 11. The power of the Mann-Kendall test with 5% signifi- cance level dashed linesfor time series with different lengths N obtained by superposing AR1noise with Z =1 on the monotone trend f t= pt + a 0 sin 2t, t 0,1, in terms of the parameter p. The continuous line represents the relative frequency of the mono- tone estimated trends by the ACD algorithm for the same time series. The couple N , takes the values 10,0asterisk marker , 10,0.9cross marker , 50,0circle marker , and 50,0.9 square marker . For the results of the ACD algorithm a standard error smaller than 2% is obtained with statistical ensembles contain- ing 1000 time series. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-12
gorithm on the parameter a. From these results the same conclusions as those in Sec. IV in connection with Fig. 5 can be drawn: the accuracy is improving when N increases or decreases and the trend removal is efficient if the signal is not dominated by noise. Furthermore, the ACD accuracy is worsening when the contribution of the harmonic component increases. Since, as seen in Sec. V, the performance of the polyno- mial fit increases for smaller degrees, we compare in Fig. 13bthe ACD accuracy only with the linear fit. For small values of parameter a the trend is almost linear, such that the accuracy of the linear fit is better in all the cases. For greater values of the parameter a and for small noise, the ACD accuracy becomes better N =50or at least comparable N =10. Only for a signal dominated by white noise =0 and Z =1the results of the linear fit are better in all the cases. These are the only time series for which the jump process has a better accuracy than the ACD algorithm see Fig. 13c. So we can conclude that, also for very short time series or nonmonotone trend, the ACD accuracy remains comparable with that of the usual trend estimation methods. Besides, this result is obtained preserving the automatic fea- ture of the ACD algorithm. VII. CONCLUSIONS The automatic ACD algorithm for monotone trend re- moval presented in this paper can be applied on a time series with stationary noise without any preliminary subjective in- spection of the time series. The analysis of its performance on artificial signals with AR1noise indicates that the quality of the estimated trend is slightly better than that of the usual nonautomatic methods. This conclusion is confirmed by the results obtained for a real time series from paleoclimatology. The preliminary tests for nonstationary 1 / f noise show that the ACD accuracy is worse than that of the polynomial fit and likely of the parametric methods in general. The ex- planation is presumably related to the fact that the realiza- tions of a nonstationary stochastic process are difficult to differentiate from a signal containing a deterministic trend component. Then the ability of the ACD method to adjust itself to the signal shape makes possible the confusion be- tween the stochastic trend of the nonstationary noise and the deterministic trend, whereas the parametrical methods which have less degrees of freedom are less affected by such con- fusions. Therefore the application of the ACD method to nonstationary noises needs a special analysis. However, the possibility to apply automatically the ACD method to a large number of time series is an advantage which may justify its utilization in some applications like DFA. The analysis of the automatic ACD algorithm presented in this paper should be extended. Further tests should be per- formed both on different real time series and artificial ones with greater lengths N 1000, other types of stationary noises, and other functional forms of the monotone trend. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 t n x n t n + 2a 0 sin(2π t n ) ↓ FIG. 12. The estimated ACD trend dashed linefor a time series continuous linewith N = 1000 values, without noise, ob- tained by superposing a sinusoid on a linear component dotted line, using dimensionless time t n = n / N. 0 0.5 1 1.5 2 0 0.5 1 1.5 a /a 0 〈θ〉 (a) 0 0.5 1 1.5 2 −1 −0.5 0 0.5 a /a 0 〈 ln(θ poly /θ ACD ) 〉 (b) 0 0.5 1 1.5 2 −0.5 0 0.5 1 1.5 2 a /a 0 〈 ln(θ jump /θ ACD ) 〉 (c) FIG. 13. aDependence of the average evaluation index of the ACD algorithm on the parameter a in the trend formula f t= t + a sin 2t, t 0,1. The time series are characterized by different values of their length N and the random component of the signals is a stationary AR1noise with standard deviation Z and different serial dependence = 0.9 continuous linesand =0 dashed lines. The couple N , takes the values 10,0.1asterisk marker , 10,1cross marker , 50,0.1circle marker , and 50,1square marker . A standard error smaller than 2% is obtained with statistical ensembles containing 1000 time series. bComparison of the ACD accuracy with that of the linear fit for the same time series as for panel a. cComparison of the ACD accuracy with that of the jump process with R =0.4 and M =100 for the same time series as for panel a. AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-13
Also it is necessary to compare the performance of the ACD algorithm with other trend removal methods besides the two algorithms presented here polynomial fit and jump process, in comparison with which the ACD accuracy was proved to be slightly better. The extension of the ACD method to nonmonotone trends can be achieved by applying the ACD method separately to the increasing, respectively, decreasing part of the signal. This is possible if the sample average in Eq. 7is computed separately for the positive, respectively, negative variations of the time series values. In this case a unique trend cannot be determined because the final result depends on the manner in which the two monotone parts of the estimated trend are combined. However, the main advantage of the ACD method is preserved, i.e., the functional form of the estimated trend need not be imposed and the design of an automatic algo- rithm for nonmonotone trend removal is possible. In essence the ACD method tries to use in a more com- plete way the information available in the transition prob- ability of a stochastic process. The most elaborate method of this type is that in nonequilibrium statistical physics, where a hydrodynamical description of a fluid is obtained by averag- ing the two-particle distribution function. In 13it was shown that a description of hydrodynamical type can be ob- tained for any corpuscular system using space-time coarse- grained averages of a kinematic description of the movement of the component particles. This approach was applied to a one-dimensional system 14and then to a time series with the space coordinate replaced by the values of a share price 15. From this point of view, the ACD method is a dis- cretized form of the space-time coarse-grained averages with disjoint space averaging intervals and the time averaging in- terval identical with the signal duration. It is possible to de- velop the ACD method for a time series with nonmonotone trends as a more complete hydrodynamical description if multiple overlapping subintervals are used for both space and time. APPENDIX A: ESTIMATION OF NOISE STANDARD DEVIATION A justification not a rigorous proofof estimation 12of the noise standard deviation for a realization x n of the sto- chastic process in Eq. 1is based on the fact that Z n does not depend on the trend values f n . Then for N large enough we have the approximate equality m x n 2  m f n 2 + m z n 2 , A1 where ·is the usual quadratic norm and, as in Sec. III, we use the notation m x = x m+n - x n . Using the property f ¯ =0, one can prove that if f n is monotone, then for m N /2 the sequence m f n is increasing m+1 f n m f n . A2 It follows that the left-hand term in Eq. A1decreases when m increases only due to the term m z n 2 . According to Eq. 13, m 0 is the number of time steps when the first decrease of m x n 2 takes place, i.e., the value of m for which the variations of x n are dominated by noise. Since we intend to estimate only the magnitude order of Z , we neglect in Eq. A1the term due to the trend. If m 0 is large enough to render Z m 0 +n and Z n independent, then Eq. A1becomes m 0 x n 2 2z n 2 . A3 Since Z n = 0, the right-hand term is proportional with the sample standard deviation of the noise ˆ Z . For the left-hand term it is not rigorously true that it is proportional to the corresponding standard deviation, but taking into account the approximations made until now we write Eq. A3as estima- tion 12. The numerical tests presented in Fig. 3 show that the obtained estimation is remarkable for an approximation only of the magnitude order of the noise standard deviation. The estimation quality remains the same for other functional forms of the monotone trend polynomial, exponential, loga- rithmic, etc.. If we want to apply estimation 12on real signals, we have to establish what information on the noise standard de- viation is supplied by this formula if the trend is not mono- tone. Then Eq. A2is not true and the negative variations of m x n 2 with respect to m can result from both right-hand terms in Eq. A1. So, even if the signal does not contain noise, estimation 12can take for noise the trend oscilla- tions. Figure 14 shows the results obtained applying Eq. 12 to a sinusoidal signal f t= sin 2t, t 0,1. For 0.25, the signal is monotone and coincides with the trend, so that the noise term in Eq. A1vanishes and according to Eq. A2 m x n is increasing. Then there is no m 0 satisfying Eq. 13and it is normal to consider that Z est =0. In Fig. 14 Z est =0 not only for 0.25, but also for 1 = 0.4. Hence according to evaluation 12, a sinusoid containing up to 1 periods is taken for a trend without noise, even if it is not monotone anymore. Consequently, the estimation 12is able to recognize nonmonotone trends as well. For 2 =1, the estimated standard deviation oscillates about an asymptotic 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 1.2 ν σ Z est σ Z σ ∞ σ Z =0 σ Z =0.1 FIG. 14. Estimated standard deviation Z est of a Gaussian white noise superposed on a sinusoidal signal with frequency , f t = sin 2t, t 0,1, discretized on N = 200 time steps. The standard deviation of the noise is Z =0 circle markersand Z = 0.1 points markers. The horizontal lines represent the nonvanishing value of the noise standard deviation and the standard deviation of an infinite sinusoid =1/ 2. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-14
value est = 1, a larger value than the standard deviation of an infinite sinusoid =1/ 2. In other words, a sinusoid con- taining more than 2 periods is erroneously identified to a noise with standard deviation est . The interval 1 , 2 represents a transition domain between the two dominant behaviors of evaluation 12. Actually we are interested first of all in the ability of estimation 12to identify a noise superposed on a nonmono- tone trend. Therefore, on the sinusoid analyzed above we superpose a white Gaussian noise with Z =0.1. In this case the results strongly depend on the signal resolution. Figure 14 shows the estimations for N = 200 and one notices that for different realizations of the noise, different values of Z est for the same can be obtained. However, if 1 , all values of Z est are clustered about the real value Z , and if 2 , most values cluster about the value est obtained for sinusoid with- out noise. Hence if the sinusoid without noise is interpreted as a trend, then the noise standard deviation is correctly es- timated and when it is interpreted as noise then the whole signal is considered as noise. If the signal resolution in- creases N 200then the correct value of Z is obtained for 2 too. For example, if N = 1000, then Z est Z for 2. This result proves that Eq. 12gives a useful estimation of the order of magnitude of a noise superposed on a non- monotone trend. A detailed analysis of the extension for non- monotone trends of estimation 12will be the subject of a future paper. As a first application of Eq. 12we estimate the mini- mum number of values N min needed lest a Gaussian white noise should be interpreted as a trend. If the number of val- ues is very small N 10, then there is a significant prob- ability that no m 0 satisfying Eq. 13should exist and the noise is interpreted as a trend. For example, for N =4 the minimum value of N for which Eq. 13is meaningful, out of 1000 realizations of a Gaussian white noise with unit stan- dard deviation 35% are not recognized as noise. The largest value of N for which we still obtain Z est =0 in 0.1% cases is N =12. However, we choose N min =14 in order to be sure that in no case Z est =0, and then Z est = 1.105 with Z est 0.38,2.26. These results are consistent with the Mann- Kendall test for monotone trends 11, which is applicable only if N 10. The estimation of the noise standard deviation is strongly improved when N increases. For N =1000 we ob- tain Z est = 1.015 with Z est 0.935 , 1.110. These results are used in Sec. IV and in Appendix B to distribute the time series values into disjoint intervals. In the case of correlated AR1noise, the value of m 0 increases when approaches 1. For example, for 1000 time series with = 0.9, Z =1 and N = 1000, the average value of m 0 is 26.9 and Z est = 0.984 with Z est 0.724 , 1.365. The average value Z est is as accurate as that for the Gaussian white noise, but the fluctuations of the individual values are greater. This behavior is due to the fact that for large values of the stochastic trend becomes important enough to be interpreted in Eq. 12as a deterministic trend. APPENDIX B: NONHOMOGENEOUS INTERVALS OF SIGNAL VALUES The signals dominated by noise or trend impose different numerical restrictions on the ACD algorithm and it is advis- able that the numerical algorithm should adapt itself to the signal type. Since the parameter 0 defined by Eq. 11can be calculated only if we know the signal components, we have to introduce a new parameter expressed by means of the values of the initial series est = ˆ x Z est -1, B1 where ˆ x is the sample standard deviation of the signal and Z est is given by estimation 12. If est 1, then ˆ x - Z est Z est and the noise standard deviation represents the largest part of the whole signal, i.e., the signal is dominated by noise. Conversely, if est 1, then the signal is dominated by trend. Equation 16applied to a signal dominated by noise est 1generates a small number of homogeneous inter- vals I s , each of them containing a large number of values. For longer signals, the number of values in I s increases and the computing time can become prohibitive. For this reason, when est 1, we split into two subintervals those intervals which by splitting generate a larger number of values than the quantity N min defined in Appendix A. For example, con- sider the interval I s = s , s+1 . By splitting it, we obtain two subintervals I s = s , s+1 and I s+1 = s+1 , s+1 , where s+1 = s + s+1 /2. The two subintervals contain in general differ- ent numbers of values N s N s+1 . If N s N min and N s+1 N min , then the splitting is validated and the number of in- tervals I s increases by one. If all S est intervals can be split generating 2S est new intervals, then the splitting process is repeated until at least one of the intervals does not satisfy anymore the condition that the two subintervals should con- tain more than N min values. In the end we obtain S * nonho- mogeneous intervales I s containing arbitrary numbers of time series values. Applied to a signal dominated by trend est 1, Eq. 16 generates a large number of homogeneous intervals I s , each of them with a small number of values inducing large fluc- tuations of the mean slope in Eq. 7. Therefore we merge the intervals with smaller length than Z est , s+1 - s Z est . Con- cretely, the successive intervals I s , I s+1 ,..., I s+m are joined in a new interval I s = s , s+m+1 , containing N s = N s + N s+1 + ¯ + N s+m values if s+m+1 - s Z est and s+m+2 - s Z est . Thus we obtain S * nonhomogeneous intervals. AUTOMATICALGORITHM FOR MONOTONE TREND REMOVAL PHYSICAL REVIEW E 75, 036705 2007 036705-15
1G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forcast- ing and Control, 2nd ed. Holden-Day, San Francisco, 1976. 2D. B. Percival and A. T. Walden, Wavelet Methods for Time Series Analysis Cambridge University Press, Cambridge, En- gland, 2000. 3C. K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stan- ley, and A. L. Goldberger, Phys. Rev. E 49, 1685 1994. 4K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E. Stanley, Phys. Rev. E 64, 011114 2001. 5S. Zhao and G. W. Wei, Comput. Stat. Data Anal. 42, 219 2003. 6P. J. Brockwell and R. A. Davies, Time Series: Theory and Methods, 2nd ed. Springer-Verlag, New York, 1991. 7N.-N. Pang, Y.-K. Yu, and T. Halpin-Healy, Phys. Rev. E 52, 3224 1995. 8E. L. Andreas and G. Trevino, J. Atmos. Ocean. Technol. 14, 554 1997. 9P. D. Jones and M. E. Mann, Rev. Geophys. 42,1 2004. 10P. D. Jones and M. E. Mann, Climate over Past Millenia, Data Contribution Series # 2004-085 2004, in NOAA/NGDC Pa- leoclimatology Program, URL http://www.ncdc.noaa.gov/ paleo/pubs/jones2004/ 11M. G. Kendall, Rank Correlation Methods, 4th ed. Griffin, London, 1975. 12R. M. Hirsch and J. R. Slack, Water Resour. Res. 20, 727 1984. 13C. Vamoş, A. Georgescu, N. Suciu, and I. Turcu, Physica A 227, 81 1996. 14C. Vamoş, N. Suciu, and A. Georgescu, Phys. Rev. E 55, 6277 1997. 15C. Vamoş, N. Suciu, and W. Blaj, Physica A 287, 461 2000. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-16