Brownian Motion in Space - Algorithms

Algorithms used

The algorithms used in the process of generating the tables of distribution functions of the suprema can be separated into two parts. Firstly there is the process of simulating the brownian motion itself, and secondly there is the process by which each of the simulated occurrences is checked to find the supremum.

2-dimensional uncorrelated

In the two-dimensional uncorrelated case, we have:

The distribution function of supremum over (s,t) in [0,1]² of |W(s,t)| is the limit as n tends to infinity of, and is approximated by,
P{ | ξ_n(s,t) - nst | / √n < x for all 0 < s,t < 1}
where ξ_n(s,t) is a poisson process on [0,1]² with expected value E[ ξ_n(s,t) ] = nst.

Simulating the Brownian Motion

The simulation of the normalised two-dimensional poisson process, which approximates the brownian motion, is achieved by first selecting the number of points, then selecting the locations within the unit square for each of these points.

This is achieved by sampling once from a poisson distribution of mean n to get the value N, then twice sampling N values from a uniform distribution on [0, 1] to get the x- and y-coordinates of each point. The R language code to achieve this is:

  n <- 5000
  N <- rpois(1, n)
  X1 <- runif(N)
  X2 <- runif(N)

The three values of n, X1 and X2, where X1 and X2 are vectors of length N, completely define the trajectory of the simulated poisson process in two dimensions on the unit square.

Now given the values of n, X1 and X2, ξ_n(s,t) can be calculated as _i=1Σ^N Ι_{{X1i≤s, X2i≤t}}
where Ι is the indicator function, and N is the length of the X1 and X2 vectors.

Calculating the Supremum

Calculating the supremum of | ξ_n(s,t) - nst | / √n now involves checking the value of the expression at each of the values (X1_i, X2_j), i=1,..N, j=1,..N. Since the value of ξ_n(s,t) is a step function, and nst defines a monotonic surface, it is not necessary to consider any values of (s,t) other than those defined by X1 and X2, except along the upper boundaries of the unit square (s=1 or t=1). However note that there may be two possible values of ξ_n(s,t) at each of these points, corresponding to the 'bottom' and 'top' of a 'step', although in practise a step will occur at approximately half of these points. Thus overall there are approximately (N+1)² values that may be calculated. This is a large number (approximately 2.5*10⁷ when n is 5000), particularly remembering that this process has to be repeated 20000 times to produce the tables provided.

The nature of the step function ξ_n(s,t) is demonstrated with an example plot which has been generated for n=10 and N=10. A plot of the corresponding function ξ_n(s,t) - nst is also shown in example plots [PDF, 108KB]

When using a high-level language such as R, the advantage of ease of programming has to be traded off against the disadvantage of the inefficiencies of execution of an interpreted language. With R though there is a distinct value to be gained in performing computations in a vectorised (or even matrix) form, rather than one value at a time.

For the calculation of ξ_n(s,t) the algorithm used is to cycle through the values of X1 (plus the boundary value of s=1) and calculate the expression for all possible values of X2 (plus the boundary value of t=1). These values are then candidates to provide the maximum value of the function of which we are searching for the supremum. These calculations are repeated for a set of points slightly closer to the origin, to provide candidates for the minimum of the target function. For each point in these two sets the value nst is then subtracted from ξ_n(s,t) and the greatest result of these in absolute value is then the supremum for this trajectory.

The R language code to achieve this is:

function(n, X1, X2) {
  N <- length(X1)
  if (N <= n) {
    maxval <- 0
    maxs <- 0
    maxt <- 0
    minval <- N - n
    mins <- 1
    mint <- 1
  } else {	# N > n
    maxval <- N - n
    maxs <- 1
    maxt <- 1
    minval <- 0
    mins <- 0
    mint <- 0
  }
  oX2 <- order(X2)
  X1oX2 <- X1[oX2]
# Add dummy points to force calculations at boundaries
  U1 <- c(X1, 1)
  U2 <- c(X2, 1)
  oU2 <- order(U2)
  rU2 <- order(oU2)	# same as rank(), but faster
  nU2 <- n*U2
  for (i in 1:(N + 1)) {
    s <- U1[i]
    nsU2 <- s*nU2
    val <- cumsum(c(X1oX2 <= s, 0))[rU2] - nsU2
    mval <- max(val)
    if (mval > maxval) {
      maxval <- mval
      maxs <- s
      maxt <- U2[which.max(val)]
    }
    val <- cumsum(c(0, X1oX2 < s))[rU2] - nsU2
    mval <- min(val)
    if (mval < minval) {
      minval <- mval
      mins <- s
      mint <- U2[which.min(val)]
    }
  }
    return(c(c(minval, maxval) / sqrt(n), mins, mint, maxs, maxt))
}

2-dimensional correlated

In the two-dimensional correlated case, we have:

The distribution function of supremum over (s,t) in [0,1]² of |W_H(s,t)| is the limit as n tends to infinity of, and is approximated by,
P{ | ξ_nH(s,t) - nH(s,t) | / √n < x for all 0< s,t < 1}
where ξ_nH(s,t) is a poisson process on [0,1]² with expected value E[ ξ_nH(s,t) ] = nH(s,t).
Here H(s,t) is the following copula function:
H(s,t) = F(F^-1(s), F^-1(t); r)
where F^-1(x) is the quantile function (inverse) of the standard normal distribution function and F(x, y; r) is the bi-variate normal distribution function with mean 0, variances 1 and correlation coefficient r.

Calculating the Supremum

Calculating the supremum of | ξ_nH(s,t) - nH(s,t) | / √n now involves checking the value of the expression at each of the values (X1_i, X2_j), i=1,..N, j=1,..N. Since the value of ξ_nH(s,t) is a step function, and nH(s,t) defines a monotonic surface, it is not necessary to consider any values of (s,t) other than those defined by X1 and X2, except along the upper boundaries of the unit square (F^-1(s)=1 or F^-1(t)=1). However note that there may be two possible values of ξ_n(s,t) at each of these points, corresponding to the 'bottom' and 'top' of a 'step', although in practise a step will occur at approximately half of these points.

The nature of the step function ξ_nH(s,t) is demonstrated with an example plot which has been generated for n=10 and N=10 and r=0.5. A plot of the corresponding function ξ_nH(s,t) - nH(s,t) is also shown in example plots [PDF, 109KB]

For the calculation of ξ_nH(s,t) the algorithm used is as for the uncorrelated case, since this is a similar step function. However, unlike the uncorrelated case, the value nH(s,t) to be subtracted from the step function values is non-trivial to compute.

H(s,t) is the bivariate standard normal distribution function with correlation coefficient r. Analytically this is a double integral which is computationally very intensive. Given that for each of the 20000 simulated trajectories, this double integral needs to be calculated at each of the (N+1)² points (s,t), an alternate algortihm is needed. What was chosen was to use linear interpolation on a non-regular grid of values, where the grid is chosen in such a way that grid values are closer together in places where the function shows a greater change in gradient (i.e. where the function is more curved). Thus to calculate a value for H(s,t), interpolation is performed on the rectangle surrounding the point (s,t). One important aspect of this interpolation is that it is a vectorised function, so that if either (or both) of s or t is a vector, a corresponding vector (or matrix) of interpolated values is returned.

Then for each point in the two sets of values of ξ_nH(s,t) the value nH(s,t) is subtracted and the greatest of these in absolute value is then the supremum for this trajectory.

Finally, for the chosen supremum, a refined value of the trajectory is recalculated using a numerical double integral. Thus the interpolated values for H(s,t) are used to identify the supremum, but then the value of this supremum is refined using a more accurate but time-consuming double integral.