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

Answer 4

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 5

A different approach could be Empirical Mode Decomposition. The R package is called EMD developed by the inventor of the method:

require(EMD)
ndata <- 3000  
tt2 <- seq(0, 9, length = ndata)  
xt2 <- sin(pi * tt2) + sin(2* pi * tt2) + sin(6 * pi * tt2) + 0.5 * tt2  
try <- emd(xt2, tt2, boundary = "wave")  
### Ploting the IMF's  
par(mfrow = c(try$nimf + 1, 1), mar=c(2,1,2,1))  
rangeimf <- range(try$imf)  
for(i in 1:try$nimf) {  
plot(tt2, try$imf[,i], type="l", xlab="", ylab="", ylim=rangeimf, main=paste(i, "-th IMF", sep="")); abline(h=0)  
}  
plot(tt2, try$residue, xlab="", ylab="", main="residue", type="l", axes=FALSE); box()

The method was branded 'Empirical' for a good reason and there is a risk that the Intrinsic Mode Functions (the individual additive components) get mixed up. On the other hand the method is very intuitive and may be helpful for a quick visual inspection of cyclicity.

Answer 6

In reference to Rob Hyndman's post above https://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:

find.freq.all=function(x){  
  f=find.freq(x);
  freqs=c(f);  
  while(f>1){
    start=1; #also try start=f;
    x=period.apply(x,seq(start,length(x),f),mean); 
    f=find.freq(x);
    freqs=c(freqs,f);
  }
  if(length(freqs)==1){ return(freqs); }
  for(i in 2:length(freqs)){
    freqs[i]=freqs[i]*freqs[i-1];
  }
  freqs[1:(length(freqs)-1)];
}
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?

Answer 7

One can also use Ljung-Box test to figure out which seasonal difference reaches to best stationarity. I was working on a different subject and I used this actually for the same purposes. Try different periods such as 3 to 24 for a monthly data. And test each of them by Ljung-Box and store Chi-Square results. And choose the period with the lowest chi-square value.

Here is a simple code to do that.

minval0 <- 5000 #assign a big number to be sure Chi values are smaller
minindex0 <- 0
periyot <- 0

for (i in 3:24) { #find optimum period by Qtests over original data

        d0D1 <- diff(a, lag=i)

        #store results
        Qtest_d0D1[[i]] <- Box.test(d0D1, lag=20, type = "Ljung-Box")

        #store Chi-Square statistics
        sira0[i] <- Qtest_d0D1[[i]][1]
}
#turn list to a data frame, then matrix
datam0 <- data.frame(matrix(unlist(sira0), nrow=length(Qtest_d0D1)-2, byrow = T))
datamtrx0 <- as.matrix(datam0[])
#get min value's index
minindex0 <- which(datamtrx0 == min(datamtrx0), arr.ind = F)
periyot <- minindex0 + 2