Period detection of a generic time series

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.

• Have you tried xts::periodicity? – Fabrício Oct 1 '16 at 23:40

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)
}

• 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
• 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
• An improved version is in the forecast package as findfrequency – Rob Hyndman Feb 18 '17 at 0:00

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.

• Looks interesting but I don't understand the output, it doesn't tell me where the period starts, and most pvalues of 1. – Herman Toothrot Nov 20 '17 at 18:30

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.

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

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.

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

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