Code guide · Time series

Higher-Order Spectral Analysis in R

A practical introduction to spectra, bispectra, polyspectra, and simple R workflows for exploring nonlinear and non-Gaussian dependence in time series.

1 · Why higher order?

The ordinary spectrum only sees second-order dependence

Second-order spectrum

Describes how variance is distributed across frequencies and is driven by autocovariance structure.

Higher-order spectra

Use third- and higher-order cumulants to detect dependence patterns that can be invisible to ordinary spectral analysis.

Key idea: two time series can have very similar second-order spectra while differing substantially in nonlinear or non-Gaussian structure.
2 · Start with the ordinary spectrum

Frequency-domain analysis begins with the Fourier transform

Time seriesx₁, x₂, …, xₙ
→
FFTcomplex frequency coefficients
→
Periodogrampower by frequency
x <- scale(my_series, center = TRUE, scale = FALSE) n <- length(x) X <- fft(x) freq <- (0:(n - 1)) / n periodogram <- Mod(X)^2 / n

This is the familiar second-order starting point. The periodogram summarizes energy at each frequency, but it does not capture phase coupling among several frequencies.

3 · The bispectrum

Third-order dependence in the frequency domain

The bispectrum is the Fourier-domain analogue of third-order cumulant structure. A common sample building block uses products of three Fourier coefficients whose frequencies satisfy a frequency-sum constraint.

B(ω₁, ω₂) ∝ X(ω₁) X(ω₂) X*(ω₁ + ω₂)
# Illustrative raw bispectral product k1 <- 5 k2 <- 9 k3 <- (k1 + k2) %% n B_raw <- X[k1 + 1] * X[k2 + 1] * Conj(X[k3 + 1]) / n

The magnitude reflects the strength of interaction, while the complex phase carries information about phase relationships among frequencies.

Important: this raw product is an instructional building block, not a smoothed production estimator. Practical bispectral estimation usually requires averaging or smoothing to reduce variance.
4 · Beyond order three

Polyspectra generalize the same idea

Order 2

Spectrum: variance and pairwise temporal dependence.

Order 3

Bispectrum: third-order dependence and frequency coupling.

Higher order

Polyspectra: richer nonlinear and non-Gaussian interaction structure.

The dimension of the frequency-domain object increases with order, so visualization and direct estimation become increasingly challenging. This is one reason summaries such as polyspectral means are useful.

5 · A practical R workflow

From raw series to higher-order summaries

1. PreprocessCenter the series; detrend or difference if scientifically appropriate.
2. TransformCompute the FFT once and reuse the Fourier coefficients.
3. CombineForm products of coefficients satisfying the required frequency constraints.
4. AverageUse smoothing, binning, or weighted integration to stabilize the estimate.
# Example: evaluate a grid of raw bispectral products B <- matrix(NA_complex_, nrow = m, ncol = m) for (i in seq_len(m)) { for (j in seq_len(m)) { k1 <- i - 1 k2 <- j - 1 k3 <- (k1 + k2) %% n B[i, j] <- X[k1 + 1] * X[k2 + 1] * Conj(X[k3 + 1]) / n } } B_mag <- Mod(B)

This code shows the computational pattern. In applied work, choose a nonredundant frequency region and an appropriate smoothing or weighting strategy rather than interpreting every raw cell directly.

6 · Polyspectral means

Summarize a high-dimensional spectrum with a scientifically chosen weight

A polyspectral mean is a weighted average or integral of a higher-order spectrum. Instead of trying to interpret an entire high-dimensional spectral surface, we target a lower-dimensional summary designed around the scientific question.

Mg = ∫ g(ω) fk(ω) dω

Why it helps

Weighted summaries are easier to estimate, compare, test, and use as features for clustering or classification.

Why the weight matters

The weight function g determines which frequency interactions are emphasized, so it should reflect the inferential goal rather than be chosen mechanically.

This perspective underlies my work on Polyspectral Mean Estimation of General Nonlinear Processes, which develops asymptotic inference for polyspectral means of general order and uses bispectral summaries for linearity testing and time-series clustering.

7 · What should you look for?

Interpret higher-order spectra as structure, not decoration

Frequency couplingLocalized peaks can indicate interactions among particular frequency combinations.
Departure from linearityHigher-order spectral features can be used to test or diagnose behavior inconsistent with a simple linear-process model.
Non-Gaussian dependenceThird- and higher-order cumulants can reveal structure absent from covariance-based summaries.
Feature constructionBispectral or polyspectral summaries can become inputs for clustering, classification, or other downstream analyses.
Do not read a raw bispectrum like a heatmap of correlations. Symmetry, redundancy, smoothing, normalization, and the null model all matter for interpretation.
Go deeper

Related work

For theory and formal inference, see the polyspectral-mean paper and the broader time-series work on the Publications page.

Polyspectral mean paper Related publications Back to Tutorials & Resources