Estimating the spectral density

Today let's understand how to turn a periodogram's jagged spikes into a smooth, reliable estimate of the spectral density hiding underneath them.

Northside Diner is a lunch restaurant that logged how many lunch covers it served every day for ten weeks, seventy days in all. Feed those seventy days into R's spec.pgram(), and it returns a periodogram: one number, called an ordinate, for every candidate frequency, saying how strongly the series moves in step with a wave of that frequency. Two ordinates tower over the rest, a real weekly cycle at frequency 1/7 with ordinate 14,997.4, and a smaller cycle at frequency 2/7, repeating every 3.5 days, with ordinate 3,576.7. Away from those two frequencies, the ordinates are just noise, jumping up and down with no real pattern.

Here is that whole raw periodogram again, all 35 ordinates together.

Two real cycles poking out of a jagged floor. Turning that floor into something you can trust at every frequency, not only the two tallest points, is what today's lesson does.

The spectral density: the curve behind a periodogram

Every raw ordinate spec.pgram() prints, big or small, stands in for something it can never show you directly: the series' true spectral density, a fixed curve saying exactly how much of the series' variance sits at each frequency. Call that curve \(f(\nu)\), a number for every frequency \(\nu\). A raw ordinate at frequency \(\nu\), written \(I(\nu)\), is only an estimate of \(f(\nu)\), and a specific kind of noisy one.

Rebuild Northside Diner's covers and its raw periodogram, and look at the top ordinates again.

RInteractive R
# Rebuild Northside Diner's covers and its raw periodogram set.seed(2024) day <- 1:70 northside_covers <- round(120 + 30 * cos(2 * pi * (day - 5) / 7) + 12 * cos(4 * pi * (day - 5) / 7) + rnorm(70, 0, 10)) pg <- spec.pgram(ts(northside_covers), taper = 0, detrend = FALSE, demean = TRUE, fast = FALSE, plot = FALSE) pg$df top5 <- order(pg$spec, decreasing = TRUE)[1:5] data.frame(frequency = round(pg$freq[top5], 6), ordinate = round(pg$spec[top5], 1)) #> [1] 2 #> frequency ordinate #> 1 0.142857 14997.4 #> 2 0.285714 3576.7 #> 3 0.442857 345.8 #> 4 0.071429 334.1 #> 5 0.271429 307.6

  

pg$df reads 2. Every one of those ordinates, the towering 14,997.4 at 1/7 included, was built from a distribution with exactly 2 degrees of freedom, and that number would still read 2 whether Northside Diner had kept 70 days of records or 7,000. For a series long enough that the usual large-sample approximations hold, twice a raw ordinate divided by the true spectral density at that frequency, \(2I(\nu)/f(\nu)\), behaves like a draw from a chi-squared distribution with 2 degrees of freedom.

That distribution has a mean of 2 and a standard deviation of 2 as well, so its spread relative to its own average, called its coefficient of variation, sits at 100%, permanently. A chi-squared distribution with exactly 2 degrees of freedom even has a simpler name, an exponential distribution, the same shape you'd get timing how long you wait for a rare event, whose standard deviation always equals its mean no matter how you set its rate.

That's the real reason a single raw ordinate never settles down: it carries a fixed 100% relative spread by construction, not because Northside Diner's 70 days happen to be too few. What can shrink that spread is averaging several ordinates together instead of reading just one.