ブラウン橋(Brownian bridge)作りたい

簡単なんだけど、やってなかったなって。ブラウン橋言ってる人もいるし、ブラウニアン橋いってる人もいてどっちが正しいのかわかんねぇけど、ここではブラウン橋でいく。

ブラウン橋は、 $W(t)$ を標準ブラウン運動として
$B_t := W_t - t W_1,\;t \in [0,1]$
と書かれる確率過程 $B(t)$ だと思っておけばいい。こいつの平均値は0で、分散は $t(1-t)$ となる。で、実際にRでやってみましょうと。

このブラウン橋を生成する関数は以下。

bb <- function(size.path, size.timestep)
{
  result <- matrix(0, nrow=size.path, ncol=size.timestep+1)
  dt <- 1/size.timestep
  for(i in 1:size.path)
  {
    z <- cumsum(rnorm(size.timestep))*sqrt(dt)
    result[i,] <- c(0, z - (1:size.timestep)/size.timestep*z[length(z)])
  }
  result
}

実際に走らせてみると、まぁ0.5時点での標準偏差がだいたい0.5になってるので、あってるだろうと。

> size.path <- 10^3
> size.timestep <- 10^2
> x <- bb(size.path, size.timestep)
> sd(x[,size.timestep/2])
[1] 0.5078946
> matplot(t(x[1:100,]), type="l")

描画させたmatplotの結果が以下。

また、このブラウン橋と同じ確率分布となる確率過程は以下のように書ける。
$dY_t := -\frac{Y_t}{1-t}dt + dW_t ,\;t \in [0,1], Y_0=0$
$B_t$ はブラウン運動 $W_t$ の作るσ加法族の $\mathcal{F}_t$ 可測とはならないけど、こちらは $\mathcal{F}_t$ 可測になる点が異なる。

この確率微分方程式をベースに、ブラウン橋を生成する関数は以下。

bb2 <- function(size.path, size.timestep)
{
  result <- matrix(0, nrow=size.path, ncol=size.timestep+1)
  for(i in 1:size.path)
  {
    dt <- 1/size.timestep
    z <- rnorm(size.timestep)*sqrt(dt)
    x <- rep(0, size.timestep+1)
    for(j in 1:size.timestep)
    {
      x[j+1] <- x[j] + -x[j]/(1-(j*dt))*dt + z[j] 
    }
    result[i,] <- x
  }
  result
}

実際に走らせてみると、こちらも0.5時点での標準偏差がだいたい0.5になってるのであってそうだ。

> size.path <- 10^3
> size.timestep <- 10^2
> y <- bb2(size.path, size.timestep)
> sd(y[,size.timestep/2])
[1] 0.5296784
> matplot(t(y[1:100,]), type="l")

描画させたmatplotの結果が以下。

この話をパッケージ使ってやるなら、ウィーン工科大の雑多な関数パッケージである「e1071」パッケージ中にあるrbridge 関数使うとよい。
パスは一本しか出ない。

> library(e1071)
> x <- rbridge()
> plot(x,type="l")

参考

Brownian bridge - Wikipedia