# Simulation of Wiener process

# Number of Schauder functions which are used for approximation 

N = 10

# For N=10 the approximation of the Wiener process is not very good
# If you choose N=256 for example, you get an approximation such that
# the realiezd paths look very similar to the one of a Wiener process.

#N = 256

# For simulation of the Wiener process on [0,2] we have to simulate two independent Wiener processes on
# [0,1], see page 63 of the Lecture notes. 

T1 <- seq(0,1,0.001)
T2 <- seq(1,2,0.001)

Y1 <- rnorm(N)
Y2 <- rnorm(N)

X1 <- rep(0,length(T1))
X2 <- rep(0,length(T1))

for(i in seq(1, N, 1)){
  X1 <- X1 + Y1[i] * SchauderVec(i,T1)
  X2 <- X2 + Y2[i] * SchauderVec(i,T1)
}

X2 <- X2 + X1[length(T1)]

plot(c(T1,T2),c(X1,X2),type='l')

# The n-th Schauder function is evaluated for a single point t 

Schauder <- function(n,t){
  if(n == 1){
    return(t)
  }
  
  m <- floor(log(n) / log(2))
  x <- 2^m
  k <- n - x
  
 
  
  if(abs(k) < 0.000001){
    if(t < (x-2)/x|| t > 1){
      return(0)
    }
    x <- 2^(m-1)
    u <- (2*x-1)/(2*x)
    y0 <- sqrt(x)
    y1 <- y0 * (t - (x - 1)/x)
    y2 <- 1/(y0 * 2) - y0 * (t - (x*2 - 1)/(x * 2))
    
    return(y1 *(t < u) + y2 * (t >= u))
  }
  
  if(t < (k-1)/x || t > k/x){
    return(0)
  }
  
  x <- 2^m
  u <- (2*k-1)/(2*x)
  y0 <- sqrt(x)
  y1 <- y0 * (t - (k - 1)/x)
  y2 <- 1/(y0 * 2) - y0 * (t - (k*2 - 1)/(x * 2))
  
  return(y1 * (t < u) + y2 * (t >= u))
}


# The n-th Schauder function is evaluated for all entries of the vector T

SchauderVec <- function(n,T){
  S <- T
  for(i in seq(1,length(T))){
    S[i] = Schauder(n,T[i])
  }
  return(S)
}