# Random Fourier series

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

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.

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))

plt.xlabel("$t$")
plt.ylabel("$W(t)$")
plt.savefig("random_fourier.png")


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')
plt.savefig("qqplot.png")


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

Here’s the  Q-Q plot:

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

## One thought on “Random Fourier series”

1. Jicun Zhong

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)

Regards!
Jicun