Automatic algorithm for monotone trend removal

7 years ago

(original), ISI/JCR, Numerical Modeling, paper, Time Series

Abstract

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 evaluation 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.

Authors

C. Vamoș
-Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Keywords

trend removal; monotone trends; time series;

Cite this paper as:

C. Vamoş, Automatic algorithm for monotone trend removal, Physical Review E 75 (2007) article id. 036705, DOI: 10.1103/PhysRevE.75.036705

References

see the expansion block below.

PDF

2007 Vamos Automatic algorithm for monotone trend removal

About this paper

Journal

Physical Review E

Publisher Name

APS

DOI

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.75.036705

Print ISSN

2470-0053

Online ISSN

2470-0045

MR

ZBL

Google Scholar citations

citations

[1] G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forcasting and Control, 2nd ed. Holden-Day, San Francisco, 1976.

[2] D. B. Percival and A. T. Walden, Wavelet Methods for Time Series Analysis, Cambridge University Press, Cambridge, England, 2000.

[3] C. K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger, Phys. Rev. E 49, 1685 1994.

[4] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E. Stanley, Phys. Rev. E 64, 011114 2001.

[5] S. Zhao and G. W. Wei, Comput. Stat. Data Anal. 42, 219 2003.

[6] P. J. Brockwell and R. A. Davies, Time Series: Theory and Methods, 2nd ed. Springer-Verlag, New York, 1991.

[7] N.-N. Pang, Y.-K. Yu, and T. Halpin-Healy, Phys. Rev. E 52, 3224 1995.

[8] E. L. Andreas and G. Trevino, J. Atmos. Ocean. Technol. 14,

[9] P. D. Jones and M. E. Mann, Rev. Geophys. 42, 1 2004.

[10] P. D. Jones and M. E. Mann, Climate over Past Millenia, Data Contribution Series # 2004-085 2004, in NOAA/NGDC Paleoclimatology Program, URL http://www.ncdc.noaa.gov/ paleo/pubs/jones2004/

[11] M. G. Kendall, Rank Correlation Methods, 4th ed. Griffin, London, 1975.

[12] R. M. Hirsch and J. R. Slack, Water Resour. Res. 20, 727 1984.

[13] C. Vamoş, A. Georgescu, N. Suciu, and I. Turcu, Physica A 227, 81 1996.

[14] C. Vamoş, N. Suciu, and A. Georgescu, Phys. Rev. E 55, 6277 1997.

[15] C. Vamoş, N. Suciu, and W. Blaj, Physica A 287, 461 2000.

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 artiﬁcial 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 ﬁtting 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 ﬁtting has the best accuracy. As a veriﬁcation 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 brieﬂy discussed. DOI: 10.1103/PhysRevE.75.036705 PACS numbers: 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 ﬁt of a parametric family of functions or smoothing by moving average to more complex as discrete differencing 1 or wavelet analysis 2. More sophisticated methods should be used only if their results are signiﬁcantly 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 signiﬁcant, 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 signiﬁcant 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 artiﬁcial sig- nals with stationary noises of type AR1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 ﬂuctuation 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/753/03670516 ©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 signiﬁcantly 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 ﬁrst one is the polynomial ﬁtting with different degrees used in DFA. The second one is “the jump process” presented in 5 which has some similarities with the method used in the ACD algorithm to dampen the ﬂuctuations. Thus I have chosen for comparison one representative from each of the two major types of trend removal methods parametric and nonparametric extensively 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 artiﬁcial 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 ﬁt 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 brieﬂy 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 t at the moments t n = n -1t, t being the sampling interval. The time takes values within the ﬁnite interval t  0,1. We assume that Z n does not depend on the trend values f n and f t has 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 justiﬁcation of the ACD method can be given. We denote by p z z the probability distribution of Z n and by p x x , n that 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 inﬁnitesi- 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 deﬁne 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 z  z, 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 deﬁnition 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. 3 becomes g = 1 N    n=1 N-1  f n+1 - f n  p z  - f n  +    . 4 The ﬁrst term is the average of the one-step variation of the trend within the neighborhood of . A similar relation is ob- tained if the ﬁnal 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 g and 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 -z and 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 = at, from Eqs. 4 and 2 it follows that g = at, i.e., we ﬁnd 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 Ft estimated 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 argument should be proportional to the average one- step displacement g = FF -1 t . 6 This relation holds only if Ft is invertible, i.e., if g pre- serves the same sign over all its domain of deﬁnition. Our goal is to ﬁnd the numerical quantity corresponding to the theoretical one deﬁned in Eq. 3 using the values of a time series x n  obtained as a realization of the stochastic process X n . The inﬁnitesimal intervals I  must be replaced with ﬁnite 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   deﬁned 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 ﬁnal 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 ﬁnal 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 deﬁnition 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 ˆ t and corresponding to the theoreti- cal one in Eq. 6. Since the initial f 1 and the ﬁnal f N values of the trend are unknown, we use instead the extreme values of the time series x n . The domain of deﬁnition t  0, T of the function F ˆ t is 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 ˆ  + nt, 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 ﬂuctuations described by   causes the quantities gˆ s to have different signs, then the estimated trend F ˆ n cannot be determined and the ﬂuctuations 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 ﬁrst n + K the last N - n + K +1 values. 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 ﬁrst 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 efﬁciency 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- ﬁcially 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 ﬁrst time the last time a 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. 5 which 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 inﬁnitesimal intervals I  with the ﬁnite 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 ﬁnal 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. 1 superposing 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 ﬁt, 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-ﬁt 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 deﬁnition domain: the minimum at t =0 and the maximum at t =1. It is more difﬁcult to remove such a trend than one with the maximum slope in the interior of the deﬁnition domain because numerical algorithms are more difﬁcult 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 AR1autoregressive of order one deﬁned 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 AR1 stochastic 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 1a contains two signals with AR1 noise, 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 7 for 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 1b shows 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–10 x 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 ﬂuc- 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 sufﬁcient diversity of the time series characteristics on which the ACD algorithm is tested. The signals in Fig. 1a correspond 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,1 represented with a thin line. The values of the trend parameter are a =2 thick line and a = 1.1 point markers. a Stationary AR1 noise with  = 0.9 and  Z =1. b Nonstationary 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 AR1 noise 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 justiﬁcation 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 12 on 1000 time series with AR1 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. 12 is 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. 13 is missing from the ﬁgure. Such situations occur when the noise standard deviation is very small see Appen- dix A and 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 12 is 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: 1 establish the values of two parameters the maximum value  of the ratio between the initial signal stan- dard deviation and the ﬁnal residual standard deviation and the maximum value K f of the semilength of the averaging interval; and 2 distribute the signal values into disjoint intervals I s . In this section we show how these tasks can be automatically accomplished. First we analyze the inﬂuence 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 ﬂuctuations 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 = 100 the 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 12 of 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.1 and the minimum one a = 2.0 of the trend parameter for time series with AR1 noise. 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 AR1 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 ﬁnal 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. 12 does 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 ﬁrst the trend is removed from the time series x n  and the estimation 12 is 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. 12 to 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. 14 can 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 8 and 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 4a shows 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 AR1 with  Z =1 and  = 0.9 con- tinuous lines or 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. 14 assures 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 , 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 . In panels a, c, and d the random component of the signals is a stationary AR1 noise with  Z =1 and  = 0.9 continuous lines and a nonstationary fractional Brownian motion with  = 1.4 dashed lines. In panel b the random component is a stationary AR1 noise with  Z =0.1 and  = 0.9 continuous lines, and  =0 dashed lines. a and b The variable parameter is the maximum value  of the ratio between the standard deviation of the initial signal and that of the ﬁnal residual, while the other two ACD parameters have random values. A log-log scale is used. c The 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. d The 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 c and 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 ﬁnal 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 ﬂuctua- tions, not the trend. As seen in Fig. 4a signals 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 14 implies a special procedure to calculate  Z est when  Z has small values. The validity of this procedure is veriﬁed in Fig. 4b for time series with AR1 noise with the minimum standard deviation  Z = 0.1. In all cases the average evaluation index  reaches a stationary value for  =  opt after signiﬁcant 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 ﬁrst parameter given by Eq. 14, we go on to analyze the parameter K f controlling efﬁciency 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 ﬂuctuations are, but at the same time the more distorted the trend is. Figure 4c contains the results for signals with AR1 noise with  Z =1 and  = 0.9 continuous lines and 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 ﬂuctuations 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 ﬂuctuations 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 ﬂuctuations 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. 16 is 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 satisﬁed 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 4d presents 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. 4d the 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. 4d one 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.1 there are values of S for which the accuracy obtained with homogeneous intervals is signiﬁ- cantly better. For all the analyzed signals with AR1 noise results not presented here the 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 veriﬁcation 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 ﬁnal automatic form of the ACD algorithm with nonhomogeneous intervals applied to the artiﬁcial 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.1 or 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 efﬁcient 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 AR1 noise and 0.8 CPU seconds for 1/ f noise. The maximum number of partial estimated trend components was 4, and it was not signiﬁcantly inﬂuenced 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. 15 of 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 ﬁt of degrees from 1 to 10. First we notice that the polynomial ﬁtting 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 ﬁnal 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- ﬁcients to be determined by least-squares ﬁt increases too and they are more strongly affected by the random compo- nent of the signal. Therefore the accuracy of the polynomial ﬁt 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. 6b and 6c one can see that, as expected, the polynomial ﬁtting gives better results if the polynomial has a larger smaller degree when the signal is dominated by trend noise. This contradictory behavior makes very difﬁcult an automatic choice of the optimum value of the polynomial degree for all the possible types of time series. From the ﬁrst three panels of Fig. 6, it is difﬁcult to draw a deﬁnite conclusion regarding the relation between the ACD accuracy and that of the polynomial ﬁt. For  Z =1 continu- ous lines in Fig. 6a, the ACD algorithm has better perfor- mance for short time series N = 100, when the ﬁtted poly- nomial follows more the noise ﬂuctuations 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. 6b  Z = 0.1 the 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 signiﬁcantly better and for degrees higher than 5 signiﬁcantly worse. For the signals dominated by noise in Fig. 6c  Z =10, the ACD is better in almost all cases. Although the accuracy of the ACD algorithm is not signiﬁcantly better than that of the polynomial ﬁtting, 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. 6a, 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 AR1 noise 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 difﬁcult 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. 6d shows 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 identiﬁcation of the trend monotony. The second comparison method is the “jump process” in 5 and it consists in an iterative weighted averaging x n i+1 = x n i + Rx n-1 i -2x n i + x n+1 i  , 17 where ratio R 0  R  0.5 and iteration parameter M i  M are user-speciﬁed constants. At boundaries a sym- metric extension of data is employed. As in the case of the polynomial ﬁtting, the averaging in Eq. 17 has 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 ﬂuctuations 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 ﬁxed. For signals with AR1 noise, excepting a few cases with  Z = 0.1 in Fig. 7b, the ACD method has better accuracy. For signals with 1/ f dashed lines in Fig. 7a, although not so bad as for the polynomial ﬁt, the results of the ACD algorithm continue to be worse than those of the comparison method. As discussed above, comparing Figs. 7b and 7c one notices that, as expected, the jump process has better results if the averaging parameter R is smaller larger when the signal is dominated by trend noise. From Fig. 7d, 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 8a shows 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 ﬁt 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 . a The random component is a stationary AR1 noise with  Z =1 and  = 0.9 continuous lines and a nonstationary fractional Brownian motion with  = 1.4 dashed lines. b The random component is a stationary AR1 noise with  Z =0.1 and  = 0.9 continuous lines and  =0 dashed lines. c The random component is a stationary AR1 noise with  Z =10 and  = 0.9 continuous lines and  =0 dashed lines. d. The random component is as in panel a but 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 9 and 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. 8b the ACD estimated trend is compared with those estimated by the other two methods. As it follows from Fig. 6d, most of the estimated trends by nonlinear polynomial ﬁt 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 ﬁxed 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 , 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 . a The random component is a stationary AR1 noise with  Z =1 and  = 0.9 continuous lines and a nonstationary fractional Brownian motion with  = 1.4 dashed lines. b The random component is a stationary AR1 noise with  Z =0.1 and  = 0.9 continuous lines and  =0 dashed lines. c The random component is a stationary AR1 noise with  Z =10 and  = 0.9 continuous lines and  =0 dashed lines. d. The random component is as in panel a but 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 9 and the trend estimated by means of ACD algorithm continuous line. In panel b the trends estimated by means of the linear ﬁt dashed line and jump process dotted line are 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. 8a is such a case. In order to verify that the estimated trend F ˆ n  is signiﬁcant, we build surrogate series xˆ n = F ˆ n +  ˆ n . Since the serial dependence of the noise has a signiﬁcant inﬂuence on the estimated trend quality, we choose the surrogate series  ˆ n  as realizations of an AR1 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. 10 we 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. 8a, 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 conﬁdence 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 “signiﬁcant,” 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 signiﬁcance 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 = NN -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 signiﬁcance 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  sign or decreasing  sign trend. The three bigger signs correspond to signiﬁcant trends with a con- ﬁdence 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 modiﬁcations. 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% signiﬁcance level applied to AR1 noise without trend, in terms of the correlation param- eter , for different lengths of the time series. The ﬁgure 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. 6d for nonlinear polynomial ﬁt and Fig. 7d for 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 = pt + a sin 2t, 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 t has an inﬂection 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 deﬁned 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 1 is obtained for white noise  =0 and longer series N =50. For correlated noise the results are worse because the probability of type I error is much larger than the signiﬁcance 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 2t, t  0,1, without superposed noise. One notices that the estimated trend is signiﬁcantly different from the “real” monotone component of the trend, which is a straight line. This example shows that for a nonmonotone trend one cannot deﬁne unequivocally a monotone component and the ACD monotone estimated trend is only one of the inﬁnitely 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 2t with variable a, such that the trend is continu- ously modiﬁed from a linear one a =0 to the nonmonotone one represented in Fig. 12 a =2a 0 . Figure 13a shows 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% signiﬁcance level applied to AR1 noise 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% signiﬁ- cance level dashed lines for time series with different lengths N obtained by superposing AR1 noise with  Z =1 on the monotone trend f t = pt + a 0 sin 2t, 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 efﬁcient 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 ﬁt increases for smaller degrees, we compare in Fig. 13b the ACD accuracy only with the linear ﬁt. For small values of parameter a the trend is almost linear, such that the accuracy of the linear ﬁt is better in all the cases. For greater values of the parameter a and for small noise, the ACD accuracy becomes better N =50 or at least comparable N =10. Only for a signal dominated by white noise  =0 and  Z =1 the results of the linear ﬁt 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. 13c. 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 artiﬁcial signals with AR1 noise indicates that the quality of the estimated trend is slightly better than that of the usual nonautomatic methods. This conclusion is conﬁrmed 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 ﬁt 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 difﬁcult 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 artiﬁcial 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 line for a time series continuous line with 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. a Dependence of the average evaluation index of the ACD algorithm on the parameter a in the trend formula f t = t + a sin 2t, 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 AR1 noise with standard deviation  Z and different serial dependence  = 0.9 continuous lines and  =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. b Comparison of the ACD accuracy with that of the linear ﬁt for the same time series as for panel a. c Comparison 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 ﬁt 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. 7 is computed separately for the positive, respectively, negative variations of the time series values. In this case a unique trend cannot be determined because the ﬁnal 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 ﬂuid is obtained by averag- ing the two-particle distribution function. In 13 it 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 14 and 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 justiﬁcation not a rigorous proof of estimation 12 of the noise standard deviation for a realization x n  of the sto- chastic process in Eq. 1 is 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. A1 decreases 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 ﬁrst 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. A1 the term due to the trend. If m 0 is large enough to render Z m 0 +n and Z n independent, then Eq. A1 becomes  m 0 x n  2  2z 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. A3 as 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 12 on 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. A2 is 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 12 can 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. A1 vanishes and according to Eq. A2   m x n  is increasing. Then there is no m 0 satisfying Eq. 13 and 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 12 is 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 markers and  Z = 0.1 points markers. The horizontal lines represent the nonvanishing value of the noise standard deviation and the standard deviation of an inﬁnite 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 inﬁnite sinusoid   =1/  2. In other words, a sinusoid con- taining more than  2 periods is erroneously identiﬁed 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 ﬁrst of all in the ability of estimation 12 to 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  200 then 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. 12 gives 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 12 will be the subject of a future paper. As a ﬁrst application of Eq. 12 we 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 signiﬁcant prob- ability that no m 0 satisfying Eq. 13 should exist and the noise is interpreted as a trend. For example, for N =4 the minimum value of N for which Eq. 13 is 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 AR1 noise, 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 ﬂuctuations 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. 12 as 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 deﬁned by Eq. 11 can 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 16 applied to a signal dominated by noise  est  1 generates 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 deﬁned 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 ﬂuc- 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

1 G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forcast- ing and Control, 2nd ed. Holden-Day, San Francisco, 1976. 2 D. B. Percival and A. T. Walden, Wavelet Methods for Time Series Analysis Cambridge University Press, Cambridge, En- gland, 2000. 3 C. K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stan- ley, and A. L. Goldberger, Phys. Rev. E 49, 1685 1994. 4 K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E. Stanley, Phys. Rev. E 64, 011114 2001. 5 S. Zhao and G. W. Wei, Comput. Stat. Data Anal. 42, 219 2003. 6 P. J. Brockwell and R. A. Davies, Time Series: Theory and Methods, 2nd ed. Springer-Verlag, New York, 1991. 7 N.-N. Pang, Y.-K. Yu, and T. Halpin-Healy, Phys. Rev. E 52, 3224 1995. 8 E. L. Andreas and G. Trevino, J. Atmos. Ocean. Technol. 14, 554 1997. 9 P. D. Jones and M. E. Mann, Rev. Geophys. 42,1 2004. 10 P. 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/ 11 M. G. Kendall, Rank Correlation Methods, 4th ed. Grifﬁn, London, 1975. 12 R. M. Hirsch and J. R. Slack, Water Resour. Res. 20, 727 1984. 13 C. Vamoş, A. Georgescu, N. Suciu, and I. Turcu, Physica A 227, 81 1996. 14 C. Vamoş, N. Suciu, and A. Georgescu, Phys. Rev. E 55, 6277 1997. 15 C. Vamoş, N. Suciu, and W. Blaj, Physica A 287, 461 2000. CĂLIN VAMOŞ PHYSICAL REVIEW E 75, 036705 2007 036705-16

soon

2007

On the convergence of the Newton-GMBACK method

February 28, 2020

Best uniform approximation of semi-Lipschitz function by extension

May 2, 2023

On a Steffensen type method

February 4, 2018

Automatic algorithm for monotone trend removal

Abstract

Authors

Keywords

Cite this paper as:

References

PDF

About this paper

Journal

Publisher Name

DOI

Print ISSN

Online ISSN

MR

ZBL

Google Scholar citations

References

Paper in html format

Related Posts