From 9923160f727e35e81e8514f364a57ac263d49321 Mon Sep 17 00:00:00 2001
From: Frank Schaefer <knarf-schaefle@web.de>
Date: Wed, 7 Feb 2024 10:33:37 -0500
Subject: [PATCH] geometric brownian bridge fix

---
 src/geometric_bm.jl | 15 ++++++++++++---
 test/bridge_test.jl | 27 ++++++++++++++++++++++-----
 2 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/src/geometric_bm.jl b/src/geometric_bm.jl
index eea313a..bb3c93b 100644
--- a/src/geometric_bm.jl
+++ b/src/geometric_bm.jl
@@ -27,14 +27,23 @@ https://math.stackexchange.com/questions/412470/conditional-distribution-in-brow
 =#
 function gbm_bridge(dW, gbm, W, W0, Wh, q, h, u, p, t, rng)
     if dW isa AbstractArray
-        return gbm.σ * sqrt((1 - q) * q * abs(h)) * wiener_randn(rng, dW) + q * Wh
+        gbb_mean = @.. log(W0) + q * (log(W0 + Wh) - log(W0))
+        gbb_std = sqrt((1 - q) * q * abs(h)) * gbm.σ
+
+        x = gbb_mean + gbb_std * wiener_randn(rng, dW)
+        return exp.(x) - W0
     else
-        return gbm.σ * sqrt((1 - q) * q * abs(h)) * wiener_randn(rng, typeof(dW)) + q * Wh
+        gbb_mean = log(W0) + q * (log(W0+Wh) - log(W0))
+        gbb_std = sqrt((1 - q) * q * abs(h)) * gbm.σ
+
+        x = gbb_mean + gbb_std * wiener_randn(rng, typeof(dW))
+        return exp(x) - W0
     end
 end
 function gbm_bridge!(rand_vec, gbm, W, W0, Wh, q, h, u, p, t, rng)
     wiener_randn!(rng, rand_vec)
-    @.. rand_vec = gbm.σ * sqrt((1 - q) * q * abs(h)) * rand_vec + q * Wh
+    @.. rand_vec = gbm.σ * sqrt((1 - q) * q * abs(h)) * rand_vec + (log(W0) + q * (log(W0 + Wh) - log(W0)))
+    @.. rand_vec = exp(rand_vec) - W0
 end
 
 @doc doc"""
diff --git a/test/bridge_test.jl b/test/bridge_test.jl
index 6cc1ebf..a9e7d02 100644
--- a/test/bridge_test.jl
+++ b/test/bridge_test.jl
@@ -19,12 +19,29 @@
     @test ≈(timestep_mean(sol, 11)[1], 1.0, atol = 1e-16)
     @test ≈(timestep_meanvar(sol, 11)[2], 0.0, atol = 1e-16)
 
-    μ = 1.2
-    σ = 2.2
-    W = GeometricBrownianBridge(μ, σ, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0)
-    prob = NoiseProblem(W, (0.0, 1.0))
+    μ = 0.10500000000000001
+    σ = 0.1
+    t0 = 0.0
+    tend = 10.0
+    GB0 = 1.0
+    GBend = 4.0
+
+    W = GeometricBrownianBridge(μ, σ, t0, tend, GB0, GBend)
+    prob = NoiseProblem(W, (t0, tend))
     ensemble_prob = EnsembleProblem(prob)
-    @time sol = solve(ensemble_prob, dt = 0.1, trajectories = 100)
+    @time sol = solve(ensemble_prob, dt=1.0, trajectories=100000)
+    ts = t0:1.0:tend
+    for i = 2:10
+        t = ts[i]
+        mean = log(GB0) + (t - t0) / (tend - t0) * (log(GBend) - log(GB0))
+        var = (t - t0) / (tend - t0) * (tend - t) * σ^2
+
+        m, v = timestep_meanvar(sol, i)
+
+        @test ≈(m, exp(mean + var / 2), atol=1e-2)
+        @test ≈(v, (exp(var) - 1) * exp(2 * mean + var), atol=1e-2)
+    end
+
 
     Random.seed!(100)
     r = 100 # should be independent of the rate, so make it crazy