Adjusted FD for relative bumps

Author: aleCombi

examples/comparisons/american.jl

@@ -36,16 +36,16 @@ steps_lsm = 100
 seed = 42
 
 # Create deterministic seeds for reproducibility
-path_seeds = rand(Xoshiro(seed), UInt64, trajectories)
-lsm_config = SimulationConfig(trajectories, steps=steps_lsm)
+seeds = Base.rand(UInt64, trajectories)
+lsm_config = SimulationConfig(trajectories, steps=steps_lsm, variance_reduction=Hedgehog2.Antithetic(), seeds=seeds)
 degree = 5  # Polynomial degree for regression
 lsm_method = LSM(dynamics, BlackScholesExact(), lsm_config, degree)
 
 # ------------------------------
 # Define lenses for Greeks
 # ------------------------------
 vol_lens = VolLens(1,1)
-spot_lens = @optic _.market_inputs.spot
+spot_lens = Hedgehog2.SpotLens()
 rate_lens = ZeroRateSpineLens(1)
 lenses = (spot_lens, vol_lens, rate_lens)
 
@@ -58,8 +58,9 @@ df_put = run_model_comparison_table(
     [crr_method, lsm_method],
     lenses;
     ad_method = ForwardAD(),
-    fd_method = FiniteDifference(1e-4),
+    fd_method = FiniteDifference(1e-2),
     analytic_method = nothing,  # No analytic for American
+    use_belapsed = true
 )
 println(df_put)
 
@@ -84,7 +85,9 @@ df_call_div = run_model_comparison_table(
     [crr_method, lsm_method],
     lenses;
     ad_method = ForwardAD(),
-    fd_method = FiniteDifference(1e-4),
-    analytic_method = nothing # No analytic for American
+    fd_method = FiniteDifference(1e-2),
+    analytic_method = nothing, # No analytic for American
+    use_belapsed = true
+
 )
 println(df_call_div)

examples/polished_examples/lsm.jl

@@ -0,0 +1,48 @@
+using Revise, Hedgehog2, BenchmarkTools, Dates
+using Accessors
+import Accessors: @optic
+using Test
+using DataFrames
+using Random
+
+# ------------------------------
+# Define payoff and pricing problem
+# ------------------------------
+strike = 100.0
+reference_date = Date(2020, 1, 1)
+expiry = reference_date + Year(1)
+rate = 0.05
+spot = 100.0
+sigma = 0.2
+
+# --- American Put Option ---
+american_put = VanillaOption(strike, expiry, American(), Put(), Spot())
+market_inputs = BlackScholesInputs(reference_date, rate, spot, sigma)
+put_prob = PricingProblem(american_put, market_inputs)
+
+# Create deterministic seeds for reproducibility
+seeds = Base.rand(UInt64, trajectories)
+lsm_config = SimulationConfig(trajectories, steps=steps_lsm, variance_reduction=Hedgehog2.Antithetic(), seeds=seeds)
+degree = 5  # Polynomial degree for regression
+lsm_method = LSM(dynamics, BlackScholesExact(), lsm_config, degree)
+
+# ------------------------------
+# Define lenses for Greeks
+# ------------------------------
+vol_lens = VolLens(1,1)
+spot_lens = Hedgehog2.SpotLens()
+rate_lens = ZeroRateSpineLens(1)
+lenses = (spot_lens, vol_lens, rate_lens)
+
+delta_prob = GreekProblem(put_prob, spot_lens)
+
+ad_sol = solve(delta_prob, ForwardAD(), lsm_method)
+fd_sol = solve(delta_prob, FiniteDifference(1), lsm_method)
+
+price_sol = solve(put_prob, lsm_method).price
+
+market_inputs_bumped = BlackScholesInputs(reference_date, rate, spot + 1, sigma)
+put_prob_bumped = PricingProblem(american_put, market_inputs_bumped)
+
+price_sol_bumped = solve(put_prob_bumped, lsm_method).price
+

src/greeks/greeks_problem.jl

@@ -46,29 +46,29 @@ end
 
 function compute_fd_derivative(::FDForward, prob, lens, ε, pricing_method)
     x₀ = lens(prob)
-    prob_up = set(prob, lens, x₀ + ε)
+    prob_up = set(prob, lens, x₀ * (1 + ε))
 
     v_up = solve(prob_up, pricing_method).price
     v₀ = solve(prob, pricing_method).price
-    return (v_up - v₀) / ε
+    return (v_up - v₀) / (x₀ * ε)
 end
 
 function compute_fd_derivative(::FDBackward, prob, lens, ε, pricing_method)
     x₀ = lens(prob)
-    prob_down = set(prob, lens, x₀ - ε)
+    prob_down = set(prob, lens, x₀ * (1 - ε))
     v_down = solve(prob_down, pricing_method).price
     v₀ = solve(prob, pricing_method).price
-    return (v₀ - v_down) / ε
+    return (v₀ - v_down) / (x₀ * ε)
 end
 
 function compute_fd_derivative(::FDCentral, prob, lens, ε, pricing_method)
     x₀ = lens(prob)
 
-    prob_up = set(prob, lens, x₀ + ε)
-    prob_down = set(prob, lens, x₀ - ε)
+    prob_up = set(prob, lens, x₀ * (1 + ε))
+    prob_down = set(prob, lens, x₀ * (1 - ε))
     v_up = solve(prob_up, pricing_method).price
     v_down = solve(prob_down, pricing_method).price
-    return (v_up - v_down) / (2ε)
+    return (v_up - v_down) / (2ε * x₀)
 end
 
 function solve(