added antithetic variates for discretized processes

Author: aleCombi

examples/antithetic_mc_heston.jl

@@ -0,0 +1,31 @@
+using DifferentialEquations
+using DiffEqNoiseProcess
+using Plots
+using Random
+using Hedgehog2
+
+
+u0 = [100.0, 0.04]  # Initial (Stock Price, Variance)
+tspan = (0.0, 1.0)  # Simulation for 1 year
+# Create and solve the original problem with saved noise
+rng = MersenneTwister(123)
+prob = Hedgehog2.HestonProblem(0.05, 2.0, 0.04, 0.3, -0.5, u0, tspan)  # Heston parameters
+sol = solve(prob, EM(), dt=1//250, save_noise=true)
+
+# Extract the 2D noise and flip it for antithetic
+W_t = sol.W.t
+W_u = sol.W.W  # This is a Matrix: each row is a Brownian component
+
+# Flip both components for antithetic version
+antithetic_noise = NoiseGrid(W_t, -W_u)
+
+# Create new problem with antithetic noise
+prob_antithetic = remake(prob; noise=antithetic_noise)
+sol_anti = solve(prob_antithetic, EM(), dt=1//250)
+
+# Plot only the first component (stock price)
+plot(sol.t, sol.u .|> x -> x[1], label="Original", linewidth=2)
+plot!(sol_anti.t, sol_anti.u .|> x -> x[1], label="Antithetic", linestyle=:dash, linewidth=2)
+xlabel!("Time")
+ylabel!("Stock Price")
+title!("Heston Process: Original vs Antithetic Paths")

examples/antithetic_montecarlo_EM.jl

@@ -0,0 +1,50 @@
+using DifferentialEquations
+using DiffEqNoiseProcess
+using Random
+using Plots
+
+# Parameters
+μ = 0.05
+σ = 0.2
+t0 = 0.0
+T = 1.0
+S0 = 100.0
+tspan = (t0, T)
+N = 250
+dt = (T - t0) / N
+
+# Define GBM drift and diffusion
+function drift!(du, u, p, t)
+    du[1] = μ * u[1]
+end
+
+function diffusion!(du, u, p, t)
+    du[1] = σ * u[1]
+end
+
+# Create noise with fixed seed
+rng = MersenneTwister(123)
+wiener = WienerProcess(t0, 0.0)
+
+# Define and solve the original GBM SDE
+u0 = [S0]
+sde = SDEProblem(drift!, diffusion!, u0, tspan, noise=wiener)
+sol = DifferentialEquations.solve(sde, dt=dt, save_noise=true)
+sol = DifferentialEquations.solve(sde, dt=dt, save_noise=true)
+plot(sol.W.t, sol.W.u)
+# Extract the used noise path
+# W_vals = sol.W[1, :]  # Extract Brownian path for the first dimension
+
+# Flip it for antithetic path
+wiener_antithetic = NoiseGrid(sol.W.t, -sol.W.W)
+
+# # Define and solve the antithetic GBM SDE
+sde_antithetic = SDEProblem(drift!, diffusion!, u0, tspan, noise=wiener_antithetic)
+sol_antithetic = solve(sde_antithetic, dt=dt)
+
+# # Plot
+plot(sol, label="Original", linewidth=2)
+plot!(sol_antithetic, label="Antithetic", linestyle=:dash, linewidth=2)
+# xlabel!("Time")
+# ylabel!("GBM Value")
+# title!("GBM with Antithetic Variates via Flipped Brownian Path")

src/pricing_methods/montecarlo.jl

@@ -59,26 +59,37 @@ function sde_problem(::LognormalDynamics, s::BlackScholesExact, market::BlackSch
     return NoiseProblem(noise, tspan; s.kwargs...)
 end
 
+function antithetic_sde_problem(d::LognormalDynamics, s::SimulationStrategy, m::BlackScholesInputs, tspan, normal_sol)
+    s_flipped = @set s.seeds = normal_sol.seeds
+    m_flipped = @set m.sigma = -m.sigma
+    return sde_problem(d, s_flipped, m_flipped, tspan)
+end
+
 function sde_problem(::HestonDynamics, ::EulerMaruyama, m::HestonInputs, tspan)
     @assert is_flat(m.rate) "Heston simulation requires flat rate curve"
     rate = zero_rate(m.rate, 0.0)
     return HestonProblem(rate, m.κ, m.θ, m.σ, m.ρ, [m.spot, m.V0], tspan)
 end
 
+function antithetic_sde_problem(d::PriceDynamics, s::EulerMaruyama, m::AbstractMarketInputs, tspan, normal_sol::CustomEnsembleSolution)
+    base_prob = original_sol.solutions[1].prob
+
+    antithetic_modify = function (_base_prob, _seed, i)
+        sol = original_sol.solutions[i]
+        flipped_noise = NoiseGrid(sol.W.t, -sol.W.W)
+        return remake(_base_prob; noise=flipped_noise)
+    end
+
+    return CustomEnsembleProblem(base_prob, normal_sol.seeds, antithetic_modify)
+end
+
 function sde_problem(::HestonDynamics, strategy::HestonBroadieKaya, m::HestonInputs, tspan)
     @assert is_flat(m.rate) "Heston simulation requires flat rate curve"
     rate = zero_rate(m.rate, 0.0)
     noise = HestonNoise(rate, m.κ, m.θ, m.σ, m.ρ, 0.0, [log(m.spot), m.V0], Z0=nothing; strategy.kwargs...)
     return NoiseProblem(noise, tspan)
 end
 
-function antithetic_sde_problem(d::PriceDynamics, s::SimulationStrategy, m::AbstractMarketInputs, tspan, seeds)
-    s_flipped = @set s.seeds = seeds
-    m_flipped = @set m.sigma = -m.sigma
-    return sde_problem(d, s_flipped, m_flipped, tspan)
-end
-
-
 # ------------------ Marginal Laws (optional) ------------------
 
 function marginal_law(::LognormalDynamics, m::BlackScholesInputs, t)
@@ -134,15 +145,14 @@ function simulate_paths(method::MonteCarlo, market_inputs::I, T) where {I <: Abs
     end
 
     # Step 2: simulate antithetic paths using same seeds, flipped sigma
-    seeds = normal_sol.seeds
-    antithetic_prob = antithetic_sde_problem(dynamics, strategy, market_inputs, tspan, seeds)
+    antithetic_prob = antithetic_sde_problem(dynamics, strategy, market_inputs, tspan, normal_sol)
     antithetic_sol = montecarlo_solution(antithetic_prob, strategy)
 
     # Step 3: combine both solutions
     combined_paths = vcat(normal_sol.solutions, antithetic_sol.solutions)
 
     # Create new EnsembleSolution with merged paths
-    return CustomEnsembleSolution(combined_paths, seeds)
+    return CustomEnsembleSolution(combined_paths, normal_sol.seeds)
 end
 
 function solve(prob::PricingProblem{VanillaOption{European, C, Spot}, I}, method::MonteCarlo) where {C, I <: AbstractMarketInputs}