diff --git a/src/simulations/Simulation.jl b/src/simulations/Simulation.jl
index b030c1c3..4d63f0d0 100644
--- a/src/simulations/Simulation.jl
+++ b/src/simulations/Simulation.jl
@@ -4,6 +4,7 @@
         model::Model
         times::AbstractVector
         X::AbstractArray
+        dZ::Union{AbstractArray, Nothing}
     end
 
 A `Simulation` object represents the result of a Monte Carlo simulation.
@@ -15,11 +16,13 @@ Elements are:
      - `N_1` is `length(m.state_alias)`,
      - `N_2` is number of Monte Carlo paths,
      - `N_3` is `length(times)`.
+  - `dZ` - Brownian motion increments.
 """
 struct Simulation
     model::Model
     times::AbstractVector
     X::AbstractArray
+    dZ::Union{AbstractArray, Nothing}
 end
 
 
@@ -30,7 +33,8 @@ end
         times::AbstractVector,
         n_paths::Int;
         with_progress_bar::Bool = true,
-        brownian_increments::Function = pseudo_brownian_increments
+        brownian_increments::Function = pseudo_brownian_increments,
+        store_brownian_increments::Bool = false,
         )
 
 A simple Monte Carlo simulation method assuming all model components are state-independent.
@@ -41,9 +45,10 @@ function simple_simulation(
     times::AbstractVector,
     n_paths::Int;
     with_progress_bar::Bool = true,
-    brownian_increments::Function = pseudo_brownian_increments
+    brownian_increments::Function = pseudo_brownian_increments,
+    store_brownian_increments::Bool = false,
     )
-    Z = brownian_increments(
+    dZ = brownian_increments(
         length(state_alias(model)),
         n_paths,
         length(times) - 1,
@@ -59,10 +64,13 @@ function simple_simulation(
         (vol, corr) = volatility_and_correlation(model,ch,times[k-1],times[k])
         L = cholesky(corr).L
         # apply diffusion
-        X_t += (sqrt(times[k] - times[k-1]) * (L .* vol)) * Z[:,:,k-1]
+        X_t += (sqrt(times[k] - times[k-1]) * (L .* vol)) * dZ[:,:,k-1]
         X = cat(X, X_t, dims=3)
     end
+    if !store_brownian_increments
+        dZ = nothing
+    end
     #
-    return Simulation(model, times, X)
+    return Simulation(model, times, X, dZ)
 end
 
diff --git a/test/unittests/paths/path.jl b/test/unittests/paths/path.jl
index d2e14c18..65825b27 100644
--- a/test/unittests/paths/path.jl
+++ b/test/unittests/paths/path.jl
@@ -69,7 +69,7 @@ end
     X = ones(length(DiffFusion.state_alias(m)), n_paths, length(times))
     X[:,:,1] = zeros(length(DiffFusion.state_alias(m)), n_paths)
     X[:,:,5] = 2.0 * ones(length(DiffFusion.state_alias(m)), n_paths)
-    sim = DiffFusion.Simulation(m, times, X)
+    sim = DiffFusion.Simulation(m, times, X, nothing)
 
     # valuation context with deterministic dividend yield
 
@@ -293,7 +293,7 @@ end
         @test_throws KeyError DiffFusion.zero_bond(p, t, T, "SXE50:OIS")
         #
         X = zeros(length(DiffFusion.state_alias(m)), n_paths, length(times))
-        sim = DiffFusion.Simulation(m, times, X)  # simplify calculations
+        sim = DiffFusion.Simulation(m, times, X, nothing)  # simplify calculations
         p = DiffFusion.path(sim, ts, context)
         @test isapprox(DiffFusion.asset(p, t, "EUR-USD"), ones(5) * 1.25 * exp(0.01*t), atol=5.0e-15)
         @test isapprox(DiffFusion.asset(p, t, "SXE50"), ones(5) * 3750.00 * exp(0.01*t), atol=5.0e-15)
@@ -368,7 +368,7 @@ end
                 ("SOFR", DiffFusion.FixingEntry("SOFR", "USD-SOFR-Fixings")),
             ]),
         )
-        det_sim = DiffFusion.Simulation(m, zeros(0), zeros(0,1,0) )
+        det_sim = DiffFusion.Simulation(m, zeros(0), zeros(0,1,0), nothing)
         p = DiffFusion.path(det_sim, ts, det_context)
         t = 2.0
         T = 5.0
diff --git a/test/unittests/simulations/simple_models.jl b/test/unittests/simulations/simple_models.jl
index 75e94fbf..6f2e65c4 100644
--- a/test/unittests/simulations/simple_models.jl
+++ b/test/unittests/simulations/simple_models.jl
@@ -154,6 +154,7 @@ using Test
         n_paths = 2^13
         sim = DiffFusion.simple_simulation(m, ch_full, times, n_paths, with_progress_bar = false)
         # println(size(sim.X))
+        @test isnothing(sim.dZ)
         @test size(sim.X) == (9,8192,6)
         # martingale test domestic numeraire
         one = mean(exp.(-sim.X[4,:,:]), dims=1)
@@ -222,6 +223,13 @@ using Test
         @test isapprox(one[4], 0.9967937694265642, atol=abs_tol)
         @test isapprox(one[5], 0.9976008565308596, atol=abs_tol)
         @test isapprox(one[6], 0.9947929666530295, atol=abs_tol)
+        #
+        # store Brownian increments
+        sim_w_dZ = DiffFusion.simple_simulation(m, ch_full, times, n_paths,
+            with_progress_bar = false, store_brownian_increments = true)
+        @test !isnothing(sim_w_dZ.dZ)
+        @test size(sim_w_dZ.dZ) == (9,8192,5)
+        @test sim_w_dZ.X == sim.X
     end
 
 end
\ No newline at end of file