Implemented antithetic variates for discretized processes

Author: aleCombi

examples/antithetic_heston_check.jl

@@ -0,0 +1,53 @@
+using Revise
+using Hedgehog2
+using Dates
+using Statistics
+using Random
+using Printf
+using Plots
+println("=== Heston Model: Monte Carlo with Antithetic Variates ===\n")
+
+# --- Market Inputs ---
+reference_date = Date(2020, 1, 1)
+
+# Heston model parameters
+S0 = 100.0        # Initial spot price
+V0 = 0.04         # Initial variance
+κ = 2.0           # Mean reversion speed
+θ = 0.04          # Long-term variance
+σ = 0.3           # Volatility of variance (vol-of-vol)
+ρ = -0.7          # Correlation between asset and variance processes
+r = 0.05          # Risk-free rate
+
+market_inputs = HestonInputs(reference_date, r, S0, V0, κ, θ, σ, ρ)
+
+# --- Payoff ---
+expiry = reference_date + Year(1)
+strike = S0  # ATM European call
+payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
+
+# --- Pricing problem ---
+prob = PricingProblem(payoff, market_inputs)
+
+# --- Reference price (using Carr-Madan Fourier method) ---
+carr_madan_method = CarrMadan(1.0, 32.0, HestonDynamics())
+carr_madan_solution = solve(prob, carr_madan_method)
+reference_price = carr_madan_solution.price
+
+println("Reference price (Carr-Madan): $reference_price\n")
+
+# --- Function to run Monte Carlo trials ---
+# Create Euler-Maruyama strategy
+trajectories = 5000
+steps = 100
+strategy = EulerMaruyama(trajectories, steps; antithetic=true)
+method = MonteCarlo(HestonDynamics(), strategy)
+sol = solve(prob, method)
+
+first_path = sol.ensemble.solutions[1]
+first_path_u = [first_path.u[i][1] for i in range(1,length(first_path))]
+plot(first_path.t, first_path_u)
+
+anti_path = sol.ensemble.solutions[5001]
+anti_path_u = [anti_path.u[i][1] for i in range(1,length(first_path))]
+plot!(first_path.t, anti_path_u)

examples/montecarlo_heston.jl

@@ -36,7 +36,7 @@ carr_madan_problem = PricingProblem(payoff, market_inputs)
 carr_madan_solution = solve(carr_madan_problem, carr_madan_method)
 
 # --- Monte Carlo Parameters ---
-trajectories = 100_000
+trajectories = 1_000
 steps = 500
 
 # --- Euler–Maruyama ---

examples/montecarlo_heston_antithetic.jl

@@ -0,0 +1,102 @@
+using Revise
+using Hedgehog2
+using Dates
+using Statistics
+using Random
+using Printf
+
+println("=== Heston Model: Monte Carlo with Antithetic Variates ===\n")
+
+# --- Market Inputs ---
+reference_date = Date(2020, 1, 1)
+
+# Heston model parameters
+S0 = 100.0        # Initial spot price
+V0 = 0.04         # Initial variance
+κ = 2.0           # Mean reversion speed
+θ = 0.04          # Long-term variance
+σ = 0.3           # Volatility of variance (vol-of-vol)
+ρ = -0.7          # Correlation between asset and variance processes
+r = 0.05          # Risk-free rate
+
+market_inputs = HestonInputs(reference_date, r, S0, V0, κ, θ, σ, ρ)
+
+# --- Payoff ---
+expiry = reference_date + Year(1)
+strike = S0  # ATM European call
+payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
+
+# --- Pricing problem ---
+prob = PricingProblem(payoff, market_inputs)
+
+# --- Reference price (using Carr-Madan Fourier method) ---
+carr_madan_method = CarrMadan(1.0, 32.0, HestonDynamics())
+carr_madan_solution = solve(prob, carr_madan_method)
+reference_price = carr_madan_solution.price
+
+println("Reference price (Carr-Madan): $reference_price\n")
+
+# --- Function to run Monte Carlo trials ---
+function run_mc_trials(n_trials, trajectories, steps; use_antithetic=false)
+    prices = Float64[]
+    times = Float64[]
+    for trial in 1:n_trials
+        # Create Euler-Maruyama strategy
+        strategy = EulerMaruyama(trajectories, steps; antithetic=use_antithetic)
+        method = MonteCarlo(HestonDynamics(), strategy)
+        
+        # Measure time and solve
+        time = @elapsed begin
+            sol = solve(prob, method)
+            push!(prices, sol.price)
+        end
+        
+        push!(times, time)
+    end
+    
+    return prices, times
+end
+
+# --- Run experiments ---
+n_trials = 20
+trajectories = 5_000
+steps = 100
+
+println("Running experiments with $n_trials trials, $trajectories paths, $steps steps...")
+
+# Standard Monte Carlo
+std_prices, std_times = run_mc_trials(n_trials, trajectories, steps, use_antithetic=false)
+
+# Antithetic Monte Carlo (same number of total paths)
+anti_prices, anti_times = run_mc_trials(n_trials, trajectories ÷ 2, steps, use_antithetic=true)
+
+# --- Compute metrics ---
+std_mean = mean(std_prices)
+std_error = std_mean - reference_price
+std_stdev = std(std_prices)
+std_rmse = sqrt(mean((std_prices .- reference_price).^2))
+std_time = mean(std_times)
+
+anti_mean = mean(anti_prices)
+anti_error = anti_mean - reference_price
+anti_stdev = std(anti_prices)
+anti_rmse = sqrt(mean((anti_prices .- reference_price).^2))
+anti_time = mean(anti_times)
+
+# --- Print results ---
+println("\n=== Results ===")
+println("Metric               | Standard MC       | Antithetic MC      | Improvement")
+println("--------------------+-------------------+--------------------+-------------")
+@printf("Mean price           | %.6f          | %.6f          | -\n", std_mean, anti_mean)
+@printf("Bias                 | %.6f          | %.6f          | %.2fx\n", 
+        std_error, anti_error, abs(std_error)/max(abs(anti_error), 1e-10))
+@printf("Standard deviation   | %.6f          | %.6f          | %.2fx\n", 
+        std_stdev, anti_stdev, std_stdev/anti_stdev)
+@printf("RMSE                 | %.6f          | %.6f          | %.2fx\n", 
+        std_rmse, anti_rmse, std_rmse/anti_rmse)
+@printf("Avg. execution time  | %.3f s           | %.3f s           | %.2fx\n",
+        std_time, anti_time, std_time/anti_time)
+
+# --- Path distribution visualization (code) ---
+println("\nNote: To visualize the distribution of prices across trials, you can plot a histogram")
+println("of the price vectors 'std_prices' and 'anti_prices' to see the variance reduction effect.")

src/market_inputs/market_inputs.jl

@@ -65,9 +65,9 @@ end
 Market data inputs for the Heston stochastic volatility model.
 
 # Fields
-- `referenceDate`: The base date for maturity calculation.
+- `referenceDate`: The base date for maturity calculation (in ticks).
 - `rate`: The risk-free interest rate (annualized).
-- `spot`: The current spot price of the underlying asset.
+- `spot`: The current spot price of the underlying.
 - `V0`: The initial variance of the underlying.
 - `κ`: The rate at which variance reverts to its long-term mean.
 - `θ`: The long-term mean of the variance.
@@ -77,7 +77,7 @@ Market data inputs for the Heston stochastic volatility model.
 Used for pricing under the Heston model and simulation of stochastic volatility paths.
 """
 struct HestonInputs <: AbstractMarketInputs
-    referenceDate
+    referenceDate::Real
     rate::RateCurve
     spot
     V0
@@ -88,12 +88,23 @@ struct HestonInputs <: AbstractMarketInputs
 end
 
 HestonInputs(
-    reference_date,
+    reference_date::TimeType,
+    rate::RateCurve,
+    spot,
+    V0,
+    κ,
+    θ,
+    σ,
+    ρ
+) = HestonInputs(to_ticks(reference_date), rate, spot, V0, κ, θ, σ, ρ)
+
+HestonInputs(
+    reference_date::TimeType,
     rate::Real,
     spot,
     V0,
     κ,
     θ,
     σ,
     ρ
-) = HestonInputs(reference_date, FlatRateCurve(rate), spot, V0, κ, θ, σ, ρ)
+) = HestonInputs(reference_date, FlatRateCurve(rate; reference_date=reference_date), spot, V0, κ, θ, σ, ρ)

src/pricing_methods/carr_madan.jl

@@ -52,6 +52,7 @@ function solve(
     method::CarrMadan
 ) where {C, I <: AbstractMarketInputs}
 
+    println("started carr madan")
     if !is_flat(prob.market.rate)
         throw(ArgumentError("Carr–Madan pricing only supports flat rate curves."))
     end

src/pricing_methods/montecarlo.jl

@@ -59,30 +59,32 @@ 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
+function get_antithetic_ensemble_problem(::PriceDynamics, ::EulerMaruyama, normal_sol::CustomEnsembleSolution)
+    base_prob = normal_sol.solutions[1].prob
 
     antithetic_modify = function (_base_prob, _seed, i)
-        sol = original_sol.solutions[i]
+        sol = normal_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 get_antithetic_ensemble_problem(d::LognormalDynamics , s::BlackScholesExact, normal_sol::CustomEnsembleSolution)
+    tspan = normal_sol.solutions[1].prob.tspan
+    s_flipped = @set s.seeds = normal_sol.seeds
+    m_flipped = @set m.sigma = -m.sigma
+    flipped_problem = sde_problem(d, s_flipped, m_flipped, tspan)
+    return CustomEnsembleProblem(flipped_problem, 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)
@@ -100,21 +102,21 @@ end
 
 function marginal_law(::HestonDynamics, m::HestonInputs, t)
     α = yearfrac(m.rate.reference_date, t)
-    return HestonDistribution(m.spot, m.V0, m.κ, m.θ, m.σ, m.ρ, m.rate, α)
+    rate = zero_rate(m.rate, t)
+    return HestonDistribution(m.spot, m.V0, m.κ, m.θ, m.σ, m.ρ, rate, α)
 end
 
 # ------------------ Ensemble Simulation Wrapper ------------------
 
-function montecarlo_solution(problem::Union{NoiseProblem, SDEProblem}, strategy::S) where {S <: SimulationStrategy}
-    dt = problem.tspan[end] / strategy.steps
+function get_ensemble_problem(problem::Union{NoiseProblem, SDEProblem}, strategy::S) where {S <: SimulationStrategy}
     N = strategy.trajectories
 
     seeds = strategy isa BlackScholesExact && strategy.seeds !== nothing ? strategy.seeds :
             Base.rand(1_000_000_000:2_000_000_000, N)
 
     modify = (p, seed, _) -> remake(p; seed=seed)
     custom_prob = CustomEnsembleProblem(problem, collect(seeds), modify)
-    return solve_custom_ensemble(custom_prob; dt=dt)
+    return custom_prob
 end
 
 # ------------------ Terminal Value Extractors ------------------
@@ -133,21 +135,23 @@ function simulate_paths(method::MonteCarlo, market_inputs::I, T) where {I <: Abs
     strategy = method.strategy
     dynamics = method.dynamics
     tspan = (0.0, T)
+    dt = T / strategy.steps
 
     antithetic = get(strategy.kwargs, :antithetic, false)
 
     # Step 1: simulate original paths
     normal_prob = sde_problem(dynamics, strategy, market_inputs, tspan)
-    normal_sol = montecarlo_solution(normal_prob, strategy)
+    ensemble_prob = get_ensemble_problem(normal_prob, strategy)
+    normal_sol = solve_custom_ensemble(ensemble_prob; dt=dt)
 
     if !antithetic
         return normal_sol
     end
 
     # Step 2: simulate antithetic paths using same seeds, flipped sigma
-    antithetic_prob = antithetic_sde_problem(dynamics, strategy, market_inputs, tspan, normal_sol)
-    antithetic_sol = montecarlo_solution(antithetic_prob, strategy)
-    
+    antithetic_ensemble_prob = get_antithetic_ensemble_problem(dynamics, strategy, normal_sol)
+    antithetic_sol = solve_custom_ensemble(antithetic_ensemble_prob; dt=dt)
+
     # Step 3: combine both solutions
     combined_paths = vcat(normal_sol.solutions, antithetic_sol.solutions)
 

src/solutions/pricing_solutions.jl

@@ -22,8 +22,9 @@ function solve_custom_ensemble(prob::CustomEnsembleProblem; dt, solver=EM(), tra
 
     Threads.@threads for i in 1:N
         seed = seeds[i]
+        pmod = remake(prob.base_problem; seed=seed)
         pmod = prob.modify(prob.base_problem, seed, i)
-        sols[i] = DifferentialEquations.solve(pmod, solver; dt=dt)
+        sols[i] = DifferentialEquations.solve(pmod, solver; dt=dt, save_noise=true)
     end
 
     return CustomEnsembleSolution(sols, seeds)