Brownian motion is a central object of probability theory. One perspective to look at its construction uses Hilbert spaces.
Let be a Hilbert basis
equipped with standard scalar product
. Hence, every
can be written in a unique way as



Using the independence of we obtain
Brownian motion is then given by (a continuous modification of) ,
. Hence, “all” we need to do is to find a Hilbert basis. The most known is the following.
Haar basis
Here we refer to http://datatreker.com/haar-basis for an introduction.
haar_mother <- function(x){
(x >0 & x <= 0.5) - (x >0.5 & x <= 1)
}
haar <- function(x,j,k){
(2^(j/2) * haar_mother(2^j*x-k))
}
set.seed(42)
j_max <- 10
n_max <- 11 # maximal resolution
delta <- 2^{-n_max}
x <- (0: 2^n_max)/2^n_max
xi_0 <- rnorm(1) # xi corresponding to constant function
xi <- list() # list of random variables xi
for (j in 0:j_max){
xi[[j+1]] <- rep(0, 2^j)
for (k in 0:(2^j-1))
{
xi[[j+1]][k+1] <-rnorm(1)
}
}
# data.frame containing approximation of the Brownian motion
df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())
for (i in x){
y <- c(1,rep(1, i*2^n_max), rep(0, 2^n_max-i*2^n_max))
alpha <- list() # wavelet coefficients for y
bm <- sum(y)/(length(y)+1) * xi_0
for (j in 0:j_max){
alpha[[j+1]] <- rep(0, 2^j)
for (k in 0:(2^j-1))
{
alpha[[j+1]][k+1] <-sum(haar(x,j,k)*y)*delta
bm <- bm + alpha[[j+1]][k+1]*xi[[j+1]][k+1]
}
df_new <- data.frame(x=i,
y=bm,
id=j)
df <- bind_rows(df, df_new)
}
}
#### Animation using gganimate
ggplot(df, aes(x=x, y=y)) +
geom_line()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = "Construction of Brownian motion using Haar basis. Step: {closest_state}") +
ease_aes('sine-in-out')

The result looks the same as with Levy’s construction which is not astonishing. The implementation above is not optimal and is much slower than the previous one.
It might be worth to consider a different Hilbert basis.
Trigonometric Hilbert basis
The trigonometric basis is defined as
together with

set.seed(42)
N_max <- 200 # number of trigonometric basis functions used
n_max <- 11 # maximal resolution
delta <- 2^{-n_max}
x <- (0: 2^n_max)/2^n_max
xi <- rnorm(N_max+1) # list of random variables xi
# data.frame containing approximation of the Brownian motion
df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())
for (i in x){
y <- c(1,rep(1, i*2^n_max), rep(0, 2^n_max-i*2^n_max))
alpha <- numeric() # wavelet coefficients for y
alpha[1] <- mean(y)
bm <- alpha[1] * xi[1]
for (n in 1:N_max){
alpha[n] <-sum(sqrt(2)*cos(x*pi*n)*y)*delta
bm <- bm + alpha[n]*xi[n]
df_new <- data.frame(x=i,
y=bm,
id=n)
df <- bind_rows(df, df_new)}
}
#### Animation using gganimate
p <- ggplot(df, aes(x=x, y=y)) +
geom_line()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = "Construction of Brownian motion using trigonometric basis. Step: {closest_state}") +
ease_aes('sine-in-out')
animate(p, nframes=1000)

The above convergence is very slow. The irregularity properties of the Brownian motion are not yet very visible.