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
Ι{X1
i≤s,
X2
i≤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*107
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 |WH(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.