Power density spectrum formula in R

Question

I am trying to understand the formula that is used when computing the power spectrum for a time series x in R.,

x <- as.ts(c(-0.060280817,0.360106897,-0.065735609,0.030000389,-0.065027960,0.124425013, -0.156030715,-0.231449154,-0.034686819,0.010966803,0.072178488,0.023459309,  0.466334338,-0.121615219,0.340923473,-0.372170475,-0.171057656,-0.487355035,-0.460152330,-0.153500145,-0.095249271,0.079752533,-0.006659656,0.429775255))

I thought that the formula would be

f <- seq(0,0.5,length.out=length(x))

Compute lag-1 autocorrelation (phi) and random noise variance (vare) for x

AR <- ar(x = x, order.max = 1,na.action=na.omit)
phi <- as.numeric(unlist(AR[2]))
varx <- as.numeric(unlist(AR[3]))
vare <- (1-phi*phi)*varx

Compute spectrum

S <- vare/(1+phi*phi-2*phi*cos(2*pi*f))

However the S that I compute like this is totally different from the one computed using the spec.ar function with order 1. Can someone help me clarifying this point?

Thanks

Answer 1

There are many ways to estimate the power spectrum of a time-series. The simplest (and quite frankly, "worst") way to do this is a simple periodogram (scaled absolute value squared of the Fourier transform of the data).

As far as code is concerned, you will usually want to zero-pad so that the frequency domain representation is "nicer" (you're forcing R to compute the value of the estimate at more frequencies). Also, since the spectrum of real data is symmetric about $f = 0$ and has period 1 (assuming sample rate of $\Delta t = 1$), we usually only plot the first $\frac{M}{2}+1$ frequency bins, where $M$ is the total length of the zero-padded series. $M$ should be factorable using small primes (2 for example) in order to take advantage of the efficiency gains of the fast Fourier transform.

N <- length(x)
M <- 2^10 # zero-pad total length
freq <- seq(0, 0.5, by = 1/M)

$\hat{S}^{(P)}(f) = \frac{1}{N}| \sum_{n=0}^{N-1} x_{n} e^{-i2\pi f n} |^{2}$

x.zp <- c(x, rep(0, M-N))
S.pgram <- (1/N)*abs(fft(x.zp)[1:(M/2+1)])^2

plot(freq, S.pgram, type='l', log='y', xlab = "Frequency", ylab = "Spectrum")

Likely, you would like a better estimate, so we use a taper which gives us a "direct spectral estimate" (the periodogram is part of this class, it's taper is the "boxcar taper").

The most optimal taper (with respect to out of band bias / leakage) is a Kaiser taper, or zeroeth order Slepian sequence (aka. direct prolate spheroidal sequence (DPSS)). You can calculate these yourself, but that's silly! There's a package in R called "multitaper" that will do it for you. Letting $h_{n}$ be our tapering sequence:

$\hat{S}^{(D)}(f) = |\sum_{n=0}^{N-1} h_{n} x_{n} e^{-i2\pi f n}|^2$

library(multitaper)
v <- dpss(n = N, nw = 1, k = 1)$v #generates the zeroeth order Slepian
xv <- x*v  # taper times data 
xv.zp <- c(xv, rep(0, M-N))  # now zero-padding

S.kaiser <- abs(fft(xv.zp)[1:(M/2+1)])^2 # the estimate!

plot(freq, S.kaiser, type='l', log='y', xlab="Frequency", ylab = "Spectrum")

Lastly! The Multitaper Spectral estimator is one of the better estimators because you're using a weighted average of $K$ direct spectral estimates using orthogonal tapers (the Slepians). The extra direct estimates gives you variance reduction and the adaptive weighting provides protection against bias.

S.mtm <- spec.mtm(x)

spec.mtm() is provided in the multitaper package and does the zeropadding, etc for you. Unfortunately, you need a bit more data for the multitaper to be useful. The direct estimate is "good enough" in many situations, but if you have more than 100 data points, there's no reason that you wouldn't use a multitaper estimate.

References:

Windowing Functions / Tapers (Wikipedia)

Spectrum Estimation and Harmonic Analysis, David Thomson, 1982, IEEE (paper about Multitaper method)

Spectral Analysis for Physical Applications (book)