cultivars.spectral.periodogram#
Nonparametric spectra: the data’s own frequency-domain story, model-free.
density reads the spectrum off a fitted
system and inherits every specification error the fit carries. This
module estimates it from the data alone, on the same Fourier grid and
with the same \(1/2\pi\) normalization, so the two can be drawn on
one plot and their disagreement read as the model’s misspecification at
each frequency. The starting point is the periodogram of a demeaned
sample,
which is unbiased for \(f(\omega)\) asymptotically but inconsistent – each ordinate is \(f(\omega_j)\chi^2_2/2\) however long the sample – so every estimator here trades resolution for stability in a stated way and reports both sides of the trade. The smoothed periodogram averages \(2m + 1\) adjacent ordinates with a modified Daniell kernel; Welch averages periodograms of half-overlapped Hann-tapered segments on the segment’s coarser grid; Thomson’s multitaper averages periodograms of the same full sample under \(K\) orthogonal Slepian tapers, which is the least leakage-prone of the three at a given variance and the default. Every ordinate of every estimator is approximately \(f(\omega)\chi^2_\nu/\nu\), and the result carries \(\nu\) so a confidence band is one method call.
Two commitments shape the surface. First, the result is a distinct type
from SpectralDensityResult on
purpose: a number from a periodogram and a number from a fitted model
answer different questions, and the package never lets one be mistaken
for the other – the two overlay, they do not interconvert. Second,
resolution and degrees of freedom are stated on the result rather than
left to the reader’s memory of the method: bandwidth says how close
two peaks can be before they merge and degrees_of_freedom says how
wide the band around each ordinate is, so the smoothing choice can be
audited from the printed summary.
Layout. MultitaperSpectrum, DaniellSpectrum and
WelchSpectrum take a panel at construction, validate the
smoothing design, and estimate through compute; all three return
SpectrumEstimate, which exposes per-series spectra, squared
coherence, periods and the chi-squared band. The numerics live in
_core: _cross_periodogram forms the \((n, k, k)\) tapered
cross-periodogram, _fourier_frequencies the grid,
_multitaper_density, _daniell_density and _welch_density the
three averages with their degrees of freedom, _ols_detrend removes
the deterministic terms, _validate_spectrum_panel coerces the input
and resolves labels with the _MIN_SPECTRUM_OBS floor of 32, and
_SPECTRAL_METHOD_TITLES names the methods in summaries.
_SummaryMixin from _internals gives the result its printed form.
The model-implied density these estimates are meant to be overlaid on
is SpectralDensity; the band-pass
filters that isolate whatever band the spectrum shows are in
band_pass.
References
Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096.
Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73.
Bloomfield, P. (2000). Fourier Analysis of Time Series: An Introduction, 2nd ed. Wiley.
Percival, D. B., & Walden, A. T. (1993). Spectral Analysis for Physical Applications. Cambridge University Press.
Example
An AR(2) with a spectral peak near 0.55 radians, estimated from the data and compared with two fitted models on the same grid. The VAR(1) is monotone by construction and its density falls inside the 95% band at less than a third of the ordinates; the VAR(2) finds the peak and sits inside the band almost everywhere:
>>> import numpy as np
>>> from cultivars.multivariate.reduced_form.vector_autoregression import VAR
>>> from cultivars.spectral.density import SpectralDensity
>>> rng = np.random.default_rng(0)
>>> y = np.zeros(1024)
>>> for t in range(2, 1024):
... y[t] = 1.4 * y[t - 1] - 0.7 * y[t - 2] + rng.standard_normal()
>>> est = MultitaperSpectrum(y).compute()
>>> lower, upper = est.confidence_interval("y")
>>> inside = {}
>>> for order in (1, 2):
... model = SpectralDensity(VAR(y, order=order).fit(), n_frequencies=513).compute()
... power = model.spectrum("y1")
... inside[order] = round(float(np.mean((lower <= power) & (power <= upper))), 2)
>>> inside
{1: 0.3, 2: 0.97}
Classes#
The periodogram smoothed with a modified Daniell kernel. |
|
Thomson's multitaper spectral estimate. |
|
A nonparametric spectral density matrix, with its resolution and degrees of freedom. |
|
Welch's averaged periodogram over half-overlapped Hann segments. |