Brownian motion via Hilbert basis

Brownian motion is a central object of probability theory. One perspective to look at its construction uses Hilbert spaces.
Let (f_i)_{i\in I} be a Hilbert basis H:=L^2([0,1]) equipped with standard scalar product (\cdot,\cdot). Hence, every f\in H can be written in a unique way as

    \[f  =  \sum_{i\in I} \alpha_{i} f_i,\]

where \alpha_{i}=(f,f_i) and \sum_{i\in I} \alpha_{i}^2=|f|_2^2<\infty. Let (\xi_i)_{i\in I} be an i.i.d. sequeence of standard gaussian random variables and consider

    \[ I: \left{ \begin{array}{ccc}  H & \longrightarrow & L^2(\Omega)\cr f=\sum_{i\in I} \alpha_{i} f_i & \longmapsto & \sum_{i\in I} \alpha_{i} \xi_i. \end{array} \right.\]

Using the independence of \xi_i we obtain

    \[ \mathbb{E} [I(f)I(g)]  =  (f,g).\]

Brownian motion is then given by (a continuous modification of) B_t:= {I}(\mathbf{1}_{[0,t]}), t\ge 0. 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

    \[f_n(x) = \sqrt{2} \cos (\pi n x), \quad n\in \mathbb{N},\]


together with f_0=1.

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.