diff --git a/src/models/rates/GaussianHjmModel.jl b/src/models/rates/GaussianHjmModel.jl
index f2ceb4cb..0486f680 100644
--- a/src/models/rates/GaussianHjmModel.jl
+++ b/src/models/rates/GaussianHjmModel.jl
@@ -14,8 +14,18 @@ struct GaussianHjmModelVolatility
     DfT::AbstractMatrix
 end
 
+function volatility(o::GaussianHjmModelVolatility, u::ModelTime)
+    # return o.HHfInv * (o.DfT .* o.sigma_f(u))  # beware DfT multiplication
+    σ = o.sigma_f(u)
+    d = length(σ)
+    return [
+        sum( o.HHfInv[i,k] * o.DfT[k,j] * σ[k] for k = 1:d )
+        for i = 1:d, j = 1:d
+    ]
+end
+
 "Calculate volatility matrix."
-(o::GaussianHjmModelVolatility)(u::ModelTime) = o.HHfInv * (o.DfT .* o.sigma_f(u))  # beware DfT multiplication
+(o::GaussianHjmModelVolatility)(u::ModelTime) = volatility(o, u)
 
 
 """
@@ -282,7 +292,8 @@ function log_zero_bond(m::GaussianHjmModel, model_alias::String, t::ModelTime, T
     G = G_hjm(m, t, T)
     y = func_y(m, t)
     X_ = @view(X.X[idx:idx+(d-1),:])
-    return X_' * G .+ 0.5 * (G'*y*G)
+    GyG = sum(G[i] * sum(y[i,j] * G[j] for j in 1:d) for i in 1:d)
+    return X_' * G .+ (0.5 * GyG)
 end
 
 """
@@ -333,7 +344,9 @@ function log_compounding_factor(
     G2 = G_hjm(m, t, T2)
     y = func_y(m, t)
     X_ = @view(X.X[idx:idx+(d-1),:])
-    return X_' * (G2 .- G1) .+ 0.5 * ((G2'*y*G2) - (G1'*y*G1))
+    G1yG1 = sum(G1[i] * sum(y[i,j] * G1[j] for j in 1:d) for i in 1:d)
+    G2yG2 = sum(G2[i] * sum(y[i,j] * G2[j] for j in 1:d) for i in 1:d)
+    return X_' * (G2 .- G1) .+ 0.5 * (G2yG2 - G1yG1)
 end
 
 
diff --git a/src/models/rates/SeparableHjmModel.jl b/src/models/rates/SeparableHjmModel.jl
index 0cdaa157..435cf5f4 100644
--- a/src/models/rates/SeparableHjmModel.jl
+++ b/src/models/rates/SeparableHjmModel.jl
@@ -36,7 +36,9 @@ end
 Diagonal entries of ``H(s,t)``.
 """
 function H_hjm(chi::AbstractVector, s::ModelTime, t::ModelTime)
-    return exp.(-chi .* (t-s))
+    # exp.(-chi .* (t-s))
+    δt = t-s
+    return [ exp(-χ * δt) for χ in chi ]
 end
 
 """
@@ -54,8 +56,10 @@ end
 Vector function ``G(s,t)``.
 """
 function G_hjm(chi::AbstractVector, s::ModelTime, t::ModelTime)
+    # (1.0 .- exp.(-chi .* (t-s))) ./ chi
+    δt = t-s
     # This is an unsafe implementation. Better use Taylor expansion.
-    return (1.0 .- exp.(-chi .* (t-s))) ./ chi
+    return [ (1.0 - exp(-χ * δt)) / χ for χ in chi ]
 end
 
 """
@@ -93,6 +97,7 @@ end
 Benchmark times volatility scaling matrix ``H [H^f]^{-1} = [H^f H^{-1}]^{-1}``.
 """
 function benchmark_times_scaling(chi::AbstractVector, delta::AbstractVector)
+    # Hf_H_inv = exp.(-delta * chi')
     Hf_H_inv = [ exp(-chi_ * delta_) for delta_ in delta, chi_ in chi ]  # beware the order of loops!
     HHfInv = inv(Hf_H_inv)
     return HHfInv
@@ -118,13 +123,22 @@ function func_y(
     s::ModelTime,
     t::ModelTime,
     )
+    # chi_i_p_chi_j = chi .+ chi'
+    # H_i_j = exp.(-chi_i_p_chi_j .* (t-s))
+    # V = sigmaT * transpose(sigmaT)
+    # G_i_j = (1. .- H_i_j) ./ chi_i_p_chi_j
+    # y1 = y0 .* H_i_j .+ V .* G_i_j
+    #
     # better exploit symmetry and update in-place
-    chi_i_p_chi_j = [ (chi_i + chi_j) for chi_i in chi, chi_j in chi ]
-    H_i_j = exp.(-chi_i_p_chi_j .* (t-s))
-    V = sigmaT * transpose(sigmaT)
     # this is unsafe, better use Taylor expansion
-    G_i_j = (1. .- H_i_j) ./ chi_i_p_chi_j
-    return y0 .* H_i_j .+ V .* G_i_j
+    d = length(chi)
+    δt = t - s
+    return [
+        y0[i,j] * exp(-(chi[i] + chi[j]) * δt) +
+        sum( sigmaT[i,k] * sigmaT[j,k] for k in 1:d ) *
+        (1.0 - exp(-(chi[i] + chi[j]) * δt)) / (chi[i] + chi[j])
+        for i in 1:d, j in 1:d
+    ]
 end
 
 
@@ -146,13 +160,21 @@ function _func_y(
     s::ModelTime,
     t::ModelTime,
     )
+    # chi_i_p_chi_j = chi .+ chi'
+    # H_i_j = exp.(-chi_i_p_chi_j .* (t-s))
+    # V = sigmaT * transpose(sigmaT)
+    # G_i_j = (1. .- H_i_j) ./ chi_i_p_chi_j
+    # y1 = 0 + V .* G_i_j
+    #
     # better exploit symmetry and update in-place
-    chi_i_p_chi_j = [ (chi_i + chi_j) for chi_i in chi, chi_j in chi ]
-    H_i_j = exp.(-chi_i_p_chi_j .* (t-s))
-    V = sigmaT * transpose(sigmaT)
     # this is unsafe, better use Taylor expansion
-    G_i_j = (1. .- H_i_j) ./ chi_i_p_chi_j
-    return V .* G_i_j
+    d = length(chi)
+    δt = t - s
+    return [
+        sum( sigmaT[i,k] * sigmaT[j,k] for k in 1:d ) *
+        (1.0 - exp(-(chi[i] + chi[j]) * δt)) / (chi[i] + chi[j])
+        for i in 1:d, j in 1:d
+    ]
 end
 
 
diff --git a/test/componenttests/scenarios/swap_benchmarks.jl b/test/componenttests/scenarios/swap_benchmarks.jl
new file mode 100644
index 00000000..e61269eb
--- /dev/null
+++ b/test/componenttests/scenarios/swap_benchmarks.jl
@@ -0,0 +1,63 @@
+
+using DiffFusion
+using OrderedCollections
+using Test
+
+
+@testset "Benchmark scenario generation." begin
+
+    @testset "Test Vanilla swaps" begin
+        results = OrderedDict[]
+        for example_string in DiffFusion.Examples.examples
+            @info "Run example " * example_string
+            serialised_example = DiffFusion.Examples.load(example_string)
+            example = DiffFusion.Examples.build(serialised_example)
+            example["config/simulation"]["simulation_times"]["step"] = 0.25
+            example["config/instruments"]["obs_times"]["step"] = 0.25
+            example["config/simulation"]["n_paths"] = 2^13
+            example["config/simulation"]["with_progress_bar"] = false
+            example["config/instruments"]["with_progress_bar"] = false
+            #
+            path_ = DiffFusion.Examples.path!(example)
+            portfolio_ = DiffFusion.Examples.portfolio!(
+                example,
+                32,  # swap
+                0,  # swaptions
+                0,  # berms
+            )
+            legs = vcat(portfolio_...)
+            #
+            config = example["config/instruments"]
+            obs_times = config["obs_times"]
+            if isa(obs_times, AbstractDict)
+                obs_times = Vector(obs_times["start"]:obs_times["step"]:obs_times["stop"])
+            end
+            with_progress_bar = config["with_progress_bar"]
+            discount_curve_key = config["discount_curve_key"]
+            #
+            for leg in legs
+                if isa(leg, DiffFusion.BermudanSwaptionLeg)
+                    DiffFusion.reset_regression!(leg, path_, leg.regression_data.make_regression)
+                end
+            end
+            #
+            GC.gc()
+            time_ = @elapsed @time scens = DiffFusion.scenarios(
+                legs,
+                obs_times,
+                path_,
+                discount_curve_key,
+                with_progress_bar=with_progress_bar
+            )
+            push!(
+                results,
+                OrderedDict(
+                    "example"   => example_string,
+                    "run_time"  => time_,
+                )
+            )    
+        end
+    end
+
+
+end
\ No newline at end of file
diff --git a/test/unittests/models/rates/gaussian_hjm_model.jl b/test/unittests/models/rates/gaussian_hjm_model.jl
index 1038cad9..33fbe27d 100644
--- a/test/unittests/models/rates/gaussian_hjm_model.jl
+++ b/test/unittests/models/rates/gaussian_hjm_model.jl
@@ -98,7 +98,7 @@ using LinearAlgebra
         ]
         @test isapprox(theta_model[1], theta_ref[1], atol=1.5e-13 )
         @test isapprox(theta_model[2], theta_ref[2], atol=7.0e-13 )
-        @test isapprox(theta_model[3], theta_ref[3], atol=3.0e-11 )
+        @test isapprox(theta_model[3], theta_ref[3], atol=4.0e-11 )
         @test isapprox(theta_model[4], theta_ref[4], rtol=7.5e-9  )
         @test isapprox(theta_model[5], theta_ref[5], rtol=2.0e-9  )
     end
@@ -250,8 +250,8 @@ using LinearAlgebra
         SX_2 = DiffFusion.model_state(X, m)
         @test string(SX_2) == string(SX)
         @test DiffFusion.log_bank_account(m, DiffFusion.alias(m), 1.0, SX) == [ 4., 8., 12. ]
-        @test_throws AssertionError DiffFusion.log_bank_account(m, "WrongAlias", 1.0, SX) 
-        @test DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 8.0, SX) == X[1:3,:]' * G .+ 0.5 * G'*y*G
+        @test_throws AssertionError DiffFusion.log_bank_account(m, "WrongAlias", 1.0, SX)
+        @test isapprox(DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 8.0, SX), X[1:3,:]' * G .+ 0.5 * G'*y*G, rtol=1.0e-14)
         @test_throws AssertionError DiffFusion.log_zero_bond(m, "WrongAlias", 4.0, 8.0, SX)
         #
         df1 = DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 8.0, SX)
@@ -262,17 +262,17 @@ using LinearAlgebra
         #
         dfT = DiffFusion.log_zero_bonds(m, DiffFusion.alias(m), 4.0, [4.0, 6.0, 8.0, 10.0], SX)
         @test size(dfT) == (3, 4)
-        @test dfT[:,1] == DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 4.0, SX)
-        @test dfT[:,2] == DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 6.0, SX)
-        @test dfT[:,3] == DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 8.0, SX)
-        @test dfT[:,4] == DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 10.0, SX)
+        @test isapprox(dfT[:,1], DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 4.0, SX), rtol=1.0e-14)
+        @test isapprox(dfT[:,2], DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 6.0, SX), rtol=1.0e-14)
+        @test isapprox(dfT[:,3], DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 8.0, SX), rtol=1.0e-14)
+        @test isapprox(dfT[:,4], DiffFusion.log_zero_bond(m, DiffFusion.alias(m), 4.0, 10.0, SX), rtol=1.0e-14)
         #
         X = [ 0., 0., 1., 2., 3., 4., 0. ] * [ 1., 2., 3.]'
         s_alias = [ "1", "2", "Theta_3F_x_1", "Theta_3F_x_2", "Theta_3F_x_3", "Theta_3F_s", "7" ]
         dict = DiffFusion.alias_dictionary(s_alias)
         SX = DiffFusion.model_state(X, dict)
         @test DiffFusion.log_bank_account(m, "Theta_3F", 1.0, SX) == [ 4., 8., 12. ]
-        @test DiffFusion.log_zero_bond(m, "Theta_3F", 4.0, 8.0, SX) == X[3:5,:]' * G .+ 0.5 * G'*y*G
+        @test isapprox(DiffFusion.log_zero_bond(m, "Theta_3F", 4.0, 8.0, SX), X[3:5,:]' * G .+ 0.5 * G'*y*G, rtol=1.0e-14)
         #
         df1 = DiffFusion.log_zero_bond(m, "Theta_3F", 4.0, 8.0, SX)
         df2 = DiffFusion.log_zero_bond(m, "Theta_3F", 4.0, 10.0, SX)