---
title: "Tutorial for BrownianMotion package"
output:
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 4
    toc: true
bibliography: refs.bib
vignette: >
  %\VignetteIndexEntry{Tutorial for BrownianMotion package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

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.

```{r setup}
library(BrownianMotion)
set.seed(212)
```

This tutorial will begin by considering one dimensional standard Brownian motion. In [Multiple Dimensions](#multi) 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$).

```{r}
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:

```{r}
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:

```{r}
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](#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: 

```{r}
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:

```{r}
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:

```{r}
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:

```{r}
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):

```{r}
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 {#peek1}

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$.

```{r}
bm[,"t"]
bm[,"W_t"]
```

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): 

```{r}
bm[t = 1, "W_t"]
```

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$):

```{r}
bm[1, "W_t"]
```

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:

```{r}
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:

```{r}
bm[,"labels"]
```

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](#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)

```{r}
bm[,"labels"]$user
```

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`:

```{r}
bm["user","t"]
```

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:

```{r}
bm |> sim(c(6, 7), label = c("special.time.1", "special.time.2"))
bm[,"labels"]
```

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`:

```{r}
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"]
```

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)$: 

```{r}
bm |>
  delete.labels(label = c("special.time.1")) |>
  delete.labels(t = c(7, 2, 3), label = c("special.time.2"))
bm[,"labels"]
```

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

```{r,error=TRUE,purl=FALSE}
bm |> add.labels(t = 6, label = "end")
```


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`:

```{r}
bm["special.time.2", "W_t"]
```

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$:

```{r}
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!](#layers)). As a result, directly using the `sim()` function (for instance, simulating the path at time $t=-2$) will result in an error:

```{r,error=TRUE,purl=FALSE}
bm |> sim(-2)
```

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:

```{r}
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](#concat))

```{r}
bm.adjust <- concat.bm(bm, t0 = 0)
plot(bm.adjust)
bm.adjust[,"labels"]
```

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,

```{r}
# Current end time of path in bm
bm["end","t"]
# Assign bm to a new variable called bm.x
bm.x <- bm
# Check end time of bm.x
bm.x["end","t"]
# 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"]
bm.x["end","t"]
```

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:

```{r}
# Current end time of path in bm
bm["end","t"]
# Use concat to trigger a copy
bm.x <- concat.bm(bm)
# Check end time of bm.x
bm.x["end","t"]
# 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"]
bm.x["end","t"]
```

## Layers, layers, everywhere! {#layers}

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 [@MCAP].
-   *Localised*: This follows (broadly) the approach outlined in [@DEV], following the implementation in [@SCALE]
-   *Intersection*: This follows (broadly) the approach outlined in [@BERN],

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:

```{r}
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:

```{r}
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]$:

```{r}
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]$:

```{r}
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:

```{r}
bm.adjust |>
  sim(15) |>
  first.passage(delta = 1)
plot(bm.adjust)
```

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:

```{r}
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:

```{r}
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]$:

```{r}
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:

```{r,error=TRUE,purl=FALSE}
bm.adjust |> first.passage(l = 80, u = 85)
```

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).

```{r}
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"`:

```{r}
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,

```{r}
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$):

```{r}
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]$):

```{r}
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:

```{r}
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](#multi).

-   `$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:

```{r,eval=FALSE}
bm.adjust[,"layers"]
```

```{r,echo=FALSE}
rmarkdown::paged_table(bm.adjust[,"layers"])
```

```{r,eval=FALSE}
bm.adjust$user.layers
```

```{r,echo=FALSE,results='asis'}
knitr::kable(bm.adjust$user.layers)
```

```{r,eval=FALSE}
bm.adjust$bb.local[,"layers"]
```

```{r,echo=FALSE,results='asis'}
knitr::kable(bm.adjust$bb.local[,"layers"])
```

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\}$:

```{r}
bm.adjust[t = seq(4.5, 5.5, by = 0.1), "layers"]
```


### 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:

```{r}
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$:

```{r}
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$:

```{r}
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,

```{r}
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"`):

```{r}
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:

```{r}
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):

```{r}
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]$:

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

and then the simulation of intermediary points in this interval:

```{r}
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):

```{r}
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):

```{r}
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:

```{r}
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]$:

```{r}
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:

```{r}
bm3 |> sim(seq(0, 10, by = 0.5))
plot(bm3)
```

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

```{r}
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]$:

```{r}
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$ :

```{r}
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:

```{r}
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:

```{r}
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:

```{r}
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)$:

```{r}
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)$:

```{r}
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:

```{r}
bm5 |> delete.skeleton(l = 5, r = 7.5, type = "label")
bm5[,"labels"]
```


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:

```{r}
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:

```{r}
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:

```{r}
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:

```{r,error=TRUE,purl=FALSE}
bm5 |> delete.skeleton()
```

## 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:

```{r}
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):

```{r}
data.frame(t        = bm6[,"t"],
           W_t      = bm6[,"W_t"],
           W_tm     = bm6[,"W_tm"],
           equality = bm6[,"W_t"] == bm6[,"W_tm"])
```

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:

```{r}
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"])
```

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$:

```{r}
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:

```{r}
bm6[,"labels"]
```

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$:

```{r}
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$:


```{r}
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, 

```{r}
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):

```{r}
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$:

```{r}
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"])
```

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).  

```{r}
bm6[,"labels"]
```

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`:

```{r}
bm6[t = 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$


```{r}
bm.concat <- concat.bm(bm6, bm5, bm6, t0 = 10)
plot(bm.concat)
```

# Multiple Dimensions {#multi}

## 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:

```{r}
bm.d <- create.bm(dim=3)
bm.d$dim
bm.d[t=0, "W_t"]
```

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,

```{r}
bm.d <- create.bm(W_t = 3, dim=3)
bm.d[t=0, "W_t"]
```

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

```{r}
bm.d <- create.bm(W_t = c(0,1,2), dim=3)
bm.d[t=0, "W_t"]
```

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)$

```{r}
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"]
```

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)$.:

```{r}
bm.d |> add.points(t = 0, W_t = c(-1, -1, -1))
bm.d[t=c(0,1,2), "W_t"]
```

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, 

```{r}
bm.d <- create.bm(dim = 3, cov = (diag(x= 0.5, nrow=3) + 0.5))
bm.d$cov
```

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: 

```{r}
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:

```{r}
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](#multidplot)).  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](#multidplot).

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: 

```{r}
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): 

```{r}
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$:

```{r}
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](#peek1) applies, with the obvious changes to multi-dimensional outputs for `W_t` and `W_tm`.

Repeating the earlier examples on a multi-dimensional path:

```{r}
bm.d[,"t"]
bm.d[,"W_t"]
bm.d[t = 0.5, "W_t"]
bm.d[4, "W_t"]
```

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:

```{r}
bm.d$Z[,"t"]
bm.d$Z[,"Z_t"]
bm.d$Z[t = 0.5, "Z_t"]
bm.d$Z[4, "Z_t"]
```

Layer information can be obtained in the standard way: 

```{r}
bm.d$Z[,"layers"]
```

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: 

```{r}
bm.d$Z.bm[[2]][]
```

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:

```{r}
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: 

```{r}
bm.d[,"layers"]
```

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 {#multidplot}

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

```{r}
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: 

```{r}
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: 

```{r}
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]$):

```{r}
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).

```{r}
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: 

```{r}
plot(bm.d, dims = 2)
```

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

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

# References
