SpectrumEstimate#

class cultivars.spectral.periodogram.SpectrumEstimate(*, names, frequencies, density, method, degrees_of_freedom, bandwidth, nobs, detrend)[source]

Bases: _SummaryMixin

A nonparametric spectral density matrix, with its resolution and degrees of freedom.

Frequencies are radians per observation on the Fourier grid \(\omega_j = 2\pi j / N\), \(j = 0, \dots, \lfloor N/2 \rfloor\) (or the segment’s grid under Welch), and the density carries the \(1/2\pi\) normalization, exactly as SpectralDensityResult, so

\[2 \int_0^{\pi} \hat f_{ii}(\omega)\, d\omega \approx \hat\sigma_i^2\]

and the two objects overlay on one set of axes. Every ordinate is approximately \(f(\omega)\,\chi^2_\nu / \nu\) with \(\nu\) the estimator’s equivalent degrees of freedom, which is what confidence_interval() inverts, and the estimate cannot separate two peaks closer than bandwidth radians.

Variables:

  • names – Series labels, indexing the matrix axes.

  • frequencies – The (n,) Fourier grid on [0, pi].

  • density – The (n, k, k) complex Hermitian estimate.

  • method – "daniell", "multitaper" or "welch".

  • degrees_of_freedom – Equivalent chi-squared degrees of freedom of each ordinate.

  • bandwidth – Half-width of the frequency window in radians, the resolution below which two peaks merge.

  • nobs – Observations the estimate used.

  • detrend – The deterministic terms removed first.

Note

This is the data’s periodogram, smoothed, and not a model’s closed form; the two are kept as distinct types so that a disagreement between them can be read as misspecification by frequency rather than lost. The price of consistency is bias: the smoothing window averages the true spectrum over bandwidth radians, so a sharp peak is flattened and the variance identity above holds only approximately (the taper’s leakage and the window’s averaging both cost a little). The summary’s “peak” is the argmax of a noisy estimate over the interior grid; for a monotone spectrum such as an AR(1)’s it lands wherever sampling noise puts it in the low band, so read it with the confidence band, not as a point.

See also

References

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

Priestley, M. B. (1981). Spectral Analysis and Time Series. Academic Press.

Example

>>> rng = np.random.default_rng(0)
>>> y = np.zeros(1024)
>>> for t in range(1, 1024):
...     y[t] = 0.7 * y[t - 1] + rng.standard_normal()
>>> est = MultitaperSpectrum(y).compute()
>>> est.density.shape, est.degrees_of_freedom, round(est.bandwidth, 4)
((513, 1, 1), 14.0, 0.0245)

The estimate integrates to the sample variance, and the true AR(1) spectrum \(\sigma^2 / 2\pi|1 - \phi e^{-i\omega}|^2\) lies inside the 95% band at about 95% of the ordinates:

>>> power = est.spectrum("y")
>>> round(float(2.0 * np.trapezoid(power, est.frequencies) / y.var()), 2)
0.99
>>> theory = 1.0 / (2.0 * np.pi * np.abs(1.0 - 0.7 * np.exp(-1j * est.frequencies)) ** 2)
>>> lower, upper = est.confidence_interval("y")
>>> round(float(np.mean((lower <= theory) & (theory <= upper))), 2)
0.97

The grid is the one a fitted model’s density can be evaluated on, so the two overlay directly:

>>> from cultivars.multivariate.reduced_form.vector_autoregression import VAR
>>> from cultivars.spectral.density import SpectralDensity
>>> model = SpectralDensity(VAR(y, order=1).fit(), n_frequencies=513).compute()
>>> bool(np.allclose(model.frequencies, est.frequencies))
True

Private methods#

_index(name)

Resolve a series label.

_summary_table()

Peak and variance per series, with the estimate's resolution.
On this page