Tutorial for BrownianMotion package

Note the pipe operator in this vignette (|>) assumes you have R 4.1.0+, allowing natural chaining of commands. The package does not require this version of R, only this vignette. We also set the seed to ensure reproducibility of this vignette.

library(BrownianMotion)
set.seed(212)

This tutorial will begin by considering one dimensional standard Brownian motion. In Multiple Dimensions we consider multi-dimensional Brownian motion with user-specified covariance matrix. Unless otherwise stated, all syntax introduced in the one-dimensional case, and the output types, matches that of multi-dimensional Brownian motion.

One Dimension

The first task is to create an object which will contain the Brownian motion path and layer information. This is accomplished using the create.bm() function. It is possible to call with all default options, but here for demonstration purposes we will force the path to start at \(W_{-1} = 100\) (the default is \(W_0 = 0\)).

bm <- create.bm(-1, 100)

Of course, we may be interested in subsequently “forcing” a path constrained to hit particular point(s) at particular time(s) (i.e., bridge simulation). For instance, we may wish the path to “hit” \(W_0 = 101\) and \(W_1 = 99\). This can be achieved using the add.points which has two required arguments, one for specifying a vector of forced times, and the second a vector of forced points:

bm |> add.points(c(0, 1), c(101, 99))

At this point we can see a plot of our nascent path, by simply using the plot() function:

plot(bm)

Here, black dots represent path observations and these are connected by grey lines. Note grey lines do not represent the actual path, although see later section on Publication Graphics for details on plotting an illustrative sample path (TODO).

Alternatively, if at the outset we wish to simulate a Brownian bridge, we can supply a vector of times and position to the create.bm function. For instance, replicating the above:

bm <- create.bm(c(-1, 0, 1), c(100, 101, 99))
plot(bm)

Standard simulation

Unconditional simulation

The simplest operation is unconditional simulation of the path at future times (beyond the end time of the existing path). To do this we use the function sim(), and supply a vector of times. Any of the times supplied which are in excess of the current end point will be simulated in order unconditionally (and all other times will be dealt with appropriately automatically, as discussed later). For instance, noting for our object bm the current end point is \(t = 1\), suppose we wish to simulate the path at the times \(t = (2,3,4)\), we have:

bm |> sim(2:4)

Note that the bm object is updated in place: there is no need to assign the result back to bm (ie you do not need to do bm <- bm |> sim(2:4), although there will be no negative side effects if you do so).

We can again visualise the path by using the plot() function:

plot(bm)

Brownian bridge simulation

The simplest step beyond unconditional simulation would be to simulate times which are intermediate to already observed time points. This can be achieved using Brownian bridge simulation. As before, we can simply use the function sim(): all time points in a vector supplied to the function which fall between existing times in the existing bm path will be appropriately simulated. The handling of the function sim() intuitively deals with a mixed vector of times which require unconditioned and conditioned simulation. For example, suppling the vector \(t = (0.5, 1.5, 2.5, 3.5, 4.5)\) (and noting that the current end point of bm is \(t = 4\)), we have:

bm |> sim(c(0.5, 1.5, 2.5, 3.5, 4.5))
plot(bm)

Further note that the first argument to sim() can be a vector of times produced in any manner, say via seq() (or even by random simulation of times):

bm |> sim(runif(5, -1, 4.5))
plot(bm)

The above will simulate at 5 times randomly chosen between the start and end times of the existing path. Note this highlights that the vector of times passed to sim() need not be sorted: (non-standard) sorting of required times is handled by the package to ensure robustness of functionals.

A peek under the hood

At this point it is instructive to look at the stored path, which can be achieved by using the standard R element indexing syntax, as though the object was a data frame. Therefore,

  • bm[,"t"] contains the times of observations, \(t\);
  • bm[,"W_t"] contains the value of the path at those times, \(W_t\).
bm[,"t"]
#>  [1] -1.0000000 -0.2347582  0.0000000  0.5000000  1.0000000  1.5000000
#>  [7]  2.0000000  2.0786540  2.3973197  2.5000000  3.0000000  3.0303288
#> [13]  3.5000000  4.0000000  4.0194883  4.5000000
bm[,"W_t"]
#>  [1] 100.00000 100.20550 101.00000  99.61020  99.00000  99.39885  98.76083
#>  [8]  99.10028  99.51900  99.17654  99.43776  98.82608  98.05434  96.99743
#> [15]  97.05768  97.87486

In this simple setting \(t\) and \(W_t\) comprise the “path skeleton”: they are a finite dimensional representation of the Brownian motion path, which can be used to reconstruct the path at any other time point of interest. The indices of the vectors t and W_t match one another, and so (for instance) the value taken by the path at time \(t = 1\) is \(W_t = 99\). To avoid the user having to do complicated matching of the indices in order to extract the value of the path at a given time, one can instead index the bm object as follows, specifying in the t argument what time (or vector of times) they require, together with the part of the object they require (here, "W_t", which must be quoted):

bm[t = 1, "W_t"]
#> [1] 99

Note that by default calls to the bm environment without specifying the t argument will instead assume the user wants the associated index. For instance, the following retrieves the path location at the first time, whatever that may be (here, the first time is \(t=-1\)):

bm[1, "W_t"]
#> [1] 100

As an illustration of using this “under the hood” information, we can simulate an Exponential(\(\lambda = W_{4.5} / 10\)) time past the end of the path, without having had to track in our code that \(t=4.5\) is the end of the path what this current end time is:

bm |> sim(max(bm[,"t"]) + rexp(1, bm[t = 4.5, "W_t"]/10))
plot(bm)

Of course, the specification of the above time to simulate is rather clunky. To aid the user the package provides default labels to the time points which constitute the path skeleton (with each point possibly having multiple labels). We further give the user the option to add in their own labels. For instance, we can see the currently defined labels as a named list of times:

bm[,"labels"]
#> $start
#> [1] -1
#> 
#> $end
#> [1] 4.563118
#> 
#> $seg.start
#> [1] -1
#> 
#> $seg.end
#> [1] 4.563118
#> 
#> $forced
#> [1] -1  0  1
#> 
#> $user
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185

Note that each list item name corresponds to the times that have been assigned a particular label (and note that some times have multiple labels, for instance \(t=1\)). For our existing path bm all existing labels are those automatically assigned by the package. All automatically assigned labels are dynamically updated as the path changes. In particular:

  • start and end are the first and last times in the existing skeleton.
  • seg.start and seg.end are the first and last times for respective “segments” (see Segments).
  • forced is applied to any point introduced by the user by specifying not just the time, but also the position of the path. Recall, when we initialised the path (at \(t = -1\)), we introduced “forced” points at \(t = 0\) and \(t = 1\).
  • user is applied to any time for which the path has been simulated by the user using the sim() function.

To access the vector of times associated with a label then this can be achieved as follows (for instance, supposing we wanted the user times)

bm[,"labels"]$user
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185

Alternatively this can be achievied by simply passing a label (for instance, user) as the index in a call to bm and requesting the associated times t:

bm["user","t"]
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185

There is considerable flexibility for the user to introduce their own labelling. Existing times associated with a user specified label are fixed (and do not get dynamically updated in the same manner as the automatic labels above). This can be done as the path is generated by using the optionally available argument label in either the add.points() or sim() function we have introduced thus far. label if specified should be a vector of length \(1\) (in which case the same label is applied to all elements of the vector suppled), or of length equal to that of the supplied vector. For instance, suppose we wanted to simulate the path at times \(t = (6, 7)\) and give them the label special.time.1 and special.time.2 respectively then:

bm |> sim(c(6, 7), label = c("special.time.1", "special.time.2"))
bm[,"labels"]
#> $start
#> [1] -1
#> 
#> $end
#> [1] 7
#> 
#> $seg.start
#> [1] -1
#> 
#> $seg.end
#> [1] 7
#> 
#> $forced
#> [1] -1  0  1
#> 
#> $user
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185  6.0000000  7.0000000
#> 
#> $special.time.1
#> [1] 6
#> 
#> $special.time.2
#> [1] 7

Note now the addition of the new labels and associated times. We may wish to specify an additional user label to existing times with existing labels. This can be done by simply supplying a vector of the existing times to the function add.labels(). For instance, we may wish the times with the label forced to be be given the label special.time.1, and those with the label user to be given the additional label special.time.2:

bm |>
  add.labels(t = bm["forced","t"], label = c("special.time.1")) |>
  add.labels(t = bm["user","t"], label = c("special.time.2"))
bm[,"labels"]
#> $start
#> [1] -1
#> 
#> $end
#> [1] 7
#> 
#> $seg.start
#> [1] -1
#> 
#> $seg.end
#> [1] 7
#> 
#> $forced
#> [1] -1  0  1
#> 
#> $user
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185  6.0000000  7.0000000
#> 
#> $special.time.1
#> [1]  6 -1  0  1
#> 
#> $special.time.2
#>  [1]  7.0000000  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000
#>  [7]  2.5000000  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540
#> [13]  2.3973197 -0.2347582  4.5631185  6.0000000

Note that the above could be achieved by using the sim() function (noting that applying sim() to an existing path does not alter the existing path).

Of course, we may wish to delete labels: this can be done using the delete.labels() function. For instance, now removing the special.time.1 label, and removing the special.time.2 label from the times \(t = (7, 2, 3)\):

bm |>
  delete.labels(label = c("special.time.1")) |>
  delete.labels(t = c(7, 2, 3), label = c("special.time.2"))
bm[,"labels"]
#> $start
#> [1] -1
#> 
#> $end
#> [1] 7
#> 
#> $seg.start
#> [1] -1
#> 
#> $seg.end
#> [1] 7
#> 
#> $forced
#> [1] -1  0  1
#> 
#> $user
#>  [1]  2.0000000  3.0000000  4.0000000  4.5000000  3.5000000  2.5000000
#>  [7]  1.5000000  0.5000000  4.0194883  3.0303288  2.0786540  2.3973197
#> [13] -0.2347582  4.5631185  6.0000000  7.0000000
#> 
#> $special.time.2
#>  [1]  4.0000000  4.5000000  3.5000000  2.5000000  1.5000000  0.5000000
#>  [7]  4.0194883  3.0303288  2.0786540  2.3973197 -0.2347582  4.5631185
#> [13]  6.0000000

Note that it is only possible to add.labels and delete.labels for user labels, and not any automatic labels. For instance,

bm |> add.labels(t = 6, label = "end")
#> Error in `assert.bmlabel()`:
#> ! one of the specified labels clashes with an internal label naming scheme: label not added, please change.

Of course, we may also wish to know what value(s) the path took at time(s) given by a specific label. Again, this can be achievied by simply passing a label (for instance, special.time.2) as the index in a call to bm:

bm["special.time.2", "W_t"]
#>  [1]  96.99743  97.87486  98.05434  99.17654  99.39885  99.61020  97.05768
#>  [8]  98.82608  99.10028  99.51900 100.20550  97.52921  96.30520

Note the flexibility labelling offers, as we can use labels themselves as arguments in function calls. For instance, returning to our earlier example which motivated our labelling, we could simulate the path at an Exponential(\(\lambda = W_{4.5} / 10\)) time past the end of the path, without having had to track in our code that \(t=7\):

bm |> sim(bm["end", "t"] + rexp(1, bm["end", "W_t"]/10))
plot(bm)

Manipulating and existing skeleton

Having initialised the path at \(W_{-1}=100\) we may naturally wish to simulate a time prior to the time the path was intisialised. Simulating prior to the initialised time is ambiguous once the path also includes layers (see Layers, layers, everywhere!). As a result, directly using the sim() function (for instance, simulating the path at time \(t=-2\)) will result in an error:

bm |> sim(-2)
#> Error in `sim.BrownianMotion()`:
#> ! cannot simulate path at times before the path was initialised (at time -1) -- to add in points prior to the beginning of path use add.points()

As indicated in the error message, one apporach to handling this is to manually add in a point. For instance, if we were unconditionally simulating the path to the left of the initial time:

bm |> add.points(-2, rnorm(1, bm["start", "W_t"], sqrt(abs(-2 - bm["start", "t"]))))
plot(bm)

Of course, we may wish to now artificially “adjust” the path by shifting the path to the right so it is initialised at \(t=0\). This is easily achieved by creating a new path (say, bm.adjust) using the concat.bm() function and specifying the time the adjusted path should be initialised (although for full details of the concat.bm() function see Path Concatenation)

bm.adjust <- concat.bm(bm, t0 = 0)
plot(bm.adjust)

bm.adjust[,"labels"]
#> $start
#> [1] 0
#> 
#> $end
#> [1] 9.087926
#> 
#> $seg.start
#>   
#> 0 
#> 
#> $seg.end
#>          
#> 9.087926 
#> 
#> $forced
#> [1] 1 2 3 0
#> 
#> $user
#>  [1] 4.000000 5.000000 6.000000 6.500000 5.500000 4.500000 3.500000 2.500000
#>  [9] 6.019488 5.030329 4.078654 4.397320 1.765242 6.563118 8.000000 9.000000
#> [17] 9.087926 0.000000
#> 
#> $special.time.2
#>  [1] 6.000000 6.500000 5.500000 4.500000 3.500000 2.500000 6.019488 5.030329
#>  [9] 4.078654 4.397320 1.765242 6.563118 8.000000

Note that for the bm.adjust path all associated labels have also been shifted.

Cautionary note on bm objects

It is important for package users to understand that all Brownian motion path objects, such as bm and bm.adjust in the above, are internally stored as R environments. This means that attempting to assign from one variable name to another does not result in a copy of the object as one might expect with idomatic R: instead, you end up with two variable names referencing the same environment so that calls on one variable appear to update both. For example,

# Current end time of path in bm
bm["end","t"]
#> [1] 7.087926
# Assign bm to a new variable called bm.x
bm.x <- bm
# Check end time of bm.x
bm.x["end","t"]
#> [1] 7.087926
# Extend bm.x path
bm.x |> sim(bm.x["end","t"] + 100)
# Observe both have updated since these are just handles to the same environment
bm["end","t"]
#> [1] 107.0879
bm.x["end","t"]
#> [1] 107.0879

If you do in fact wish to create a copy rather than an additional handle to the same path, then a simple call to concat.bm() without providing any other path to concatenate with will trigger a copy:

# Current end time of path in bm
bm["end","t"]
#> [1] 107.0879
# Use concat to trigger a copy
bm.x <- concat.bm(bm)
# Check end time of bm.x
bm.x["end","t"]
#> [1] 107.0879
# Extend bm.x path
bm.x |> sim(bm.x["end","t"] + 100)
# Observe this time that both have not updated since we obtained a copy
bm["end","t"]
#> [1] 107.0879
bm.x["end","t"]
#> [1] 207.0879

Layers, layers, everywhere!

In many settings we want to simulate Brownian motion together with “layer” information. Layer information determines a compact (spatial) interval which constrains the path over a given time interval. This can be achieved using this package in three broad ways:

  • Bessel: This follows (broadly) the approach outlined in (Beskos et al. 2008).
  • Localised: This follows (broadly) the approach outlined in (Devroye 2009), following the implementation in (Pollock et al. 2020)
  • Intersection: This follows (broadly) the approach outlined in (Pollock et al. 2016),

There are merits to using each of these approaches, and so the user has the flexibility (when possible) of specifying these.

If a layer is simulated then this forms part of the path skeleton. Once a layer is simulated, then simulating further points of the path is more involved: again, these fall into the category of unconditional and conditional simulation.

Bessel

We begin by specifying that we want to simulate a Bessel layer over the interval \(t = [10,12]\) for our path bm.adjust. This can be achieved by using the (simply named!) function layers(). Note that \(t = 10\) is beyond the existing terminal time of the path, and so the path will first be (unconditionally) simulated to this time:

bm.adjust |> layers(10, 12)
plot(bm.adjust)

Purple is used to illustrate a Bessel layer (which will be useful later when we consider simultaneously different layer types). The solid lines are the layer indicating the spatial interval which constrains the path (the path falls entirely between the upper and lower solid lines). The dashed lines indicate a spatial interval which does not constrain the path (the path at some point will fall out-with the interval). As the end-points of the interval co-incide with the dashed lines we further know that the path will fall at some point both above the upper dashed line and below the lower dashed line.

Note we can now simulate at a time falling within this Bessel layer specification and the above constraints will be honoured and result in an updated layer specification. This is achieved using the same sim() function used to date: it will automatically detect the surrounding layer information and conduct conditional Bessel layer simulation. For example, simulating at \(t = 11\) gives:

bm.adjust |> sim(11)
plot(bm.adjust)

In particular, note that the Bessel layer specification has now updated, resulting in the single layer becoming two Bessel layers either side of the newly simulated time. In this case (though not always) this has resulted in two further variants of the Bessel layer. This is shown in the left hand interval by two dotted lines replacing the interior dashed lines of the earlier figure, and in the right hand interval by one dotted, one dashed lines replacing the interior dashed lines.

The dotted line indicates a more nuanced layer information, where the layer does not coincide with the end-points for the interval, but again indicates a spatial interval which does not constrain the path (the path at some point fall out-with the interval).

More precisely, in the left hand interval with two dotted lines there are three possible cases for the (unrealised) path. In particular, either:

  1. the path has a minimum which falls between the lower dotted and lower solid line, and, a maximum which falls between the upper dotted and upper solid line; or
  2. the path has a minimum which falls between the lower dotted and lower solid line, and, a maximum which falls between the maximum of the end-points of the path and the upper dotted line; or
  3. the path has a minimum which falls between the minimum of the end-points of the path and the lower dotted line, and, a maximum which falls between the upper dotted and upper solid line.

The right hand interval has both a dashed (lower) and dotted (upper) line. These follow the same convention as the two earlier examples. As one of the end-points co-incides with the dashed line (lower) we know the path at some point falls below this level. This eliminates the possibility of case 3 in the list above, leaving case 1 and 2. It is possible that the dashed and dotted lines were in the reverse order (dashed being the upper, and dotted being the lower), in which case we would either have case 1 or 3.

When using layers() there is no need for the start and end times to co-incide with the current skeleton. It can be some future interval (as above), an interval of time pre-ceding the current skeleton, or an interval which partially includes the current skeleton. If there are pre-existing layers then a call to layers() only introduces layers for intervals which do not already have layers. For instance, consider the introduction of layers over the interval \(t\in[5,14]\):

bm.adjust |> layers(5, 14)
plot(bm.adjust)

Note that as the bm.adjust object becomes larger, calling the plot command can be less informative. You may instead wish to only plot a specific time interval for clarity. This is possible using the t.lim argument within the plot command. For instance, suppose you were interested only in visualising the interval \([5,10]\):

plot(bm.adjust, t.lim = c(5, 10))

Localised

First passage times

The second layer type we will consider is based around first passage times of Brownian motion.

First passage times are the basis of the “localised layers” we introduce fully in the next subsection. Here, we consider the path initialiased at its current terminal value, and determine the first time it leaves a prescribed interval. The interval can be prescribed by the user in a number of different ways, by specifying among a choice of arguments:

  • delta: specifying a symmetric interval around the current terminal value;
  • delta.l and delta.u: specifying a asymmetric interval which is delta.l below and delta.u above the current terminal value;
  • l and u: specifying the lower and upper boundary values directly.

To better illustrate this we begin by simulating the existing path unconditionally at time \(t=15\) by again calling sim() before then simulating a first passage to a symmetric interval \(\pm 1\) from \(W_{10}\). Notice that commands can be chained together using the pipe |> operator:

bm.adjust |>
  sim(15) |>
  first.passage(delta = 1)
plot(bm.adjust)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the BrownianMotion package.
#>   Please report the issue to the authors.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

The green lines indicate the bounds specified by You, the user. The green line does not form part of the realised Brownian motion, and is used for visualisation. As such, it can be turned off using the hide.user argument in the plot command:

plot(bm.adjust, hide.user = TRUE)

Note that the green lines were over-plotting solid red lines. Again, solid lines form part of the Brownian motion object and indicate the interval which constrains the path over the specified interval. We term this a “localised layer” as the end-point coincides with layer. Red is used as the colour to plot such layer types. The dashed lines follow the earlier convention, although note there will only be one. In the case above we in effect have the upper dashed line coinciding with the left hand end-point, and both the dashed and solid lines coinciding with the right hand end-point.

If we instead we want asymmetric layers and use the arguments delta.l and delta.u then these are the lower and upper offsets (respectively) of the current terminal end-point:

bm.adjust |> first.passage(delta.l = -1, delta.u = 3)
plot(bm.adjust)

The output in this case is slightly different as this is achieved algorithmically by decomposing the requested layers into iterative calls of the symmetric variant of the function.

We may not be interested in the first passage time of a level offset by the current terminal time, but instead some specified level. This can be achieved by using the arguments l and u where these are the lower and upper layers requested. For instance, we may want the first passage time of the (spatial) interval \([90,100]\):

bm.adjust |> first.passage(l = 90, u = 100)
plot(bm.adjust)

Note that if using the arguments l and u then the current terminal point must be contained within the interval, as otherwise it is not possible and results in an error:

bm.adjust |> first.passage(l = 80, u = 85)
#> Error in `first.passage()`:
#> ! u must be greater than the current endpoint state of the Brownian motion.

Again, time points within a localised layer can be simulated using sim() in the usual manner (and the package will automatically detect the surrounding layer information to perform the correct conditioning).

bm.adjust |> sim(c(40, 60))
plot(bm.adjust)

Note conditional simulations within a localised layer can give rise to a change in layer type, in this instance resulting in the first part of the localised layer becoming of Bessel type (recalling that for localised layers the layer coincides with the the end-point).

Further note that first passage times can only be requested for one-dimensional Brownian motion. For multi-dimensional Brownian motion it is more natural to consider the approach outlined in the next section:

Localised layer over a specified interval

Simulating localised layers by means of first passage times can be advantageous in some settings, although the time the first passage occurs is random, and this randomness can introduce its own complications.

Instead, we may wish to simulate localised layers over a specified interval (in the same manner as we did for Bessel layers). For instance, suppose we want a localised layer in the interval \(t\in[23,29]\) this can again be achieved using the layers() function, but with the additional argument type = "localised":

bm.adjust |> layers(90, 140, type = "localised")
plot(bm.adjust)

Note that the resulting skeleton for the interval will comprise zero or more (red) localised layers, followed by a single (purple) Bessel layer: this is because the end-point will no longer coincide with the layer. By default delta is chosen to be \(\sqrt{t-s}/2\), although this can be altered by the user using the argument mult which multiplies the default layer size. For instance,

bm.adjust |> layers(180, 210, type = "localised", mult = 5)
plot(bm.adjust)

We may wish to add a bridge analogue of a localised layer: ie a localised layer over one or more existing intervals in the skeleton. This can be achieved again by using the layers() function with the argument type = "localised". For instance, we may wish to introduce localised intervals over every interval between times \([140,230]\) (where \(t=230\) exceeds the current end-time of \(t=210\)):

bm.adjust |> layers(140, 230, type = "localised")
plot(bm.adjust)

Note that the path has been extended (unconditionally) beyond \(t=140\) as per above. For each interval within the current skeleton, which did not already have a layer, a localised layer has been introduced. The “bridge” localised layers are simulated by means of an appropriate scaling of unconditional Brownian motions, and so have a gradient of \((W_t-W_s)/(t-s)\). The interior dashed, dotted and solid lines, together with the associated Bessel layer all follow the earlier conventions (albeit now with a gradient). The outer horizontal solid lines are for the use of You (and don’t form part of the skeleton) to easily indicate a spatial interval which constrains the path over a given (time) interval. We can again more clearly see the localised bridge layers by specifying a shorter interval to plot the path (for instance \([140,180]\)):

plot(bm.adjust, t.lim = c(140, 180))

As the horizontal layers do not strictly form part of the path skeleton, these can be toggled off when plotting by again using the hide.user argument:

plot(bm.adjust, t.lim = c(140, 180), hide.user = TRUE)

(Another) peek under the hood

In addition to the earlier discussion on the slots $t and $W_t, there are slots which contain the layering information, both for the main path and those in the transformed space just discussed above.

  • bm[,"layers"] is a data frame (tibble) containing the layering information of the main path, including variables:

    • type to indicate the layer type, one of:

      • bessel
      • localised
      • intersection (see forthcoming section);
      • bessel-bb
      • localised-bb
      • intersection-bb (see forthcoming section);

      The suffix -bb indicates that this layer type results from simulating the localised layer bridge of the last section.

    • t.l and t.u giving the temporal interval (left and right time points) to which the layer applies;

    • Ld corresponds to the lower solid line;

    • Uu corresponds to the upper solid line;

    • Lu corresponds to the lower interior dashed/dotted line;

    • Ud corresponds to the upper interior dashed/dotted line;

    • Lu.hard indicating whether the lower interior line coincides with the end point;

    • Ud.hard indicating whether the upper interior line coincides with the end point;

  • Note: The output in multiple dimensions will vary, as discussed in Multiple Dimensions.

  • $user.layers is a data frame containing the information on which first passage layers were requested by You (resulting in the green lines on the plot). This contains

    • t.l and t.u giving the temporal interval (left and right time points) to which the layer applies;
    • L and U which are the boundaries specified by You (perhaps indirectly via delta or delta.l/delta.u)

  • $bb.local is a nested Brownian motion object used in the construction of the localised layer bridges. Note this will only contain intervals pertaining to the localised layer bridge intervals of $layers (those ending -bb), but otherwise follows the same form.

  • Note: $user.layers and $bb.local are only specified for one-dimensional Browian motion

For example, the path produced to date in this vignette has the following internal layering representation:

bm.adjust[,"layers"]
bm.adjust$user.layers
t.l t.u L U
15.00000 16.29056 92.54972 94.54972
16.29056 23.91984 91.54972 95.54972
23.91984 80.43416 90.00000 100.00000 
bm.adjust$bb.local[,"layers"]
type t.l t.u Ld Uu Lu Ud Lu.hard Ud.hard
localised 140.0000 144.6401 -3.162278 3.162278 0.000000 3.162278 TRUE TRUE
bessel 144.6401 180.0000 0.000000 6.324555 3.162278 5.341517 TRUE TRUE

Again, we might naturally be more interested in the layers associated with a specific time. A given time may not lie within the current skeleton of the path, and if it does (for reasons discussed later), then of the two possible layers (to the left and right of that point) we are naturally more interested in the right hand layer. Note that there may be no layer associated with a given time (in which case it is supressed). For instance, suppose we were interested with path layers at the times \(t = \{4.5, 4.6 \dots, 5.5\}\):

bm.adjust[t = seq(4.5, 5.5, by = 0.1), "layers"]
#> # A tibble: 3 × 9
#>   type     t.l   t.u    Ld    Uu    Lu    Ud Lu.hard Ud.hard
#>   <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>   <lgl>  
#> 1 bessel  5     5.03  98.7  99.6  98.8  99.4 TRUE    TRUE   
#> 2 bessel  5.03  5.5   97.4  99.5  98.1  98.8 TRUE    TRUE   
#> 3 bessel  5.5   6     96.3  98.8  97.0  98.1 TRUE    TRUE

Intersection

Now suppose we want an Intersection layer over the interval \(t = [50,60]\), this can be done by calling the function layers() with the argument type = "intersection". The structure of introducing Intersection layers follows that of Bessel layers. Note that \(t = 280\) is beyond the current terminal point of the path:

bm.adjust |> layers(280, 350, type = "intersection")
plot(bm.adjust)

Blue is used to illustrate an Intersection layer. Again, solid lines are the layer indicating the spatial interval which constrains the path. The interior lines for Intersection layers are always dashed. Dashed lines have the same interpretation as for Bessel layers, but need not coincide with the end-points for an interval. More precisely, we know the path has a maximum which falls between the upper dashed and solid lines, and a minimum which falls between the lower dashed and solid lines.

The sim() function can be used once again to simulate intermediary points of intersection layers. For instance, consider simulating the path at \(t=325\):

bm.adjust |> sim(325)
plot(bm.adjust)

Note that on simulating an intermediary point(s) a number of Bessel (purple) or Intersection layers may arise. Intersection layers have the same interpretation as the parent Intersection layer, and Bessel layers are the same as in the earlier section. By default, the package will if possible attempt to “convert” Intersection layers to Bessel layers: we will address this default in the next subsection. Any further simulation will be handled appropriately (for the given layer type) by using the sim() function.

(Yet Another) peek under the hood

As demonstrated in the previous section, by default the package has a preference for using Bessel layers over Intersection layers wherever possible. It is typical that the user will want to only consider one layer type, and the Bessel layer type has been selected as default as this is most likely to be the one practitioners will use. We do however want the flexibility to use multiple layer types, or change the default to Intersection.

To do this we can change the default argument from prefer = "bessel" to prefer = "intersection" within the functions we have considered. For instance, we may consider a new path (bm2) initialised at \(W_0=0\):

set.seed(42)
bm2 <- create.bm(prefer = "intersection")

This has the consequence that any call to sim() (with no further argument as the preference is inherited from the creation of the path) will where possible convert Bessel layers to Intersection layers. Note that although it is only occasionally possible to convert Intersection layers to Bessel layers, it is always possible to change Bessel layers to Intersection layers.

Now suppose we want to simulate Bessel layers over \(t\in[1,3]\). This can be achieved by again using the layers() function, but noting that as the default type is “intersection” (as inherited on the creation of the path with the argument prefer = "interesection"), the additional argument type = "bessel" has to be passed to the layers() function to simulate a Bessel layer, In particular,

bm2 |> layers(1, 3, type = "bessel")
plot(bm2)

then if we want to simulate the path at \(t = 2\) and have resultant Bessel layers then this needs to be specified by the user as it differs from the default on the paths creation (prefer = "intersection"):

bm2 |> sim(2, prefer = "bessel")
plot(bm2)

If however we don’t include the argument prefer = "bessel", then further simulation of the path (say at \(t=1.25\) and \(t=2.5\)) results in both Bessel layers being converted to Intersection layers:

bm2 |> sim(c(1.25, 2.5))
plot(bm2)

Note that we do not need to specify we have a preference for Intersection layers, as this preference is an attribute assigned to the path on its creation using create.bm().

It is not possible to convert from, or to, Localised layers from other layer types. However, once the path has been created with the attribute prefer = "intersection" then conditional simulation of localised layers will result in the appearance of Intersection layers (as opposed to Bessel layers). For instance, note that when simulating a Localised layer over the interval \(t\in[4,5]\) that the final interval is an Intersection layer (and not a Bessel layer as before):

bm2 |> layers(4, 5, type = "localised")
plot(bm2)

The preference for Intersection layers inherited from the creation of the path can be changed when calling any of the other functions we have introduce. For instance, consider simulating localised layers (of the Bessel sort) in the interval \(t\in[6,8]\):

bm2 |> layers(6, 8, type = "localised", prefer = "bessel")
plot(bm2)

and then the simulation of intermediary points in this interval:

bm2 |> sim(c(6.5, 7, 7.5), prefer = "bessel")
plot(bm2)

We can change the layer type (if possible) for every interval in the path by simply calling layers() together with arguments for the end-points for the path. For instance, suppose we now wanted Intersection layers (we choose to add type = "intersection" for clarity, but it is strictly not needed since this is the prefer default set upon the creation of the path):

bm2 |> layers(0, 8, type = "intersection")
plot(bm2)

This results in any Bessel layers being converted to Intersection layers. Similarly we could attempt to convert back to Bessel layers (although this may not always be possible for every interval):

bm2 |> layers(0, 8, type = "bessel")
plot(bm2)

Note that converting the layers from Intersection to Bessel or Bessel to Intersection does not change the preference when calling the sim() function. For instance:

bm2 |> sim(3.75)
plot(bm2)

Refinement

When simulating layers the package will by default control the size of the layer so it scales appropriately with the length of the interval being considered. In particular, by default, the package will attempt to ensure the difference between the upper solid and upper dashed/dotted lines AND the difference between the lower solid and upper dashed/dotted lines is no greater than \(\sqrt{t-s}\) (where \(s\) and \(t\) are the left and right hand end-points of the interval respectively). This control on the layer size can be removed by changing the refine = TRUE attribute to refine = FALSE on creation of the path. Note that this is not recommended as it can introduce numerical instability, but to show its affect consider the following Bessel layer over the interval \([0,10]\):

set.seed(42^2)
bm3 <- create.bm(refine = FALSE)
bm3 |> layers(0, 10)
plot(bm3)

Now, if we simulate a number of intermediary points (say a grid of \(20\)), then note that although we again have Bessel layers that the size of the Bessel layers is “stretched” when contrasted with earlier simulations:

bm3 |> sim(seq(0, 10, by = 0.5))
plot(bm3)

Contrast this with a similar path where we haven’t set refine = FALSE:

bm4 <- create.bm()
bm4 |>
  layers(0, 10) |>
  sim(seq(0, 10, by = 0.5))
plot(bm4)

Whether the package does refinement by default is inherited from the attribute when the path is created. We can change this at some future point by appropriately specifying within the other functionals in the same manner as the prefer attribute discussed earlier. For instance, consider the introduction of mid-points to the existing path in the interval \([0,5]\):

bm3 |> sim(seq(0, 5, by = 0.25), refine = TRUE)
plot(bm3)

Note that the intervals in the first half of the path have been refined, but the second half of the path remains the same.

Now, in the case where we are refining by default (refine = TRUE), then we may want to alter the default layer size. This is done using the mult attribute, which was described earlier with Localised layers. Recall that by default when refining, the difference between the upper solid and upper dashed/dotted lines AND the difference between the lower solid and upper dashed/dotted lines is no greater than \(\sqrt{t-s}\). The mult attribute multiplies this tolerance by a user-specified amount. By default, mult = 1 and is specified at the creation of the path. For instance, consider a new path (bm5) together with Bessel layers in the interval \([0,10]\), but with the default layer size multiplied by \(3\) :

bm5 <- create.bm(mult = 3)
bm5 |> layers(0, 10)
plot(bm5)

The simulation of the path at additional points will inherit mult = 3. For instance, contrast the following with bm4 above:

bm5 |> sim(seq(0, 10, by = 0.5))
plot(bm5)

If we introduce an additional points (say mid-points to the existing path in the interval \([0,5]\)), but we want increased increased refinement (say of mult = 1) this can be done by changing the default argument in any of the earlier functions as before. For instance:

bm5 |> sim(seq(0, 5, by = 0.25), mult = 1)
plot(bm5)

Note that the intervals in the first half of the path are far more resolved. It is possible to set the layer size to some constant level (not dependent on the length of the interval) by specifying mult to be at a level of constant multiplied by the reciprocal of \(\sqrt{t-s}\). This however is not recommended, as it can introduce numerical instability.

Finally, we may wish to refine the existing intervals without the introduction of additional points. This can be achieved using the refine function, which can increase the resolution on any (or all) intervals. For instance:

bm5 |> refine(5, 10, mult = 1)
plot(bm5)

The mult = 1 argument has increased the resolution of existing layers in the interval \([5,10]\) by \(3\times\), and will be comparable to the first half of the path.

Deletion

It is possible that the user wants to delete both the path observations and layers, or only the layers, of the skeleton over some specified interval. This can be achieved within the package by using the delete.skeleton() function, and specifying an open interval which you wish to delete from the skeleton. The left and right hand ends of the open interval are specified using the arguments l and r respectively. For instance, considering bm5, we may wish to delete the entire skeleton in the interval \((2.5,5)\):

bm5 |> delete.skeleton(l = 2.5, r = 5)
plot(bm5)

By default the package will delete both path observations and layers within the open interval (by default the argument type = "all"). However, deletion of only the layers (retaining the path observations) is possible by setting type = "layer". For instance, considering deleting only the skeleton layers in the interval \((5,7.5)\):

bm5 |> delete.skeleton(l = 5, r = 7.5, type = "layer")
plot(bm5)

It is also possible to delete all user labels (mirroring the functionality of delete.labels discussed earlier, automatic labels are not deleted). For instance, when running the following code, note that the automatically generated labels between l and r are retained:

bm5 |> delete.skeleton(l = 5, r = 7.5, type = "label")
bm5[,"labels"]
#> $start
#> [1] 0
#> 
#> $end
#> [1] 10
#> 
#> $seg.start
#>   
#> 0 
#> 
#> $seg.end
#>    
#> 10 
#> 
#> $forced
#> [1] 0
#> 
#> $user
#>  [1] 10.00  5.00  7.50  2.50  9.00  9.50  6.50  7.00  1.50  2.00  8.50  8.00
#> [13]  6.00  5.50  1.00  0.50  2.25  1.75  1.25  0.75  0.25

Note that as a safety precaution, l and r must co-incide with time points within the skeleton. As a consequence, if You want to delete an interval which does not co-incide with time points within the skeleton then You must first introduce those points using the sim function. For instance, consider deleting the skeleton entirely in the interval \((2.0001,7.9999)\) (noting these times do not currently appear in the skeleton) then this can be achieved as follows:

bm5 |>
  sim(c(2.0001, 7.9999)) |>
  delete.skeleton(l = 2.0001, r = 7.9999)
plot(bm5)

In some scenarios we may wish to delete the path entirely to the left or right of some point. This is possible by not specifying the right or left point (respectively). For instance, suppose we wish to delete the skeleton entirely to the right of \(t = 2\) (fixing the path to the left of \(t = 2\)) then:

bm5 |> delete.skeleton(l = 2)
plot(bm5)

Alternatively, deletion entirely to the left of \(t=1\) is possible in a similar manner by not specifying the l argument:

bm5 |> delete.skeleton(r = 1)
plot(bm5)

Note that (as a fail-safe) it is not possible to delete the entire path by leaving both arguments un-specified. Intead You should do this by removing the bm5 object from your workspace using the base R rm command. If you attempt to do this You will get an (informative) error:

bm5 |> delete.skeleton()
#> Error in `delete.skeleton.BrownianMotion()`:
#> ! at least one of l or r must be specified

Segments

For many practical uses of the package, the user may wish to add in path discontinuities. In such settings we will have a number of path “segments” (the term we give for the continuous paths separated by the path discontinuities). To aid this, in addition to "t" and "W_t", the stored path actually contains a "W_tm" object:

  • bm[,"W_tm"] contains the value of the path at those times, \(W_{t-}\).

Note that when considering the simulation of any intermediary time for an existing path (using the sim() function), we treat the process as cadlag: for instance, consider simulation at time \(t'\) for an existing path (using the sim() function), where \(t_1 < t' < t_2\) and \(t_1\) and \(t_2\) are the closest realisations of the path on either side, then the simulation will be done conditional on \((t_1, W_{t_1})\) and \((t_2, W_{t_2-})\).

To illustrate "W_tm", we begin by considering a new Brownian motion object, bm6, realised at a collection of times:

bm6 <- create.bm()
bm6 |> sim(seq(0, 10, length = 11))
plot(bm6)

Note that if we look at the associated $t, $W_t and $W_tm for the bm6 Brownian motion object, that (as suggested by our plot above) we have a single continuous path (\(W_t = W_tm\) at all times the path has been realised):

data.frame(t        = bm6[,"t"],
           W_t      = bm6[,"W_t"],
           W_tm     = bm6[,"W_tm"],
           equality = bm6[,"W_t"] == bm6[,"W_tm"])
#>     t         W_t        W_tm equality
#> 1   0  0.00000000  0.00000000     TRUE
#> 2   1 -0.60461072 -0.60461072     TRUE
#> 3   2  0.15414496  0.15414496     TRUE
#> 4   3 -0.02528936 -0.02528936     TRUE
#> 5   4 -0.41498855 -0.41498855     TRUE
#> 6   5 -1.38261959 -1.38261959     TRUE
#> 7   6 -0.78879135 -0.78879135     TRUE
#> 8   7 -0.81582735 -0.81582735     TRUE
#> 9   8 -1.55678461 -1.55678461     TRUE
#> 10  9 -1.32747495 -1.32747495     TRUE
#> 11 10 -0.12975571 -0.12975571     TRUE

Now, to add in a path discontinuity we can use the add.segment() function, which moves by some user specified amount the path between two times which appear in the path skeleton , or the entire path to the right of a given time, or the entire path to the left of a given time. There are a number of ways the user can specifiy the movement.

For instance, suppose we wish to move the path between times in the interval \(t\in[1,3)\) by a constant \(1\), then this can be achieved using the optional delta argument. We may wish to label this and so we can additionally pass in our own labelling argument:

bm6 |> add.segment(l = 1, r = 3, delta = 1, label = "path.move")
plot(bm6)

data.frame(t        = bm6[,"t"],
           W_t      = bm6[,"W_t"],
           W_tm     = bm6[,"W_tm"],
           equality = bm6[,"W_t"] == bm6[,"W_tm"])
#>     t         W_t       W_tm equality
#> 1   0  0.00000000  0.0000000     TRUE
#> 2   1  0.39538928 -0.6046107    FALSE
#> 3   2  1.15414496  1.1541450     TRUE
#> 4   3 -0.02528936  0.9747106    FALSE
#> 5   4 -0.41498855 -0.4149885     TRUE
#> 6   5 -1.38261959 -1.3826196     TRUE
#> 7   6 -0.78879135 -0.7887913     TRUE
#> 8   7 -0.81582735 -0.8158274     TRUE
#> 9   8 -1.55678461 -1.5567846     TRUE
#> 10  9 -1.32747495 -1.3274749     TRUE
#> 11 10 -0.12975571 -0.1297557     TRUE

Note that the plot of the path has changed. Vertical, dotted, black lines at \(t = 1\) and \(t = 3\) indicate there are path discontinues at those times. As per convention with cadlag processes, at \(t = 1\) (for instance) there is an open circle indicating the value the path take at the time (\(W_{1-}\)), and (as per the package convention) a block dot at \(t = 1\) represented the value the path at \(W_1\). This can be verified in our data frame. As discussed above, conditional simulation of the path at intermediary times is appropriately handled by the package. For instance, consider a sequence of times between \(t = 0\) and \(t = 4\):

bm6 |> sim(seq(0, 4, length = 21))
plot(bm6)

Further note that the automatic labels we introduced for user convenience are now additionally identifying the segment start and end times (seg.start and seg.end respectively), and those appear in addition to the labels we’ve passed into the add.segment() function:

bm6[,"labels"]
#> $start
#> [1] 0
#> 
#> $end
#> [1] 10
#> 
#> $seg.start
#> [1] 0 1 3
#> 
#> $seg.end
#> [1]  1  3 10
#> 
#> $forced
#> [1] 0
#> 
#> $user
#>  [1]  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0  3.2  3.4  3.6  3.8  2.2
#> [16]  2.4  2.6  2.8  1.2  1.4  1.6  1.8  0.2  0.4  0.6  0.8
#> 
#> $path.move
#> [1] 1 3

Of course, shifting the path by some fixed amount may not be natural in all settings. Suppose instead we want to at \(t = 10\) “jump” to some given point (for instance, \(W_{10} = 0\)), then this can again be done by “moving” the path fixed at the left by the current end point (i,e. \([10,\infty)\)), to \(0\), by means of the W_t optional argument in add.segment(). To better illustrate this, following our “jump” to \(0\) we will simulate the path to \(t = 11\):

bm6 |>
  add.segment(l = 10, W_t = 0) |>
  sim(11)
plot(bm6)

Using only the argument r we can shift the path fixed at some right point in the same fashion. For instance, everything in the interval \(t\in[0,5)\) being shifted such that \(W_0 = 4\):

bm6 |> add.segment(r = 5, W_t = 4) 
plot(bm6)

A cautionary note is that if we shift the entire path by specifying l and r to be the start and end of the path respectively, then that will introduce path discontinuities at those time since the specification of l and r is of a half-open interval (ie \(t = 0-\) and \(t = 11\) are both times that appear in the existing path skeleton, and we would be moving the path in the interval \(t\in[0,11)\)). For instance,

bm6 |> add.segment(l = bm6["start", "t"], r = bm6["end", "t"], delta = -1)
plot(bm6)

Instead, we should omit both the l and r arguments (we will first undo the previous path move):

bm6 |>
  add.segment(l = bm6["start", "t"], r = bm6["end", "t"], delta = 1) |>
  add.segment(delta = -1)
plot(bm6)

Further note that the syntax of the delete.skeleton() function matches that of add.segment(). When dealing with path segments aligned with the specified l and r in then delete.skeleton() function then everything in the interval \([l, r)\) will be deleted: \(W_{l}\) will be deleted, and within the bm object \(W_{l-}\) will be copied so \(W_{l-}=W_l\); \(W_{r-}\) will be deleted and replaced with \(W_{r}\) within the bm object. For instance, lets consider removing the path segment we introduced between \(t = 1\) and \(t = 3\):

bm6 |> delete.skeleton(l = 1, r = 3)
plot(bm6)

data.frame(t        = bm6[,"t"],
           W_t      = bm6[,"W_t"],
           W_tm     = bm6[,"W_tm"],
           equality = bm6[,"W_t"] == bm6[,"W_tm"])
#>       t        W_t       W_tm equality
#> 1   0.0  3.0000000  3.0000000     TRUE
#> 2   0.2  3.0927744  3.0927744     TRUE
#> 3   0.4  3.2950951  3.2950951     TRUE
#> 4   0.6  2.4342801  2.4342801     TRUE
#> 5   0.8  2.4721542  2.4721542     TRUE
#> 6   1.0  2.3953893  2.3953893     TRUE
#> 7   3.0  2.9747106  2.9747106     TRUE
#> 8   3.2  2.9679436  2.9679436     TRUE
#> 9   3.4  2.6541898  2.6541898     TRUE
#> 10  3.6  2.7234927  2.7234927     TRUE
#> 11  3.8  3.2457690  3.2457690     TRUE
#> 12  4.0  2.5850115  2.5850115     TRUE
#> 13  5.0 -2.3826196  1.6173804    FALSE
#> 14  6.0 -1.7887913 -1.7887913     TRUE
#> 15  7.0 -1.8158274 -1.8158274     TRUE
#> 16  8.0 -2.5567846 -2.5567846     TRUE
#> 17  9.0 -2.3274749 -2.3274749     TRUE
#> 18 10.0 -1.0000000 -1.1297557    FALSE
#> 19 11.0  0.8149376  0.8149376     TRUE

Note that (as in the above example) the deletion of part of a path may result in a change to the number of segments: this will result in an updating of the associated automatic labelling (but not the user specified labelling).

bm6[,"labels"]
#> $start
#> [1] 0
#> 
#> $end
#> [1] 11
#> 
#> $seg.start
#> [1]  0  5 10
#> 
#> $seg.end
#> [1]  5 10 11
#> 
#> $forced
#> [1] 0
#> 
#> $user
#>  [1]  1.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0  3.2  3.4  3.6  3.8  0.2  0.4
#> [16]  0.6  0.8 11.0
#> 
#> $path.move
#> [1] 1 3

We can access the vector associated with W_tm by calling the environment in the same way as previously. Alternatively, by not specifying any argument then calling the environment returns t, W_t, W_tm and layers:

bm6[t = 1]
#> $t
#> [1] 1
#> 
#> $W_t
#> [1] 2.395389
#> 
#> $W_tm
#> [1] 2.395389
#> 
#> $layers
#> # A tibble: 0 × 9
#> # ℹ 9 variables: type <fct>, t.l <dbl>, t.u <dbl>, Ld <dbl>, Uu <dbl>,
#> #   Lu <dbl>, Ud <dbl>, Lu.hard <lgl>, Ud.hard <lgl>
#> 
#> $labels
#> $labels$user
#> [1] 1
#> 
#> $labels$path.move
#> [1] 1

Path Concatenation

In some settings it may be useful to concatenate path skeletons. Seperate path skeletons can be considered as segments of a newly created Brownian motion object, and if the paths are align then be considered as single segments. This can be done using the bm.concat() function that we introduced earlier: but in addition to passing in a single bm together with the required t0 argument specifying the start of the newly created path, any number of paths in the order they are to be concatenated. For instance, consider the concatenation of bm6 together with bm5, and then again bm6, into a newly created bm.concat object initialised at \(t = 10\)

bm.concat <- concat.bm(bm6, bm5, bm6, t0 = 10)
#> Warning in concat.bm.BrownianMotion(bm6, bm5, bm6, t0 = 10): Not all
#> BrownianMotion objects have the same mult parameter ... using setting of first
#> supplied object
plot(bm.concat)

Multiple Dimensions

Creation

Creating a multi-dimensional Brownian motion object is again accomplished using the create.bm() function, by specifying a number of different (optional) arguments. In addition to time t (defaulting to \(t=0\)) and position W_t (defaulting to vector \(0\)), we have dim which specifies the number of dimensions (default dim \(=1\)) and covariance matrix cov (defaulting to the identity matrix). For instance, a standard three-dimensional Brownian motion may be created as follows:

bm.d <- create.bm(dim=3)
bm.d$dim
#> [1] 3
bm.d[t=0, "W_t"]
#>      [,1] [,2] [,3]
#> [1,]    0    0    0

Some care has to be made in specifying W_t in multiple dimensions. If the specified vector is of length one, then that position is used to specify all coordinates of the path at the initial time. For instance,

bm.d <- create.bm(W_t = 3, dim=3)
bm.d[t=0, "W_t"]
#>      [,1] [,2] [,3]
#> [1,]    3    3    3

Whereas to initialize the path at the position \((0,1,2)\) we simply provide a vector to the W_t argument:

bm.d <- create.bm(W_t = c(0,1,2), dim=3)
bm.d[t=0, "W_t"]
#>      [,1] [,2] [,3]
#> [1,]    0    1    2

As before, we can specify on the path creation multiple points in the path. Here we supply a vector to t and a matrix of dimension the length of t times dim. For instance, consider the path at \(t=1\) at \(W_1=(0,1,2)\) and at \(t=2\) at \(W_2 = (2,0,1)\)

bm.d <- create.bm(t = c(1,2), W_t = matrix(c(0,1,2,2,0,1),nrow = 2, ncol = 3, byrow = TRUE), dim=3)
bm.d[t=c(1,2), "W_t"]
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    2    0    1

The function add.points can be used directly in multiple dimensions. For instance, if we want to add the point \(W_0 = (-1,-1,-1)\).:

bm.d |> add.points(t = 0, W_t = c(-1, -1, -1))
bm.d[t=c(0,1,2), "W_t"]
#>      [,1] [,2] [,3]
#> [1,]   -1   -1   -1
#> [2,]    0    1    2
#> [3,]    2    0    1

Similar care has to be taken in specifiying the covariance matrix cov. cov if specified as a scalar is interpreted as a multiplier of the identity matrix, if it is specified as a vector of length dim then it is interpreted as diagonal with the vector as the diagonal elements, or it can be specified as a matrix. For instance,

bm.d <- create.bm(dim = 3, cov = (diag(x= 0.5, nrow=3) + 0.5))
bm.d$cov
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.5  0.5
#> [2,]  0.5  1.0  0.5
#> [3,]  0.5  0.5  1.0

The final arguments in create.bm – namely refine, mult and prefer – have precisely the same interpretation as in the one dimensional setting (with the same defaults). The one difference is they can alternatively be specified as vectors of length dim if You wish different specification in each dimension.

Other Functions

Having created a BrownianMotionNd object, simulation of the path and layers is possible using exactly the same syntax as for the sim and layers functions before. For instance, suppose we simulate the path at times \(t=(0.5,5)\) and add a layers over the interval \(t=[0,10]\) then this is achieved simply by:

bm.d |> 
  sim(c(0.5,5)) |>
  layers(0,10)

Note that additional arguments for layers operate in exactly the same fashion in multiple dimensions. Note however first.passage is not possible in multiple dimensions as it does not have the same clear interpretation as in one dimension (although "localised" can be specified for the argument type in layers).

At this point we can see a plot of our nascent path, by simply using the plot() function:

plot(bm.d)

In multiple dimensions the plotting is somewhat different. By default a ‘pairwise’ plot is produced: each component is plotted along the diagonal and labelled appropriately. Note in the diagonal plots there only exists “black” layers in addition to the plotted path: this indicates intervals in which the entire path is constrained, but due to the covariance structure of the path, this does not have the more nuanced interpretation possible in the one-dimensional case (this is discussed more in the later plotting section). The off-diagonal plot on row \(i\) column \(j\) simply shows coordinate \(i\) plotted against coordinate \(j\) (with no axis for time). The red lines demarcate a region of the two dimensional space that the entire path is constrained in for the entire interval. Whereas the black demarcates a region parallel to the axis which constrains the red region. Again, more discussion on plotting in multiple dimensions is given in the later plotting section.

Refinement and Deletion also use the same functions (refine and delete.skeleton respectively) and syntax as in one dimension. For instance, suppose we wish to refine the entire path using the mult argument with mult \(=0.1\), and delete the entire path to the right of \(t=5\), then we can do this as follows:

bm.d |> 
  refine(s=0, t= 10, mult = 0.1) |>
  delete.skeleton(l = 5)
plot(bm.d)

A little more care has to be taken in adding in segments using the add.segment function. Although the arguments are the same as in one dimension, specifying the argument W_t instead follows the same conventions as in in the create.bm function discussed earlier. In particular, if W_t is a scalar then then point is repeated in every component, else a vector has to be specified of length dim. For instance, consider adding a new segment at \(t=5\) by jumping to the point \(W_5 = (10,11,12)\) can be achieved as follows (extended to \(t=10\) with layers to more clearly show the discontinuity):

bm.d |> 
  add.segment(l=5, W_t = c(12,10,10)) |>
  layers(5, 10)
plot(bm.d)

Note that the discontinuity in the diagonal pairwise plots is indicated using a dashed line (as in the one-dimensional setting). Similarly, in the off-diagonal plots the point of discontinuity (from \(W_{t-}\) to \(W_t\)) is joined using a dashed line.

Path concatenation is also possible, but care has to be taken to concatenate paths of both the same dimension and covariance structure (there are a number of checks carried out in this scenario to ensure this is the case). For instance, consider concatenating bm.d with itself into a new object named concat.bm.d, initialised at \(t = 5\):

concat.bm.d <- concat.bm(bm.d, bm.d, t0 = 5)
plot(concat.bm.d)

Of course, the off-diagonal plots don’t appear to have changed due to overplotting.

A final peek under the hood

There were some examples above of indexing into the Brownian Motion object which show that the interface in multiple dimensions is in harmony with that of univariate paths. Therefore, all the standard indexing behavior described in the earlier section applies, with the obvious changes to multi-dimensional outputs for W_t and W_tm.

Repeating the earlier examples on a multi-dimensional path:

bm.d[,"t"]
#> [1]  0.0  0.5  5.0 10.0
bm.d[,"W_t"]
#>           [,1]       [,2]     [,3]
#> [1,]  0.000000  0.0000000  0.00000
#> [2,]  0.775882 -0.7164934  0.16457
#> [3,] 12.000000 10.0000000 10.00000
#> [4,] 12.740161 13.0331808  9.02344
bm.d[t = 0.5, "W_t"]
#>          [,1]       [,2]    [,3]
#> [1,] 0.775882 -0.7164934 0.16457
bm.d[4, "W_t"]
#>          [,1]     [,2]    [,3]
#> [1,] 12.74016 13.03318 9.02344

The simulation of a path with non-identity covariance matrices relies upon transforming the path to one with identity covariance matrix, simulating in this transformed space, and then doing the reverse transformation as required. In particular, we consider the transformation \(Z_t = W_t (\Lambda^{-1/2})'\), where \(W_t\) is our original path and \(\Lambda\) is our specified covariance matrix, cov.

It is possible to access the transformed (standard) Brownian motion path by accessing the $Z member and using this in the usual way (ie to access times, positions etc), noting that the position is now denoted Z_t rather than W_t in this space:

bm.d$Z[,"t"]
#> [1]  0.0  0.5  5.0 10.0
bm.d$Z[,"Z_t"]
#>          [,1]       [,2]       [,3]
#> [1,] 0.000000  0.0000000  0.0000000
#> [2,] 1.122364 -0.8945207  0.2015562
#> [3,] 2.144014  7.4645225 12.2474487
#> [4,] 1.531647 11.3656170 11.0514117
bm.d$Z[t = 0.5, "Z_t"]
#>          [,1]       [,2]      [,3]
#> [1,] 1.122364 -0.8945207 0.2015562
bm.d$Z[4, "Z_t"]
#>          [,1]     [,2]     [,3]
#> [1,] 1.531647 11.36562 11.05141

Layer information can be obtained in the standard way:

bm.d$Z[,"layers"]
#> # A tibble: 3 × 5
#>     t.l   t.u outer.L[,1]  [,2]   [,3] outer.U[,1]   [,2]   [,3] inner.cube   
#>   <dbl> <dbl>       <dbl> <dbl>  <dbl>       <dbl>  <dbl>  <dbl> <list>       
#> 1   0     0.5     -0.0884 -1.16 -0.530        1.21  0.265  0.732 <dbl [8 × 3]>
#> 2   0.5   5       -2.48   -6.68 -2.43         1.92  0.166  0.732 <dbl [8 × 3]>
#> 3   5    10       -0.704   5.23  8.82         4.38 13.6   14.5   <dbl [8 × 3]>

Again, each row of the tibble specifies layers information for an interval of time (between t.l and t.u). Note however that the formatting of the layer information in each row, in multiple dimensions, is slightly different. Here, for ease of representation, we only supply the outer most layer for each component (Ld and Uu). This is then composed in vectors of length dim in outer.L (the lower bounds) and outer.U (the upper bounds) respectively. For instance, element 2 of outer.L is Ld in component 2 of \(Z\). inner.cube is a matrix of size \(2^{\text{dim}} \times \text{dim}\) specifying the vertices of the hypercube obtained by composing for each component of \(Z\) the vectors outer.L and outer.U.

To access the full detail (including full layer information) for each component of \(Z\) (for instance, component \(2\)), this can be achieved as follows:

bm.d$Z.bm[[2]][]
#> $t
#> [1]  0.0  0.5  5.0 10.0
#> 
#> $W_t
#> [1]  0.0000000 -0.8945207  7.4645225 11.3656170
#> 
#> $W_tm
#> [1]  0.0000000 -0.8945207 -5.6146243 11.3656170
#> 
#> $layers
#> # A tibble: 3 × 9
#>   type     t.l   t.u    Ld     Uu    Lu      Ud Lu.hard Ud.hard
#>   <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl>   <dbl> <lgl>   <lgl>  
#> 1 bessel   0     0.5 -1.16  0.265 -1.12  0.221  FALSE   FALSE  
#> 2 bessel   0.5   5   -6.68  0.166 -6.54  0.0336 FALSE   FALSE  
#> 3 bessel   5    10    5.23 13.6    7.46 11.4    TRUE    TRUE   
#> 
#> $labels
#> $labels$start
#> [1] 0
#> 
#> $labels$seg.start
#> [1] 0 5
#> 
#> $labels$forced
#> [1] 0
#> 
#> $labels$user
#> [1]  0.5  5.0 10.0
#> 
#> $labels$end
#> [1] 10
#> 
#> $labels$seg.end
#> [1]  5 10

Note that this is in exactly the same format as the discussion in the one-dimensional setting (indeed, it is a one-dimensional Brownian motion).

The full \(Z\) object can be plotted as follows:

plot(bm.d$Z)

Notice that the diagonal plots are now demarked \(Z\) (as opposed to \(W\)), and again are coloured (following the same convention as the one dimensional plotting). The off-diagonal plots have the same interpretation as previously discussed in the \(W\) setting.

Now, returning to the original (untransformed) path, \(W\). In addition to using indexing to obtain t, W_t, and labels, we can also obtain layer information in the usual way:

bm.d[,"layers"]
#> # A tibble: 3 × 5
#>     t.l   t.u outer.L[,1]  [,2]   [,3] outer.U[,1]   [,2]   [,3] inner.cube   
#>   <dbl> <dbl>       <dbl> <dbl>  <dbl>       <dbl>  <dbl>  <dbl> <list>       
#> 1   0     0.5      -0.933 -1.16 -0.433        1.71  0.441  0.598 <dbl [8 × 3]>
#> 2   0.5   5        -7.04  -6.48 -1.99         2.37  0.355  0.598 <dbl [8 × 3]>
#> 3   5    10         6.32   7.07  7.20        18.4  16.0   11.8   <dbl [8 × 3]>

Here, the layer convention follows that discussed for \(Z\) above. The inner.cube is obtained by doing the reverse transformation of each of the \(2^{\text{dim}} \times \text{dim}\) vertices of the inner.cube in the transformed \(Z\) object. As such, it specifies a region of \(W\) space which contains the entire path over an interval of time. outer.L (resp. outer.U) is then composed by taking the minimum (resp. maximum) over each coordinate of the vertices. As such, outer.L and outer.U specify an outer bounding hyper-cube, parallel to the axis in the \(W\) space.

Plotting

We return to plotting now we more fully understand the layers of the Brownian motion object. Consider the \(W\) plot:

plot(bm.d)

Here, in the off-diagonal plots to construct the region demarcated by the red lines, we have simply taken the convex hull of all vertices given by the inner cube over the specified time window, and projected the hyper cube onto the relevant two dimensional components. The interpretation is that the entire path is constrained in the inner red region for the entire time requested. The region given by the black lines is simply an outer hypercube parallel to all axis which bounds the red region (again, composed by taking the minimums and maximums over every component), which in practice will be somewhat easier to use.

Of course, we may wish to only plot over a shorter period of time, for instance \(t=[0,0.5]\), which can be achieved using the t.lim argument as before:

plot(bm.d, t.lim = c(0,0.5))

Notice, if we plot for \(Z\) the same time interval, then as the time interval only comprises a single layer, then both red and black demarked regions overlap, and they are parallel to the axis:

plot(bm.d$Z, t.lim = c(0,0.5))

Of course, we may wish a more refined look at the layer information. This can be achieved using the argument layers.2d and specifying "convex" (the default, described above), "detailed", or "union". "detailed" simply plots each layer seperately. For instance (returning to the entire interval \(t=[0,10]\)):

plot(bm.d, layers.2d = "detailed")

Note here red and black demarked regions are drawn for every layer in the requested time interval. Visually it is clear in this instance that there are two regions in which the path is constrained – corresponding to the interval prior to the addition of a path segment, and after.

In cases of paths with lots of layer information, "detailed" can become cluttered. Instead, here we can use "union" which takes the union of all red demarked regions. Note that (i) this depends on an external package, sf, which is very large; and (ii) is more computationally expensive to compute. Hence, although it is a more accurate portrayal of the layer information, it is not the default option (and sf is only a “Suggested” not required package).

plot(bm.d, layers.2d = "union")

Here, we have identified the two red demarked regions corresponding to the union of two dimensional projections of the hypercubes in every interval prior to the discontinuity, and every interval after. As before, the black demarked region identifies a bounding hypercube parallel to the axis.

The remaining plotting options are:

  • pairs (TRUE/FALSE) to indicate whether to produce the full pairwise plots (TRUE), or only off-diagonal plots (FALSE).
  • dims which enables listing only a subset of dimensions to plot.

Together, pairs and dims allow You to plot any part of the pairwise plots. For instance, if we wanted to plot only the one-dimensional component 2, then this achieved by:

plot(bm.d, dims = 2)

Whereas, if we wanted a "union" plot of component 3 against component 2, then this is possible as follows:

plot(bm.d, layers.2d = "union", dims = c(3,2), pairs = FALSE)

References

Beskos, A., O. Papaspiliopoulos, and G. O. Roberts. 2008. “A factorisation of diffusion measure and finite sample path constructions.” Methodology and Computing in Applied Probability 10: 85–104.
Devroye, L. 2009. “On exact simulation algorithms for some distributions related to Jacobi theta functions.” Statistics and Probability Letters 79: 2251–59.
Pollock, M., P. Fearnhead, A. M. Johansen, and G. O. Roberts. 2020. “Quasi-stationary Monte Carlo methods and the ScaLE algorithm (with discussion).” Journal of the Royal Statistical Society, Series B (Statistical Methodology) 82: 1–59.
Pollock, M., A. M. Johansen, and G. O. Roberts. 2016. “On the exact and \(\varepsilon\)-strong simulation of (jump) diffusions.” Bernoulli 22 (2): 794–856.