Period detection of a generic time series

Question

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.

Answer 1

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)
            if(length(j)>0)
            {
                nextmax <- j[1] + which.max(spec$spec[j[1]:500])
                period <- round(1/spec$freq[nextmax])
            }
            else
                period <- 1
        }
    }
    else
        period <- 1
    return(period)
}

Answer 2

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.

Answer 3

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.

Answer 4

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.