Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

This post is the continuation of another post related to a generic method for outlier detection in time series. Basically, at this point I'm interested in a robust way to discover the periodicity/seasonality of a generic time series affected by a lot of noise. From a developer point of view, I would like a simple interface such as:

unsigned int discover_period(vector v);

Where v is the array containing the samples, and the return value is the period of the signal. The main point is that, again, I can't make any assumption regarding the analyzed signal. I already tried an approach based on the signal autocorrelation (detecting the peaks of a correlogram), but it's not robust as I would like.

share|improve this question

5 Answers 5

If you really have no idea what the periodicity is, probably the best approach is to find the frequency corresponding to the maximum of the spectral density. However, the spectrum at low frequencies will be affected by trend, so you need to detrend the series first. The following R function should do the job for most series. It is far from perfect, but I've tested it on a few dozen examples and it seems to work ok. It will return 1 for data that have no strong periodicity, and the length of period otherwise.

Update: Version 2 of function. This is much faster and seems to be more robust.

find.freq <- function(x)
    n <- length(x)
    spec <- spec.ar(c(x),plot=FALSE)
    if(max(spec$spec)>10) # Arbitrary threshold chosen by trial and error.
        period <- round(1/spec$freq[which.max(spec$spec)])
        if(period==Inf) # Find next local maximum
            j <- which(diff(spec$spec)>0)
                nextmax <- j[1] + which.max(spec$spec[j[1]:500])
                period <- round(1/spec$freq[nextmax])
                period <- 1
        period <- 1
share|improve this answer
Thank you. Again, I will try this approach as soon as possible and will write here the final results. –  gianluca Aug 4 '10 at 14:49
Your idea is quite good, but in my case, it fails to detect the periodicity of a really simple (and not so noisy) time series like dl.dropbox.com/u/540394/chart.png. With my "empirical" approach (based on the autocorrelation), the simple algorithm I wrote returns an exact period of 1008 (having a sample every 10 minute, this means 1008/24/6 = 7, so a weekly periodicity). My main problems are: 1) It's too slow to converge (it requires a lot of historical data) and I need a reactive, online approach; 2) It's inefficient as hell from a memory usage point of view; 3) It's not robust at all; –  gianluca Aug 4 '10 at 18:14
Try v2......... –  Rob Hyndman Aug 5 '10 at 4:04
Thank you. Unfortunately, this still doesn't work as I would expect. For the same time series of the previous comment it returns 166, which is only partially right (from my point of view, the evident weekly period is more interesting). And using a very noisy time series, like this one dl.dropbox.com/u/540394/chart2.png (a TCP receiver window analysis), the function returns 10, while I would expect 1 (I can't see any obvious periodicity). BTW I know that it will be really difficult to find what I'm looking for, since I'm dealing with too different signals. –  gianluca Aug 5 '10 at 16:17
166 is not a bad estimate of 168. If you know the data are observed hourly with a weekly pattern, then why estimate the frequency at all? –  Rob Hyndman Aug 6 '10 at 3:01

If you expect the process to be stationary -- the periodicity/seasonality will not change over time -- then something like a Chi-square periodogram (see e.g. Sokolove and Bushell, 1978) might be a good choice. It's commonly used in analysis of circadian data which can have extremely large amounts of noise in it, but is expected to have very stable periodicities.

This approach makes no assumption about the shape of the waveform (other than that it is consistent from cycle to cycle), but does require that any noise be of constant mean and uncorrelated to the signal.

chisq.pd <- function(x, min.period, max.period, alpha) {
N <- length(x)
variances = NULL
periods = seq(min.period, max.period)
rowlist = NULL
for(lc in periods){
    ncol = lc
    nrow = floor(N/ncol)
    rowlist = c(rowlist, nrow)
    x.trunc = x[1:(ncol*nrow)]
    x.reshape = t(array(x.trunc, c(ncol, nrow)))
    variances = c(variances, var(colMeans(x.reshape)))
Qp = (rowlist * periods * variances) / var(x)
df = periods - 1
pvals = 1-pchisq(Qp, df)
pass.periods = periods[pvals<alpha]
pass.pvals = pvals[pvals<alpha]
#return(cbind(pass.periods, pass.pvals))
return(cbind(periods[pvals==min(pvals)], pvals[pvals==min(pvals)]))

x = cos( (2*pi/37) * (1:1000))+rnorm(1000)
chisq.pd(x, 2, 72, .05)

The last two lines are just an example, showing that it can identify the period of a pure trigonometric function, even with lots of additive noise.

As written, the last argument (alpha) in the call is superfluous, the function simply returns the 'best' period it can find; uncomment the first 'return' statement and comment out the second to have it return a list of all periods significant at the level alpha.

This function doesn't do any sort of sanity checking to make sure that you've put in identifiable periods, nor does it (can it) work with fractional periods, nor is there any sort of multiple comparison control built in if you decide to look at multiple periods. But other than that it should be reasonably robust.

share|improve this answer

You may want to define what you want more clearly (to yourself, if not here). If what you're looking for is the most statistically significant stationary period contained in your noisy data, there's essentially two routes to take:

1) compute a robust autocorrelation estimate, and take the maximum coefficient
2) compute a robust power spectral density estimate, and take the maximum of the spectrum

The problem with #2 is that for any noisy time series, you will get a large amount of power in low frequencies, making it difficult to distinguish. There are some techniques for resolving this problem (i.e. pre-whiten, then estimate the PSD), but if the true period from your data is long enough, automatic detection will be iffy.

Your best bet is probably to implement a robust autocorrelation routine such as can be found in chapter 8.6, 8.7 in Robust Statistics - Theory and Methods by Maronna, Martin and Yohai. Searching Google for "robust durbin-levinson" will also yield some results.

If you're just looking for a simple answer, I'm not sure that one exists. Period detection in time series can be complicated, and asking for an automated routine that can perform magic may be too much.

share|improve this answer
Thank you for your precious informations, I'll look at that book for sure. –  gianluca Aug 6 '10 at 3:57

You could use the Hilbert Transformation from DSP theory to measure the instantaneous frequency of your data. The site http://ta-lib.org/ has open source code for measuring the dominant cycle period of financial data; the relevant function is called HT_DCPERIOD; you might be able to use this or adapt the code to your purposes.

share|improve this answer
A really helpful link, thank you. –  gianluca Aug 11 '10 at 14:59

In reference to Rob Hyndman's post above http://stats.stackexchange.com/a/1214/70282

The find.freq function works brilliantly. On the daily data set I am using, it correctly worked out the frequency to be 7.

When I tried it on only the week days, it mentioned the frequency is 23, which is remarkably close to 21.42857=29.6*5/7 which is the average number of work days in a month. (Or conversely 23*7/5 is 32.)

Looking back at my daily data, I experimented with a hunch of taking the first period, averaging by that and then finding the next period, etc. See below:

    start=1; #also try start=f;
  if(length(freqs)==1){ return(freqs); }
  for(i in 2:length(freqs)){
find.freq.all(dailyts)  #using daily data

The above gives (7,28) or (7,35) depending on if the seq starts with 1 or f. (See comment above.)

Which would imply that the seasonal periods for msts(...) should be (7,28) or (7,35).

The logic appears sensitive to initial conditions given the sensitivity of the algorithm parameters. The mean of 28 and 35 is 31.5 which is close to the average length of a month.

I suspect I reinvented the wheel, what is the name of this algorithm? Is there a better implementation in R somewhere?

Later, I ran the above code in trying all starts of 1 through 7 and I got 35,35,28,28,28,28,28 for the second period. The average works out to 30 which is the average number of days in a month. Interesting...

Any thoughts or comments?

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.