In rchan26/layeredBB: Simulates layered Brownian bridges using C++ and Rcpp
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(layeredBB)
Simulating Brownian Motion
To simulate Brownian motion sample paths, we use the Brownian_motion_path_sampler function. This function takes in the start point of the Brownian motion as the first argument and a vector of times to simulate the path as the second argument. For instance, we can simulate a standard Brownian motion path starting at $x=0$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ using the following code:
set.seed(2021)
bm_paths <- lapply(1:3, function(i) Brownian_motion_path_sampler(0, seq(0, 1, 0.01)))
# plotting the simulated paths
plot(x = bm_paths[[1]]['time',], y = bm_paths[[1]]['X',], type = 'l', ylim = c(-2, 2),
xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n')
lines(x = bm_paths[[2]]['time',], y = bm_paths[[2]]['X',], lwd = 3, lty = 2)
lines(x = bm_paths[[3]]['time',], y = bm_paths[[3]]['X',], lwd = 3, lty = 3)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5)
axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5)
We can also verify that our code is correct by simulating many Brownian motion paths and checking the distribution of the paths at a certain time. For instance, if we simulate a Brownian motion path starting at $x$ from time $s$ to $t$, then by using the properties of Brownian motion, we know that the distribution of the path $X$ at time $t$ should be $\mathcal{N}(x,|t-s|)$.
First setting $x=0$, $s=0$ and $t=2$, we can simulate many paths (we choose to simulate $10000$ paths here):
x <- 0
s <- 0
t <- 2
replicates <- 10000
test <- lapply(1:replicates, function(i) Brownian_motion_path_sampler(x, seq(s, t, 0.01)))
First plotting the simulated sample Brownian motion paths:
plot(x = test[[1]]['time',], y = test[[1]]['X',], type = 'l',
ylim = c(-4*sqrt(t-s), 4*sqrt(t-s)),
xlab = '', ylab = '', lwd = 0.5, xaxt = 'n', yaxt = 'n')
for (i in 2:replicates) {
lines(x = test[[i]]['time',], y = test[[i]]['X',], lwd = 0.5)
}
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(s, t, 0.2), labels=seq(s, t, 0.2), font = 2, cex = 1.5)
axis(1, at=seq(s, t, 0.1), labels=rep("", length(seq(s, t, 0.1))), lwd.ticks = 0.5)
axis(2, at=seq(floor(min(sapply(test, function(path) path['X',]))),
ceiling(max(sapply(test, function(path) path['X',]))), 1),
labels=seq(floor(min(sapply(test, function(path) path['X',]))),
ceiling(max(sapply(test, function(path) path['X',]))), 1),
font = 2, cex = 1.5)
Then we can simply look at the end points of all these paths at $t=2$ and look at the density of these end points. Theoretically, they should be from a $\mathcal{N}(0,2)$ distribution:
# select the points at the end time t
end_points <- sapply(1:replicates, function(i) test[[i]]['X', which(seq(s, t, 0.01)==t)])
# plot distribution of the simulated points and the theoretical distribution
curve(dnorm(x, 0, sqrt(t-s)), -4*sqrt(t-s), 4*sqrt(t-s), lwd = 2,
ylab = 'density')
lines(density(end_points), col = 'red', lty = 2, lwd = 2)
legend('topleft',
legend = c('theoretical density', 'simulated density'),
col = c('black', 'red'),
lty = c(1, 2),
lwd = c(2, 2))
Multivariate Brownian Motion (dimension = 2)
set.seed(2021)
multi_bm_path <- multi_brownian_motion(dim = 2, x = c(0,0), times = seq(0, 1, 0.01))
plot(x = multi_bm_path[1,], y = multi_bm_path[2,], type = 'l',
xlim = c(-2,2), ylim = c(-2,2), xlab = '', ylab = '', lwd = 3)
mtext(expression(X[1]), 1, 2.75, font = 2, cex = 1.5)
mtext(expression(X[2]), 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5)
axis(1, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5)
axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5)
axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5)
points(x = c(multi_bm_path[1,1], multi_bm_path[1,ncol(multi_bm_path)]),
y = c(multi_bm_path[2,1], multi_bm_path[2,ncol(multi_bm_path)]), pch = 19, cex = 1)
Simualting Brownian Bridges
To simulate Brownian bridge sample paths (which are simply Brownian motion paths conditioned on an end point), we use the Brownian_bridge_path_sampler function. We need to pass in the start and end points into x and y respectively and the start and end times of the bridge at s and t respectively. Lastly, we need to pass in the times for which to simulate the path for into the times argument. For instance, we can simulate a Brownian bridge path starting from $x=0$ at time $s=0$, to $y=0$ and time $t=1$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ using the following code:
set.seed(2021)
bb_paths <- lapply(1:3, function(i) {
Brownian_bridge_path_sampler(0, 0, 0, 1, seq(0, 1, 0.01))$full_path})
# plotting the simulated paths
plot(x = bb_paths[[1]]['time',], y = bb_paths[[1]]['X',], type = 'l', ylim = c(-2, 2),
xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n')
lines(x = bb_paths[[2]]['time',], y = bb_paths[[2]]['X',], lwd = 3, lty = 2)
lines(x = bb_paths[[3]]['time',], y = bb_paths[[3]]['X',], lwd = 3, lty = 3)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5)
axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5)
points(x = c(0,1), y = c(0,0), pch = 19, cex = 1)
As with Brownian motion, we can verify that our code is correct by simulating many Brownian bridge sample paths and checking the distribution of the paths at a certain time. For instance, if we simulate a Brownian bridge path starting from $x$ at time $s$ to $y$ at time $t$, then it can be shown using the properties of Brownian motion, that the distribution of the path at an intermediate time $q \in (s,t)$ is:
$$
\mathcal{N} \left( \frac{(t-q)x + (q-s)y}{(t-s)}, \frac{(t-q)(q-s)}{(t-s)} \right).
$$
x <- 0
y <- 0
s <- 0
t <- 1
q <- 0.5
replicates <- 10000
paths <- lapply(1:replicates, function(i) {
Brownian_bridge_path_sampler(x = x,
y = y,
s = s,
t = t,
times = seq(s, t, 0.01))$full_path})
plot(x = paths[[1]]['time',], y = paths[[1]]['X',], type = 'l',
ylim = c(-4*sqrt(t-s), 4*sqrt(t-s)),
xlab = '', ylab = '', lwd = 0.5, xaxt = 'n', yaxt = 'n')
for (i in 2:replicates) {
lines(x = paths[[i]]['time',], y = paths[[i]]['X',], lwd = 0.5)
}
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(s, t, 0.2), labels=seq(s, t, 0.2), font = 2, cex = 1.5)
axis(1, at=seq(s, t, 0.1), labels=rep("", length(seq(s, t, 0.1))), lwd.ticks = 0.5)
axis(2, at=seq(floor(min(sapply(paths, function(path) path['X',]))),
ceiling(max(sapply(paths, function(path) path['X',]))), 1),
labels=seq(floor(min(sapply(paths, function(path) path['X',]))),
ceiling(max(sapply(paths, function(path) path['X',]))), 1),
font = 2, cex = 1.5)
Then we can simply look at the end points of all these paths at $q=0.5$ and look at the density of these end points. Theoretically, they should be from a $\mathcal{N} \left( \frac{(t-q)x + (q-s)y}{(t-s)}, \frac{(t-q)(q-s)}{(t-s)} \right) = \mathcal{N} \left( 0, 0.25 \right)$ distribution:
# select the points at the specified time q
index <- which(seq(s, t, 0.01)==q)
simulated_points <- sapply(1:replicates, function(i) paths[[i]]['X', index])
# calculate the theoretical mean and standard deviation of the simulated points at time q
theoretical_mean <- x + (q-s)*(y-x)/(t-s)
theoretical_sd <- sqrt((t-q)*(q-s)/(t-s))
# plot distribution of the simulated points and the theoretical distribution
curve(dnorm(x, theoretical_mean, theoretical_sd),
theoretical_mean-4*theoretical_sd,
theoretical_mean+4*theoretical_sd, lwd = 2,
ylab = 'density')
lines(density(simulated_points), col = 'red', lty = 2, lwd = 2)
legend('topleft',
legend = c('theoretical density', 'simulated density'),
col = c('black', 'red'),
lty = c(1, 2),
lwd = c(2, 2))
Simulating minimum and maximum values of a Brownian bridge
We can also simulate minimum and maximum values of a Brownian bridge using the min_sampler (for minima simulation) and max_sampler (for maxima simulation) functions. For these functions, we must provide a lower and upper bound for where the minimum or maximum could occur. For instance, if we set our bounds to be very large (in essence almost an unconstrained case), we can use the following code:
set.seed(2021)
replicates <- 10000
x <- 0
y <- 0
s <- 0
t <- 1
min_low_bound <- min(x,y)-1000
min_up_bound <- min(x,y)
max_low_bound <- max(x,y)
max_up_bound <- max(x,y)+1000
sim_min <- sapply(1:replicates, function(i) min_sampler(x = x,
y = y,
s = s,
t = t,
low_bound = min_low_bound,
up_bound = min_up_bound))
sim_max <- sapply(1:replicates, function(i) max_sampler(x = x,
y = y,
s = s,
t = t,
low_bound = max_low_bound,
up_bound = max_up_bound))
# plotting minima and maxima
plot(x = sim_min[2,], y = sim_min[1,], ylim = c(-2.5, 2.5),
pch = 20, lwd = 0.5, cex = 0.25,
xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
points(x = sim_max[2,], y = sim_max[1,], pch = 20, lwd = 0.5, cex = 0.25)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5)
axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5)
If we want to be more restrictive of where the minimum or maximum occurs, we can simply pass in more restrictive bounds for where they can occur through the low_bound and up_bound arguments.
set.seed(2021)
replicates <- 10000
x <- 0
y <- 0
s <- 0
t <- 1
min_low_bound <- min(x,y)-1.5
min_up_bound <- min(x,y)-0.5
max_low_bound <- max(x,y)+1
max_up_bound <- max(x,y)+1.5
sim_min_rest <- sapply(1:replicates, function(i) min_sampler(x = x,
y = y,
s = s,
t = t,
low_bound = min_low_bound,
up_bound = min_up_bound))
sim_max_rest <- sapply(1:replicates, function(i) max_sampler(x = x,
y = y,
s = s,
t = t,
low_bound = max_low_bound,
up_bound = max_up_bound))
# plotting minima and maxima
plot(x = sim_min_rest[2,], y = sim_min_rest[1,], ylim = c(-2.5, 2.5),
pch = 20, lwd = 0.5, cex = 0.25,
xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
points(x = sim_max_rest[2,], y = sim_max_rest[1,], pch = 20, lwd = 0.5, cex = 0.25)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5)
axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5)
abline(h=c(min_low_bound, min_up_bound, max_low_bound, max_up_bound), lwd = 2)
Simulating a Bessel bridge
A Bessel bridge is a Brownian bridge conditional on the minimum or maximum of the path. We can simulate Bessel bridges using min_Bessel_bridge_path_sampler (for Bessel bridges conditional on the minimum) and max_Bessel_bridge_path_sampler (for Bessel bridges conditional on the maximum). For each method, we must pass in the usual start and end points of the Brownian bridge with x, y, s and t denoting the start point, end point, start time, end time of the bridge respectively. We then pass the minimum (or maximum) at the time which it occurs in m and tau. Lastly, we need to pass in the times for which to simulate the path for into the times argument. For instance, we can simulate a Bessel bridge path starting from $x=0$ at time $s=0$, to $y=0$ and time $t=1$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ conditional on having a minimum $m=-2$ and time $\tau=0.5$ using the following code:
set.seed(2021)
x <- 0
y <- 0
s <- 0
t <- 1
min <- c('min' = -2, 'tau' = 0.5)
besb_paths <- lapply(1:3, function(i) {
return(min_Bessel_bridge_path_sampler(x = x,
y = y,
s = s,
t = t,
m = min['min'],
tau = min['tau'],
times = seq(s, t, 0.01))$full_path)
})
# plotting simulated paths
plot(x = besb_paths[[1]]['time',], y = besb_paths[[1]]['X',], type = 'l', ylim = c(-3, 3),
xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n')
lines(x = besb_paths[[2]]['time',], y = besb_paths[[2]]['X',], lwd = 3, lty = 2)
lines(x = besb_paths[[3]]['time',], y = besb_paths[[3]]['X',], lwd = 3, lty = 3)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5)
axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5)
points(x = c(s,t), y = c(x,y), pch = 19, cex = 1)
points(x = min['tau'], y = min['min'], pch = 19, cex = 1.5)
If instead we condition the bridge on achieving a maximum value of $m=2.5$ and time $\tau=0.75$, we can use:
set.seed(2021)
x <- 0
y <- 0
s <- 0
t <- 1
max <- c('max' = 2.5, 'tau' = 0.75)
besb_paths <- lapply(1:3, function(i) {
return(max_Bessel_bridge_path_sampler(x = x,
y = y,
s = s,
t = t,
m = max['max'],
tau = max['tau'],
times = seq(s, t, 0.01))$full_path)
})
# plotting simulated paths
plot(x = besb_paths[[1]]['time',], y = besb_paths[[1]]['X',], type = 'l', ylim = c(-3, 3),
xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n')
lines(x = besb_paths[[2]]['time',], y = besb_paths[[2]]['X',], lwd = 3, lty = 2)
lines(x = besb_paths[[3]]['time',], y = besb_paths[[3]]['X',], lwd = 3, lty = 3)
mtext('Time', 1, 2.75, font = 2, cex = 1.5)
mtext('W', 2, 2.75, font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5)
axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5)
axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5)
axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5)
points(x = c(s,t), y = c(x,y), pch = 19, cex = 1)
points(x = max['tau'], y = max['max'], pch = 19, cex = 1.5)
rchan26/layeredBB documentation built on March 25, 2022, 3:44 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(layeredBB)
Simulating Brownian Motion
To simulate Brownian motion sample paths, we use the Brownian_motion_path_sampler function. This function takes in the start point of the Brownian motion as the first argument and a vector of times to simulate the path as the second argument. For instance, we can simulate a standard Brownian motion path starting at $x=0$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ using the following code:
set.seed(2021) bm_paths <- lapply(1:3, function(i) Brownian_motion_path_sampler(0, seq(0, 1, 0.01))) # plotting the simulated paths plot(x = bm_paths[[1]]['time',], y = bm_paths[[1]]['X',], type = 'l', ylim = c(-2, 2), xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n') lines(x = bm_paths[[2]]['time',], y = bm_paths[[2]]['X',], lwd = 3, lty = 2) lines(x = bm_paths[[3]]['time',], y = bm_paths[[3]]['X',], lwd = 3, lty = 3) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5) axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5)
We can also verify that our code is correct by simulating many Brownian motion paths and checking the distribution of the paths at a certain time. For instance, if we simulate a Brownian motion path starting at $x$ from time $s$ to $t$, then by using the properties of Brownian motion, we know that the distribution of the path $X$ at time $t$ should be $\mathcal{N}(x,|t-s|)$.
First setting $x=0$, $s=0$ and $t=2$, we can simulate many paths (we choose to simulate $10000$ paths here):
x <- 0 s <- 0 t <- 2 replicates <- 10000 test <- lapply(1:replicates, function(i) Brownian_motion_path_sampler(x, seq(s, t, 0.01)))
First plotting the simulated sample Brownian motion paths:
plot(x = test[[1]]['time',], y = test[[1]]['X',], type = 'l', ylim = c(-4*sqrt(t-s), 4*sqrt(t-s)), xlab = '', ylab = '', lwd = 0.5, xaxt = 'n', yaxt = 'n') for (i in 2:replicates) { lines(x = test[[i]]['time',], y = test[[i]]['X',], lwd = 0.5) } mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(s, t, 0.2), labels=seq(s, t, 0.2), font = 2, cex = 1.5) axis(1, at=seq(s, t, 0.1), labels=rep("", length(seq(s, t, 0.1))), lwd.ticks = 0.5) axis(2, at=seq(floor(min(sapply(test, function(path) path['X',]))), ceiling(max(sapply(test, function(path) path['X',]))), 1), labels=seq(floor(min(sapply(test, function(path) path['X',]))), ceiling(max(sapply(test, function(path) path['X',]))), 1), font = 2, cex = 1.5)
Then we can simply look at the end points of all these paths at $t=2$ and look at the density of these end points. Theoretically, they should be from a $\mathcal{N}(0,2)$ distribution:
# select the points at the end time t end_points <- sapply(1:replicates, function(i) test[[i]]['X', which(seq(s, t, 0.01)==t)]) # plot distribution of the simulated points and the theoretical distribution curve(dnorm(x, 0, sqrt(t-s)), -4*sqrt(t-s), 4*sqrt(t-s), lwd = 2, ylab = 'density') lines(density(end_points), col = 'red', lty = 2, lwd = 2) legend('topleft', legend = c('theoretical density', 'simulated density'), col = c('black', 'red'), lty = c(1, 2), lwd = c(2, 2))
Multivariate Brownian Motion (dimension = 2)
set.seed(2021) multi_bm_path <- multi_brownian_motion(dim = 2, x = c(0,0), times = seq(0, 1, 0.01)) plot(x = multi_bm_path[1,], y = multi_bm_path[2,], type = 'l', xlim = c(-2,2), ylim = c(-2,2), xlab = '', ylab = '', lwd = 3) mtext(expression(X[1]), 1, 2.75, font = 2, cex = 1.5) mtext(expression(X[2]), 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5) axis(1, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5) axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5) axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5) points(x = c(multi_bm_path[1,1], multi_bm_path[1,ncol(multi_bm_path)]), y = c(multi_bm_path[2,1], multi_bm_path[2,ncol(multi_bm_path)]), pch = 19, cex = 1)
Simualting Brownian Bridges
To simulate Brownian bridge sample paths (which are simply Brownian motion paths conditioned on an end point), we use the Brownian_bridge_path_sampler function. We need to pass in the start and end points into x and y respectively and the start and end times of the bridge at s and t respectively. Lastly, we need to pass in the times for which to simulate the path for into the times argument. For instance, we can simulate a Brownian bridge path starting from $x=0$ at time $s=0$, to $y=0$ and time $t=1$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ using the following code:
set.seed(2021) bb_paths <- lapply(1:3, function(i) { Brownian_bridge_path_sampler(0, 0, 0, 1, seq(0, 1, 0.01))$full_path}) # plotting the simulated paths plot(x = bb_paths[[1]]['time',], y = bb_paths[[1]]['X',], type = 'l', ylim = c(-2, 2), xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n') lines(x = bb_paths[[2]]['time',], y = bb_paths[[2]]['X',], lwd = 3, lty = 2) lines(x = bb_paths[[3]]['time',], y = bb_paths[[3]]['X',], lwd = 3, lty = 3) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-2, 2, 1), labels=seq(-2, 2, 1), font = 2, cex = 1.5) axis(2, at=seq(-2, 2, 0.5), labels=rep("", 9), lwd.ticks = 0.5) points(x = c(0,1), y = c(0,0), pch = 19, cex = 1)
As with Brownian motion, we can verify that our code is correct by simulating many Brownian bridge sample paths and checking the distribution of the paths at a certain time. For instance, if we simulate a Brownian bridge path starting from $x$ at time $s$ to $y$ at time $t$, then it can be shown using the properties of Brownian motion, that the distribution of the path at an intermediate time $q \in (s,t)$ is: $$ \mathcal{N} \left( \frac{(t-q)x + (q-s)y}{(t-s)}, \frac{(t-q)(q-s)}{(t-s)} \right). $$
x <- 0 y <- 0 s <- 0 t <- 1 q <- 0.5 replicates <- 10000 paths <- lapply(1:replicates, function(i) { Brownian_bridge_path_sampler(x = x, y = y, s = s, t = t, times = seq(s, t, 0.01))$full_path})
plot(x = paths[[1]]['time',], y = paths[[1]]['X',], type = 'l', ylim = c(-4*sqrt(t-s), 4*sqrt(t-s)), xlab = '', ylab = '', lwd = 0.5, xaxt = 'n', yaxt = 'n') for (i in 2:replicates) { lines(x = paths[[i]]['time',], y = paths[[i]]['X',], lwd = 0.5) } mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(s, t, 0.2), labels=seq(s, t, 0.2), font = 2, cex = 1.5) axis(1, at=seq(s, t, 0.1), labels=rep("", length(seq(s, t, 0.1))), lwd.ticks = 0.5) axis(2, at=seq(floor(min(sapply(paths, function(path) path['X',]))), ceiling(max(sapply(paths, function(path) path['X',]))), 1), labels=seq(floor(min(sapply(paths, function(path) path['X',]))), ceiling(max(sapply(paths, function(path) path['X',]))), 1), font = 2, cex = 1.5)
Then we can simply look at the end points of all these paths at $q=0.5$ and look at the density of these end points. Theoretically, they should be from a $\mathcal{N} \left( \frac{(t-q)x + (q-s)y}{(t-s)}, \frac{(t-q)(q-s)}{(t-s)} \right) = \mathcal{N} \left( 0, 0.25 \right)$ distribution:
# select the points at the specified time q index <- which(seq(s, t, 0.01)==q) simulated_points <- sapply(1:replicates, function(i) paths[[i]]['X', index]) # calculate the theoretical mean and standard deviation of the simulated points at time q theoretical_mean <- x + (q-s)*(y-x)/(t-s) theoretical_sd <- sqrt((t-q)*(q-s)/(t-s)) # plot distribution of the simulated points and the theoretical distribution curve(dnorm(x, theoretical_mean, theoretical_sd), theoretical_mean-4*theoretical_sd, theoretical_mean+4*theoretical_sd, lwd = 2, ylab = 'density') lines(density(simulated_points), col = 'red', lty = 2, lwd = 2) legend('topleft', legend = c('theoretical density', 'simulated density'), col = c('black', 'red'), lty = c(1, 2), lwd = c(2, 2))
Simulating minimum and maximum values of a Brownian bridge
We can also simulate minimum and maximum values of a Brownian bridge using the min_sampler (for minima simulation) and max_sampler (for maxima simulation) functions. For these functions, we must provide a lower and upper bound for where the minimum or maximum could occur. For instance, if we set our bounds to be very large (in essence almost an unconstrained case), we can use the following code:
set.seed(2021) replicates <- 10000 x <- 0 y <- 0 s <- 0 t <- 1 min_low_bound <- min(x,y)-1000 min_up_bound <- min(x,y) max_low_bound <- max(x,y) max_up_bound <- max(x,y)+1000 sim_min <- sapply(1:replicates, function(i) min_sampler(x = x, y = y, s = s, t = t, low_bound = min_low_bound, up_bound = min_up_bound)) sim_max <- sapply(1:replicates, function(i) max_sampler(x = x, y = y, s = s, t = t, low_bound = max_low_bound, up_bound = max_up_bound)) # plotting minima and maxima plot(x = sim_min[2,], y = sim_min[1,], ylim = c(-2.5, 2.5), pch = 20, lwd = 0.5, cex = 0.25, xlab = '', ylab = '', xaxt = 'n', yaxt = 'n') points(x = sim_max[2,], y = sim_max[1,], pch = 20, lwd = 0.5, cex = 0.25) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5) axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5)
If we want to be more restrictive of where the minimum or maximum occurs, we can simply pass in more restrictive bounds for where they can occur through the low_bound and up_bound arguments.
set.seed(2021) replicates <- 10000 x <- 0 y <- 0 s <- 0 t <- 1 min_low_bound <- min(x,y)-1.5 min_up_bound <- min(x,y)-0.5 max_low_bound <- max(x,y)+1 max_up_bound <- max(x,y)+1.5 sim_min_rest <- sapply(1:replicates, function(i) min_sampler(x = x, y = y, s = s, t = t, low_bound = min_low_bound, up_bound = min_up_bound)) sim_max_rest <- sapply(1:replicates, function(i) max_sampler(x = x, y = y, s = s, t = t, low_bound = max_low_bound, up_bound = max_up_bound)) # plotting minima and maxima plot(x = sim_min_rest[2,], y = sim_min_rest[1,], ylim = c(-2.5, 2.5), pch = 20, lwd = 0.5, cex = 0.25, xlab = '', ylab = '', xaxt = 'n', yaxt = 'n') points(x = sim_max_rest[2,], y = sim_max_rest[1,], pch = 20, lwd = 0.5, cex = 0.25) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5) axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5) abline(h=c(min_low_bound, min_up_bound, max_low_bound, max_up_bound), lwd = 2)
Simulating a Bessel bridge
A Bessel bridge is a Brownian bridge conditional on the minimum or maximum of the path. We can simulate Bessel bridges using min_Bessel_bridge_path_sampler (for Bessel bridges conditional on the minimum) and max_Bessel_bridge_path_sampler (for Bessel bridges conditional on the maximum). For each method, we must pass in the usual start and end points of the Brownian bridge with x, y, s and t denoting the start point, end point, start time, end time of the bridge respectively. We then pass the minimum (or maximum) at the time which it occurs in m and tau. Lastly, we need to pass in the times for which to simulate the path for into the times argument. For instance, we can simulate a Bessel bridge path starting from $x=0$ at time $s=0$, to $y=0$ and time $t=1$ at times $t={0,0.01,0.02,\dots,0.98,0.99,1}$ conditional on having a minimum $m=-2$ and time $\tau=0.5$ using the following code:
set.seed(2021) x <- 0 y <- 0 s <- 0 t <- 1 min <- c('min' = -2, 'tau' = 0.5) besb_paths <- lapply(1:3, function(i) { return(min_Bessel_bridge_path_sampler(x = x, y = y, s = s, t = t, m = min['min'], tau = min['tau'], times = seq(s, t, 0.01))$full_path) }) # plotting simulated paths plot(x = besb_paths[[1]]['time',], y = besb_paths[[1]]['X',], type = 'l', ylim = c(-3, 3), xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n') lines(x = besb_paths[[2]]['time',], y = besb_paths[[2]]['X',], lwd = 3, lty = 2) lines(x = besb_paths[[3]]['time',], y = besb_paths[[3]]['X',], lwd = 3, lty = 3) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5) axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5) points(x = c(s,t), y = c(x,y), pch = 19, cex = 1) points(x = min['tau'], y = min['min'], pch = 19, cex = 1.5)
If instead we condition the bridge on achieving a maximum value of $m=2.5$ and time $\tau=0.75$, we can use:
set.seed(2021) x <- 0 y <- 0 s <- 0 t <- 1 max <- c('max' = 2.5, 'tau' = 0.75) besb_paths <- lapply(1:3, function(i) { return(max_Bessel_bridge_path_sampler(x = x, y = y, s = s, t = t, m = max['max'], tau = max['tau'], times = seq(s, t, 0.01))$full_path) }) # plotting simulated paths plot(x = besb_paths[[1]]['time',], y = besb_paths[[1]]['X',], type = 'l', ylim = c(-3, 3), xlab = '', ylab = '', lwd = 3, xaxt = 'n', yaxt = 'n') lines(x = besb_paths[[2]]['time',], y = besb_paths[[2]]['X',], lwd = 3, lty = 2) lines(x = besb_paths[[3]]['time',], y = besb_paths[[3]]['X',], lwd = 3, lty = 3) mtext('Time', 1, 2.75, font = 2, cex = 1.5) mtext('W', 2, 2.75, font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.2), labels=c("0.0", seq(0.2, 0.8, 0.2), "1.0"), font = 2, cex = 1.5) axis(1, at=seq(0, 1, 0.1), labels=rep("", 11), lwd.ticks = 0.5) axis(2, at=seq(-3, 3, 1), labels=seq(-3, 3, 1), font = 2, cex = 1.5) axis(2, at=seq(-3, 3, 0.5), labels=rep("", 13), lwd.ticks = 0.5) points(x = c(s,t), y = c(x,y), pch = 19, cex = 1) points(x = max['tau'], y = max['max'], pch = 19, cex = 1.5)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.