Added test for antithetic variates

Author: aleCombi

test/antithetic_variates.jl

@@ -0,0 +1,227 @@
+using Hedgehog2
+using Dates
+using Random
+using Statistics
+using Test
+using Printf
+
+function analyze_path_correlation(paths, n_original)
+    # Calculate correlations between original and antithetic paths
+    correlations = Float64[]
+    
+    for i in 1:n_original
+        original_path = paths[i]
+        antithetic_path = paths[i + n_original]
+        
+        # Get price values at each time step
+        if eltype(original_path.u) <: AbstractVector
+            original_prices = [u[1] for u in original_path.u]
+            antithetic_prices = [u[1] for u in antithetic_path.u]
+        else
+            original_prices = original_path.u
+            antithetic_prices = antithetic_path.u
+        end
+        
+        # Calculate correlation
+        path_correlation = cor(original_prices, antithetic_prices)
+        push!(correlations, path_correlation)
+    end
+    
+    return correlations
+end
+
+function calculate_return_correlation(paths, n_original)
+    # Calculate correlations between log returns of original and antithetic paths
+    return_correlations = Float64[]
+    
+    for i in 1:n_original
+        original_path = paths[i]
+        antithetic_path = paths[i + n_original]
+        
+        # Get price values at each time step
+        if eltype(original_path.u) <: AbstractVector
+            original_prices = [u[1] for u in original_path.u]
+            antithetic_prices = [u[1] for u in antithetic_path.u]
+        else
+            original_prices = original_path.u
+            antithetic_prices = antithetic_path.u
+        end
+        
+        # Calculate log returns
+        original_returns = diff(log.(original_prices))
+        antithetic_returns = diff(log.(antithetic_prices))
+        
+        # Skip if either array has issues
+        if any(isnan.(original_returns)) || any(isnan.(antithetic_returns)) ||
+           any(isinf.(original_returns)) || any(isinf.(antithetic_returns))
+            continue
+        end
+        
+        # Calculate correlation of returns
+        return_corr = cor(original_returns, antithetic_returns)
+        push!(return_correlations, return_corr)
+    end
+    
+    return return_correlations
+end
+
+@testset "Antithetic Path Correlation Tests" begin
+    # Common parameters
+    Random.seed!(42)
+    trajectories = 100
+    steps = 100
+    reference_date = Date(2020, 1, 1)
+    expiry = reference_date + Year(1)
+    spot = 100.0
+    strike = 100.0
+
+    # Create European call option payoff
+    payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
+    
+    @testset "Black-Scholes Model" begin
+        # BS parameters
+        rate_bs = 0.05
+        sigma_bs = 0.20
+        
+        # Create BS market inputs
+        bs_market = BlackScholesInputs(reference_date, rate_bs, spot, sigma_bs)
+        bs_prob = PricingProblem(payoff, bs_market)
+        
+        # Create Monte Carlo with antithetic sampling for BS
+        bs_seeds = rand(1:10^9, trajectories)
+        bs_strategy = BlackScholesExact(trajectories ÷ 2, steps, seeds=bs_seeds[1:trajectories÷2], antithetic=true)
+        bs_method = MonteCarlo(LognormalDynamics(), bs_strategy)
+        
+        # Solve and get paths
+        bs_solution = solve(bs_prob, bs_method)
+        bs_paths = bs_solution.ensemble.solutions
+        
+        # Analyze correlations
+        bs_correlations = analyze_path_correlation(bs_paths, trajectories ÷ 2)
+        bs_return_correlations = calculate_return_correlation(bs_paths, trajectories ÷ 2)
+        
+        # Report statistics for diagnostic purposes
+        println("Black-Scholes path correlations:")
+        println("  Mean: ", mean(bs_correlations))
+        println("  Min: ", minimum(bs_correlations))
+        println("  Max: ", maximum(bs_correlations))
+        
+        println("\nBlack-Scholes return correlations:")
+        println("  Mean: ", mean(bs_return_correlations))
+        println("  Min: ", minimum(bs_return_correlations))
+        println("  Max: ", maximum(bs_return_correlations))
+        
+        # Test that price paths show strong negative correlation
+        @test mean(bs_correlations) < -0.7
+        
+        # Test that log returns show near-perfect negative correlation
+        @test mean(bs_return_correlations) < -0.95
+        
+        # Test product of terminal values is consistent with theoretical value
+        bs_original = bs_paths[1]
+        bs_antithetic = bs_paths[trajectories ÷ 2 + 1]
+        
+        if eltype(bs_original.u) <: AbstractVector
+            bs_orig_terminal = bs_original.u[end][1]
+            bs_anti_terminal = bs_antithetic.u[end][1]
+        else
+            bs_orig_terminal = bs_original.u[end]
+            bs_anti_terminal = bs_antithetic.u[end]
+        end
+        
+        terminal_product = bs_orig_terminal * bs_anti_terminal
+        theoretical_product = spot^2 * exp(2 * rate_bs * yearfrac(reference_date, expiry))
+        
+        # The product should be roughly close to the theoretical value
+        # Using a large tolerance due to randomness
+        @test isapprox(terminal_product, theoretical_product, rtol=0.2)
+    end
+    
+    @testset "Heston Model" begin
+        # Heston parameters
+        rate_heston = 0.03
+        V0 = 0.04
+        κ = 2.0
+        θ = 0.04
+        σ = 0.3
+        ρ = -0.7
+        
+        # Create Heston market inputs
+        heston_market = HestonInputs(reference_date, rate_heston, spot, V0, κ, θ, σ, ρ)
+        heston_prob = PricingProblem(payoff, heston_market)
+        
+        # Create Monte Carlo with antithetic sampling for Heston
+        heston_seeds = rand(1:10^9, trajectories)
+        heston_strategy = EulerMaruyama(trajectories ÷ 2, steps, seeds=heston_seeds[1:trajectories÷2], antithetic=true)
+        heston_method = MonteCarlo(HestonDynamics(), heston_strategy)
+        
+        # Solve and get paths
+        heston_solution = solve(heston_prob, heston_method)
+        heston_paths = heston_solution.ensemble.solutions
+        
+        # Analyze correlations
+        heston_correlations = analyze_path_correlation(heston_paths, trajectories ÷ 2)
+        heston_return_correlations = calculate_return_correlation(heston_paths, trajectories ÷ 2)
+        
+        # Report statistics for diagnostic purposes
+        println("\nHeston path correlations:")
+        println("  Mean: ", mean(heston_correlations))
+        println("  Min: ", minimum(heston_correlations))
+        println("  Max: ", maximum(heston_correlations))
+        
+        println("\nHeston return correlations:")
+        println("  Mean: ", mean(heston_return_correlations))
+        println("  Min: ", minimum(heston_return_correlations))
+        println("  Max: ", maximum(heston_return_correlations))
+        
+        # For Heston, the correlation might not be as strong due to stochastic volatility
+        @test mean(heston_correlations) < -0.5
+        
+        # For return correlations, we still expect strong negative correlation
+        @test mean(heston_return_correlations) < -0.7
+        
+        # Test percentage of negative correlations
+        # In Heston, we expect most but not necessarily all paths to show negative correlation
+        percent_negative = sum(heston_correlations .< 0) / length(heston_correlations) * 100
+        println("  Percentage of negative correlations: $(percent_negative)%")
+        @test percent_negative > 90  # At least 90% should be negative
+    end
+    
+    @testset "Correlation Comparison" begin
+        # Compare correlations between BS and Heston to verify both are working
+        # but may have different characteristics
+        
+        # Repeat setup briefly
+        rate_bs = 0.05
+        sigma_bs = 0.20
+        bs_market = BlackScholesInputs(reference_date, rate_bs, spot, sigma_bs)
+        bs_prob = PricingProblem(payoff, bs_market)
+        bs_seeds = rand(1:10^9, trajectories)
+        bs_strategy = BlackScholesExact(trajectories ÷ 2, steps, seeds=bs_seeds[1:trajectories÷2], antithetic=true)
+        bs_method = MonteCarlo(LognormalDynamics(), bs_strategy)
+        bs_solution = solve(bs_prob, bs_method)
+        bs_paths = bs_solution.ensemble.solutions
+        
+        rate_heston = 0.03
+        V0 = 0.04
+        κ = 2.0
+        θ = 0.04
+        σ = 0.3
+        ρ = -0.7
+        heston_market = HestonInputs(reference_date, rate_heston, spot, V0, κ, θ, σ, ρ)
+        heston_prob = PricingProblem(payoff, heston_market)
+        heston_seeds = rand(1:10^9, trajectories)
+        heston_strategy = EulerMaruyama(trajectories ÷ 2, steps, seeds=heston_seeds[1:trajectories÷2], antithetic=true)
+        heston_method = MonteCarlo(HestonDynamics(), heston_strategy)
+        heston_solution = solve(heston_prob, heston_method)
+        heston_paths = heston_solution.ensemble.solutions
+        
+        # Get path correlations
+        bs_correlations = analyze_path_correlation(bs_paths, trajectories ÷ 2)
+        heston_correlations = analyze_path_correlation(heston_paths, trajectories ÷ 2)
+        
+        # Black-Scholes should typically have stronger negative correlation
+        # due to simpler model dynamics
+        @test mean(bs_correlations) < mean(heston_correlations)
+    end
+end

test/montecarlo_heston.jl

@@ -0,0 +1,134 @@
+using Test
+using Hedgehog2
+using Dates
+using Random
+using Statistics
+using Printf
+
+@testset "Heston Model Monte Carlo Test" begin
+    # --------------------------
+    # Setup common test parameters
+    # --------------------------
+    # Heston model parameters
+    spot = 100.0
+    strike = 100.0
+    reference_date = Date(2020, 1, 1)
+    expiry = reference_date + Year(1)
+    r = 0.03             # Risk-free rate
+    V0 = 0.04            # Initial variance
+    κ = 2.0              # Mean reversion speed
+    θ = 0.04             # Long-term variance
+    σ = 0.3              # Vol-of-vol
+    ρ = -0.7             # Correlation
+    
+    # Print parameter configuration for debugging
+    println("Test Parameters:")
+    println("  Spot: $spot")
+    println("  Strike: $strike")
+    println("  Rate: $r")
+    println("  Initial variance (V0): $V0")
+    println("  Mean reversion (κ): $κ")
+    println("  Long-term variance (θ): $θ")
+    println("  Vol-of-vol (σ): $σ")
+    println("  Correlation (ρ): $ρ")
+    println("  Time to maturity: 1 year")
+    println("  Option type: European Call")
+    
+    # Define payoff and market inputs
+    payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
+    market_inputs = HestonInputs(reference_date, r, spot, V0, κ, θ, σ, ρ)
+    prob = PricingProblem(payoff, market_inputs)
+    
+    # Define test parameters - increase trajectories for more reliable results
+    trajectories = 5_000
+    steps = 100
+    
+    # Reference price using Carr-Madan (Fourier) method
+    carr_madan_method = CarrMadan(1.0, 32.0, HestonDynamics())
+    carr_madan_sol = solve(prob, carr_madan_method)
+    reference_price = carr_madan_sol.price
+    
+    # Print reference price
+    println("Reference price (Carr-Madan): $reference_price")
+    
+    # Test scenarios for Euler-Maruyama
+    scenarios = [
+        ("Euler-Maruyama without antithetic", 
+         EulerMaruyama(trajectories, steps, seeds=nothing), 
+         HestonDynamics()),
+         
+        ("Euler-Maruyama with antithetic", 
+         EulerMaruyama(trajectories ÷ 2, steps, antithetic=true, seeds=nothing), 
+         HestonDynamics())
+    ]
+    
+    results = Dict()
+
+    # Run all scenarios
+    for (scenario_name, strategy, dynamics) in scenarios
+        # Print current scenario
+        println("Running scenario: ", scenario_name)
+        
+        # Execute multiple trials to measure variance
+        prices = Float64[]
+        for trial in 1:5
+            # Create a new RNG with a trial-specific seed
+            trial_rng = MersenneTwister(42 + trial)
+            
+            # Generate random seeds for each path
+            trial_seeds = rand(trial_rng, 1:10^9, trajectories)
+
+            # Create modified strategy with trial seeds
+            modified_strategy = @set strategy.seeds = trial_seeds
+
+            # Create Monte Carlo method
+            mc_method = MonteCarlo(dynamics, modified_strategy)
+            
+            # Solve with current seed
+            sol = solve(prob, mc_method)
+            push!(prices, sol.price)
+            
+            # Print diagnostic info
+            println("  Trial $trial: Price = $(sol.price)")
+        end
+        
+        # Calculate statistics
+        mean_price = mean(prices)
+        price_variance = var(prices)
+        error = abs(mean_price - reference_price)
+        rel_error = error / reference_price
+        
+        # Store results
+        results[scenario_name] = (
+            mean_price = mean_price,
+            reference_price = reference_price,
+            error = error,
+            rel_error = rel_error,
+            variance = price_variance
+        )
+        
+        # Test the result with a generous tolerance (Heston MC typically needs more paths)
+        @test isapprox(mean_price, reference_price, rtol=0.05)
+    end
+    
+    # Compare variance reduction
+    if length(scenarios) >= 2
+        std_var = results["Euler-Maruyama without antithetic"].variance
+        anti_var = results["Euler-Maruyama with antithetic"].variance
+        var_reduction = std_var / anti_var
+        
+        @info "Variance reduction (Euler-Maruyama): $(var_reduction)×"
+        @test var_reduction > 1.0  # Antithetic should reduce variance
+    end
+    
+    # Print summary
+    println("\n=== Heston Model Monte Carlo Test Results ===")
+    println("Reference price (Carr-Madan): $reference_price")
+    println("Scenario                    | Price      | Rel Error | Variance")
+    println("----------------------------|------------|-----------|----------")
+    
+    for (name, res) in results
+        @printf("%-28s | %.6f | %.6f | %.2e\n", 
+                name, res.mean_price, res.rel_error, res.variance)
+    end
+end

test/runtests.jl

@@ -7,4 +7,6 @@ include("least_squares_montecarlo.jl")
 include("implied_vol.jl")
 include("rate_curve.jl")
 include("greeks.jl")
-include("montecarlo_black_scholes.jl")
+include("montecarlo_black_scholes.jl")
+include("montecarlo_heston.jl")
+include("antithetic_variates.jl")