Random Fourier series

A theorem by Paley and Wiener says that a Fourier series with random coefficients produces Brownian motion on [0, 2π]. Specifically,

W(t) = \frac{Z_0}{\sqrt{2\pi}}\,t + \frac{2}{\sqrt{\pi}} \sum_{n=1}^\infty \frac{Z_n}{n} \sin\left(\frac{nt}{2}\right)

produces Brownian motion on [0, 2π]. Here the Zs are standard normal (Gaussian) random variables.

Here is a plot of 10 instances of this process.

Plot of 10 random Fourier series

Here’s the Python code that produced the plots.

    import matplotlib.pyplot as plt
    import numpy as np
    def W(t, N):
        z = np.random.normal(0, 1, N)
        s = z[0]*t/(2*np.pi)**0.5
        s += sum(np.sin(0.5*n*t)*z[n]/n for n in range(1, N))*2/np.pi**0.5
        return s
    N = 1000
    t = np.linspace(0, 2*np.pi, 100)
    for _ in range(10):
        plt.plot(t, W(t, N))

Note that you must call the function W(t, N) with a vector t of all the time points you want to see. Each call creates a random Fourier series and samples it at the points given in t. If you were to call W with one time point in each call, each call would be sampling a different Fourier series.

The plots look like Brownian motion. Let’s demonstrate that the axioms for Brownian motion are satisfied. In this post I give three axioms as

  1. W(0) = 0.
  2. For s > t ≥ 0. W(s) – W(t) has distribution N(0, st).
  3. For vust ≥ 0. W(s) – W(t) and W(v) – W(u) are independent.

The first axiom is obviously satisfied.

To demonstrate the second axiom, let t = 0.3 and s = 0.5, just to pick two arbitrary points. We expect that if we sample W(s) – W(t) a large number of times, the mean of the samples will be near 0 and the sample variance will be around 0.2. Also, we expect the samples should have a normal distribution, and so a quantile-quantile plot would be nearly a straight 45° line.

To demonstrate the third axiom, let u = 1 and v = √7, two arbitrary points in [0, 2π] larger than the first two points we picked. The correlation between our samples from W(s) – W(t) and our samples from W(v) – W(u) should be small.

First we generate our samples.

    N = 1000
    diff1 = np.zeros(N)
    diff2 = np.zeros(N)
    x = [0.3, 0.5, 1, 7**0.5]
    for n in range(N):
        w = W(x)
        diff1[n] = w[1] - w[0]
        diff2[n] = w[3] - w[2]

Now we test axiom 2.

    from statsmodels.api import qqplot
    from scipy.stats import norm

    print(f"diff1 mean = {diff1.mean()}, var = {diff1.var()}.")
    qqplot(diff1, norm(0, 0.2**0.5), line='45')

When I ran this the mean was 0.0094 and the variance was 0.2017.

Here’s the  Q-Q plot:

quantile-quantile plot to test normal distribution

And we finish by calculating the correlation between diff1 and diff2. This was 0.0226.

Related posts

One thought on “Random Fourier series

  1. Hi Johan,
    A slight change of the code to make it runable:
    x = [0.3, 0.5, 1, 7**0.5] => x = np.array([0.3, 0.5, 1, 7**0.5])
    w = W(x) => w = W(x, N)


Comments are closed.