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, s - t).
3. For v >= u >= s >= t >= 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.

Related posts

• Brownian bridge (interpolating Brownian motion)
• The fractal nature of Brownian motion
• Brown noise and other colors

The post Random Fourier series first appeared on John D. Cook.