MultitaperSpectrum#

class cultivars.spectral.periodogram.MultitaperSpectrum(data, *, bandwidth=4.0, n_tapers=None, detrend='c', names=None)[source]#

Thomson’s multitaper spectral estimate.

\(K\) discrete prolate spheroidal (Slepian) sequences of time-bandwidth product \(NW\) are the orthogonal tapers that concentrate the most energy inside \(\pm W\) cycles per observation; each gives a nearly leakage-free periodogram \(\hat f_k(\omega)\), and the estimate is their concentration-weighted average

\[\hat f(\omega) = \frac{\sum_{k=1}^{K} \lambda_k\, \hat f_k(\omega)} {\sum_{k=1}^{K} \lambda_k},\]

with \(\lambda_k\) the tapers’ concentration ratios, close to one for \(K \le 2NW - 1\) and falling off after. The average has \(2K\) equivalent degrees of freedom and resolution \(2\pi NW / N\) radians – two peaks closer than that merge – so the trade is explicit: \(NW\) buys smoothness at the cost of resolution, \(K\) buys degrees of freedom up to the number of well-concentrated tapers. The usual \(NW = 4, K = 7\) is the default; raise both for a smoother estimate, lower both to separate close peaks. Among the three estimators here it is the least biased at a given variance, and the default choice.

Variables:

  • _panel – The (nobs, k) panel after detrending.

  • _names – The series labels.

  • _detrend – The validated detrending choice.

  • _bandwidth – The time-bandwidth product \(NW\).

  • _tapers – The taper count \(K\).

Parameters:

  • data (ArrayLike) – A (nobs, k) panel or one-dimensional series, at least 32 observations.

  • bandwidth (float) – Time-bandwidth product \(NW\), at least 1 and below nobs / 2; need not be an integer.

  • n_tapers (int | None) – Tapers \(K\); None takes \(2NW - 1\). At most \(2NW\).

  • detrend (str) – "n", "c" (demean) or "ct" (linear detrend), removed by least squares before tapering.

  • names (tuple[str, ...] | None) – Series labels; None reads them from a frame or uses y1 .. yk (y for a single series).

Raises:

  • SpecificationError – If the bandwidth is below 1, the taper count is not an integer of at least 1 or exceeds \(2NW\), the detrend choice is unknown, the labels do not match the columns, or the panel is too short.

  • DimensionError – If the panel is malformed.

  • NumericalError – If the panel is not finite.

Note

The \(2K\) degrees of freedom are exact for an unweighted average of independent periodograms and slightly optimistic for the weighted one when the last tapers are poorly concentrated: at \(NW = 4\) the eighth ratio is 0.70, and the effective count \(2 / \sum w_k^2\) is 15.8 rather than 16. Keep \(K \le 2NW - 1\) unless the extra taper is worth the leakage it admits. Detrending matters more here than for a model-implied density: an unremoved linear trend puts power at the origin that leaks through the tapers’ sidelobes, and "ct" reduces the zero-frequency ordinate of a trending noise series by three orders of magnitude.

See also

References

Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096.

Percival, D. B., & Walden, A. T. (1993). Spectral Analysis for Physical Applications, chapter 7. Cambridge University Press.

Example

>>> rng = np.random.default_rng(0)
>>> y = np.zeros(400)
>>> for t in range(1, 400):
...     y[t] = 0.7 * y[t - 1] + rng.standard_normal()
>>> est = MultitaperSpectrum(y).compute()
>>> est.density.shape, bool(est.spectrum("y")[1] > est.spectrum("y")[-1])
((201, 1, 1), True)

Degrees of freedom are \(2K\) and the resolution is \(2\pi NW / N\); halving \(NW\) halves both:

>>> wide = MultitaperSpectrum(y, bandwidth=4).compute()
>>> narrow = MultitaperSpectrum(y, bandwidth=2).compute()
>>> wide.degrees_of_freedom, narrow.degrees_of_freedom
(14.0, 6.0)
>>> round(wide.bandwidth / narrow.bandwidth, 1), round(400 * wide.bandwidth / 2 / np.pi, 1)
(2.0, 4.0)

Methods#

compute()

Estimate the density matrix.