# 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]^{2}
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]^{2} 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)

^{2}values that may be calculated. This is a large number (approximately 2.5*10

^{7}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]^{2}
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]^{2} 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)^{2} 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.