Like what you see? Follow Adrian on Twitter to be notified of new content.

The Wiener-Khinchin theorem provides a clever way of building Gaussian processes for regression. I’ll show you this theorem in action with some Python code and how to use it to model a process.

The Wiener-Khinchin theorem states that an autocovariance function of a weakly stationary process is a function of the power spectral density and vice versa. These two are called Fourier duals and can be written as1: $$ \begin{align} \kappa(h) &= \int S(\omega) e^{2 \pi i \omega h} d \omega \label{1}\tag{1} \\ S(\omega) &= \int \kappa(h) e^{-2 \pi i \omega h} d h \\ \end{align} $$ Where \(\kappa(h)\) is an autocovariance function of a process at lag \(h\) and \(S(\omega)\) is the power spectral density of the same process at frequency \(\omega\). They are Fourier transforms of each other—hence the term Fourier duals.

The idea here is that if you want to model a process, you could estimate its power spectrum (\(S\)), calculate the autocovariance function using equation (\(\ref{1}\)) and construct a Gaussian process.

The problem is that equation (\(\ref{1}\)) is in continuous time and we generally work with discretely sampled time series. Luckily, there is a method known as the periodogram which is a discrete estimate of the power spectral density. SciPy provides this estimate under scipy.signal.periodogram. Now, equation (\(\ref{1}\)) can be treated discretely as:

$$ \kappa(h) = \sum_{\omega} S(\omega) e^{2 \pi i \omega h} d \omega \label{2}\tag{2} $$

Let’s try this out with some Python code.

Estimating the autocovariance function

Generate a time series made up of two frequencies and some noise:

from numpy import sin, pi, cumsum, linspace
from numpy.random import seed, randn


N = 500  # Number of samples
x = linspace(0, 1, N)  # One interval
y = (
    sin(x * 2 * pi * 16)    # Oscillate 16 times per interval
    + sin(x * 2 * pi * 30)  # + 30 times per interval
    + randn(N)              # + some noise

Which looks like this:

We can get the periodogram with:

from scipy.signal import periodogram
frequency, power = periodogram(y)

The array power is a noisy power spectrum with two spikes. One spike at \(\frac{16}{N}\) and the other at \(\frac{30}{N}\). It looks like:

We can create a function that uses periodogram and equation (\(\ref{2}\)) to give us the autocovariance at each \(h\) from \(0\) to \(N\):

import numpy as np
from scipy.signal import periodogram

def spectral_autocov(y):
    N = len(y)
    frequency, power = periodogram(y, nfft=N)
    d = np.mean(np.diff(frequency))
    autocov = np.real([
        np.sum(np.exp(2 * np.pi * 1j*frequency*h) * power * d)
        for h in range(N)

    return autocov

Running the time series y through the spectral_autocov function gives us an array where the \(h^{\text{th}}\) element is the autocovariance at lag \(h\):

autocov = spectral_autocov(y)

Building a Gaussian process

We can use this autocovariance function to construct a covariance matrix and make a prediction with a conditional Gaussian. I’ve talked about the maths of conditional Gaussians in the appendix of a previous article. Knocking this together looks like:

from scipy.linalg import toeplitz

# Forecast this many steps ahead
horizon = 100

# Turn the autocovariance series into a
# covariance matrix plus a noise term
cov = toeplitz(autocov) + np.eye(N) * 20

# Make a forecast using a conditional
# Gaussian process
AA = cov[:(N-horizon), :(N-horizon)]
BA = cov[(N-horizon):, :(N-horizon)]
forecast = BA @ np.linalg.inv(AA) @ y[horizon:]

Which looks like this:


You can calculate the autocovariance function of a process from its power spectral density. Using this relationship, here we explored estimating the power spectral density of a noisy time series, calculating its autocovariance matrix and then making a prediction using a conditional Gaussian.

Like what you see? Follow Adrian on Twitter to be notified of new content.


References & Notes

  1. Andrew Gordon Wilson, Ryan Prescott Adams: "Gaussian Process Kernels for Pattern Discovery and Extrapolation" International Conference on Machine Learning (ICML) 28(3) pp. 1067-1075 (2013) ↩︎


If you see any mistakes or room for improvement, please reach out to me on Twitter @DrAdrian.


Please cite this work as:

Letchford (2023), "Wiener–Khinchin theorem and Gaussian processes", OS Quant.

In BibTeX please use:

    author = {Letchford, Adrian},
    title = {Wiener–Khinchin theorem and Gaussian processes},
    journal = {OS Quant},
    year = {2023},
    note = {},