# Testing significance of peaks in spectral density

We sometimes use spectral density plot to analyze periodicity in time series. Normally we analyze the plot by visual inspection and then try to draw a conclusion about the periodicity. But has the statisticians developed any test to check whether any spikes in the plot are statistically different from white noise? Has the R-experts developed any package for spectral density analysis and for doing that kind of test? Great if someone could help.

Regards,
P.

-
Pressed by @Wesley, I deleted my quick thoughts about autocorrelation functions and periodogram (may be he indeed is a frequency domain analysis guru, but I personally don't think Bartlett, while working with autocorrelations in time domain), but still think that my second suggestion about bootspecdens may be helpful. –  Dmitrij Celov Jun 22 '11 at 23:09
I'm basing my assumption about people's response to 'what is an autocorrelation?' on literature appearances, in which almost all instances where an autocorrelation is used being of the standard, time-domain computed, Barlett autocorrelation. And, unfortunately, this is bad! :) I do appreciate the suggestion of bootspecdens from Dmitrij; looking forward to checking it out. –  Wesley Burr Jun 23 '11 at 0:21

You should be aware that estimating power spectra using a periodogram is not recommended, and in fact has been bad practice since ~ 1896. It is an inconsistent estimator for anything less than millions of data samples (and even then ...), and generally biased. The exact same thing applies to using standard estimates of autocorrelations (i.e. Bartlett), as they are Fourier transform pairs. Provided you are using a consistent estimator, there are some options available to you.

The best of these is a multiple window (or taper) estimate of the power spectra. In this case, by using the coefficients of each window at a frequency of interest, you can compute a Harmonic F Statistic against a null hypothesis of white noise. This is an excellent tool for detection of line components in noise, and is highly recommended. It is the default choice in the signal-processing community for detection of periodicities in noise under assumption of stationarity.

You can access both the multitaper method of spectrum estimation and the associated F-test via the multitaper package in R (available via CRAN). The documentation that comes with the package should be enough to get you going; the F-test is a simple option in the function call for spec.mtm.

The original reference that defines both of these techniques and gives the algorithms for them is Spectrum Estimation and Harmonic Analysis, D.J. Thomson, Proceedings of the IEEE, vol. 70, pg. 1055-1096, 1982.

Here is an example using the included data set with the multitaper package.

require(multitaper);
data(willamette);
resSpec <- spec.mtm(willamette, k=10, nw=5.0, nFFT = "default",
centreWithSlepians = TRUE, Ftest = TRUE,
jackknife = FALSE, maxAdaptiveIterations = 100,
plot = TRUE, na.action = na.fail)


The parameters you should be aware of are k and nw: these are the number of windows (set to 10 above) and the time-bandwidth product (5.0 above). You can easily leave these at these quasi-default values for most applications. The centreWithSlepians command removes a robust estimate of the mean of the time series using a projection onto Slepian windows -- this is also recommended, as leaving the mean in produces a lot of power at the low frequencies.

I would also recommend plotting the spectrum output from 'spec.mtm' on a log scale, as it cleans things up significantly. If you need more information, just post and I'm happy to provide it.

-
To Burr, Silva and Celov - thanks a lot for your interesting answers and suggestions. I look forward to test these estimators. Best regards –  Pantera Jun 22 '11 at 21:53
(+1) this night I thought carefully about your suggestions, and decided that time domain indeed is the last thing (due to lag truncation and weak properties in small samples) to try searching for cycling behavior. What I am personally concerned about are the assumptions for F statistics and the small sample size properties of the suggested scheme. Well and probably it is good to start a separate question regarding the optimal window selection, because there are indeed many. –  Dmitrij Celov Jun 23 '11 at 7:44
There are indeed many window choices, although the two most common are the Discrete Prolate Spheroidal Sequences (or Slepians) and the sine tapers. If you are looking for maximal concentration of energy in a local bandwidth, the Slepians have been proved to be optimal, and in fact are the output from the integral equation form of the spectral density (see the paper I mentioned for full details). As far as the F statistics go, there certainly are some issues with degrees of freedom, but on the whole they work quite well, with ~ 2k-2 dof available. –  Wesley Burr Jun 23 '11 at 19:56

Ronald Fisher proposed an exact test of the maximum periodogram coordinate in R.A. Fisher, Proc. R. Soc A (1929) 125:54. The test is based on the g-statistic. Specifically, the null hypothesis of Gaussian white noise is rejected if g is significantly large, that is, if one of the values ​​of $f(\omega_k)$ is significantly higher than the average of all signals.

You can get more details about the test in MB Priestley, Spectral Analysis and Time Series, Academic Press, London, 1981, page 406.

In R, the package GeneCycle contains the function fisher.g.test():

library(GeneCycle)
?fisher.g.test


Hope this helps.

-