Skip to content

Commit f09a462

Browse files
aleCombialeCombi
committed
refining antithetic tests
1 parent d6bb75d commit f09a462

9 files changed

+679
-328
lines changed

‎examples/antithetic_analysis_bk_FINAL.jl

Lines changed: 205 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,205 @@
1+
using Hedgehog2
2+
using Dates
3+
using Random
4+
using Statistics
5+
using Plots
6+
7+
function analyze_bk_antithetic_paths(num_paths=10, seed=42)
8+
# Set random seed for reproducibility
9+
Random.seed!(seed)
10+
11+
# Set up common parameters
12+
reference_date = Date(2020, 1, 1)
13+
expiry = reference_date + Year(1)
14+
spot = 100.0
15+
strike = 100.0
16+
17+
# Create payoff
18+
payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
19+
20+
# Heston parameters - can be adjusted to see different behaviors
21+
rate_heston = 0.03
22+
V0 = 0.04
23+
κ = 2.0
24+
θ = 0.04
25+
σ = 0.3
26+
ρ = -0.7
27+
28+
# Create the market inputs
29+
heston_market = HestonInputs(reference_date, rate_heston, spot, V0, κ, θ, σ, ρ)
30+
prob = PricingProblem(payoff, heston_market)
31+
32+
# Generate seeds for paths
33+
path_seeds = rand(1:10^9, num_paths ÷ 2)
34+
35+
# Create the Broadie-Kaya method with antithetic sampling
36+
bk_strategy = HestonBroadieKaya(num_paths ÷ 2, steps=20)
37+
bk_method = MonteCarlo(HestonDynamics(), bk_strategy)
38+
39+
# Update with seeds and antithetic flag
40+
antithetic_method = @set bk_method.strategy.seeds = path_seeds
41+
antithetic_method = @set antithetic_method.strategy.kwargs = merge(antithetic_method.strategy.kwargs, (antithetic=true,))
42+
43+
# Solve the problem
44+
solution = solve(prob, antithetic_method)
45+
46+
# Return the solution for further analysis
47+
return solution, prob
48+
end
49+
50+
function plot_bk_path_pair(solution, prob, path_index=1; apply_exp=true)
51+
half_paths = length(solution.ensemble.solutions) ÷ 2
52+
53+
# Get the original and antithetic paths
54+
original_path = solution.ensemble.solutions[path_index]
55+
antithetic_path = solution.ensemble.solutions[half_paths + path_index]
56+
57+
# Extract time points
58+
time_points = original_path.t
59+
60+
# Extract stock prices (log to normal scale if requested)
61+
if apply_exp
62+
original_prices = [exp(u[1]) for u in original_path.u]
63+
antithetic_prices = [exp(u[1]) for u in antithetic_path.u]
64+
ylabel_text = "Stock Price"
65+
else
66+
original_prices = [u[1] for u in original_path.u]
67+
antithetic_prices = [u[1] for u in antithetic_path.u]
68+
ylabel_text = "Log Price"
69+
end
70+
71+
# Create price plot
72+
p1 = plot(
73+
time_points, original_prices,
74+
label="Original Path",
75+
linewidth=2,
76+
title="Heston Broadie-Kaya: $(apply_exp ? "Stock Price" : "Log Price") Paths",
77+
xlabel="Time (years)",
78+
ylabel=ylabel_text,
79+
legend=:topleft
80+
)
81+
82+
plot!(
83+
p1,
84+
time_points, antithetic_prices,
85+
label="Antithetic Path",
86+
linewidth=2,
87+
linestyle=:dash,
88+
color=:red
89+
)
90+
91+
# Calculate price correlation
92+
price_corr = cor(original_prices, antithetic_prices)
93+
annotate!(p1, [(0.5, minimum(original_prices),
94+
text("Path correlation: $(round(price_corr, digits=4))", 10, :left))])
95+
96+
return p1
97+
end
98+
99+
function calculate_bk_path_statistics(solution, prob)
100+
half_paths = length(solution.ensemble.solutions) ÷ 2
101+
102+
price_correlations = Float64[]
103+
terminal_products = Float64[]
104+
discount_factor = df(prob.market.rate, prob.payoff.expiry)
105+
106+
original_terminals = Float64[]
107+
antithetic_terminals = Float64[]
108+
109+
# Calculate statistics for each path pair
110+
for i in 1:half_paths
111+
original_path = solution.ensemble.solutions[i]
112+
antithetic_path = solution.ensemble.solutions[half_paths + i]
113+
114+
# Extract full paths for correlation
115+
original_prices = [exp(u[1]) for u in original_path.u]
116+
antithetic_prices = [exp(u[1]) for u in antithetic_path.u]
117+
118+
# Calculate correlation between full paths
119+
path_corr = cor(original_prices, antithetic_prices)
120+
push!(price_correlations, path_corr)
121+
122+
# Extract terminal values
123+
original_terminal = exp(original_path.u[end][1])
124+
antithetic_terminal = exp(antithetic_path.u[end][1])
125+
126+
# Calculate product of terminal values
127+
push!(terminal_products, original_terminal * antithetic_terminal)
128+
129+
# Store terminal values for payoff calculation
130+
push!(original_terminals, original_terminal)
131+
push!(antithetic_terminals, antithetic_terminal)
132+
end
133+
134+
# Calculate payoffs
135+
original_payoffs = prob.payoff.(original_terminals)
136+
antithetic_payoffs = prob.payoff.(antithetic_terminals)
137+
138+
# Calculate standard MC estimator variance (using original paths)
139+
standard_estimator = discount_factor * original_payoffs
140+
standard_variance = var(standard_estimator)
141+
142+
# Calculate antithetic MC estimator variance
143+
antithetic_pairs = [(original_payoffs[i] + antithetic_payoffs[i]) / 2 for i in 1:half_paths]
144+
antithetic_estimator = discount_factor * antithetic_pairs
145+
antithetic_variance = var(antithetic_estimator)
146+
147+
# Calculate payoff correlation
148+
payoff_correlation = cor(original_payoffs, antithetic_payoffs)
149+
150+
# Calculate variance reduction
151+
var_reduction_ratio = standard_variance / antithetic_variance
152+
153+
return (
154+
mean_path_correlation = mean(price_correlations),
155+
path_correlations = price_correlations,
156+
mean_terminal_product = mean(terminal_products),
157+
payoff_correlation = payoff_correlation,
158+
standard_variance = standard_variance,
159+
antithetic_variance = antithetic_variance,
160+
var_reduction_ratio = var_reduction_ratio
161+
)
162+
end
163+
164+
# Run analysis and generate plots
165+
num_paths = 100 # Use an even number
166+
bk_solution, bk_prob = analyze_bk_antithetic_paths(num_paths, 42)
167+
168+
# Plot several path pairs to observe different behaviors
169+
p1 = plot_bk_path_pair(bk_solution, bk_prob, 1)
170+
p2 = plot_bk_path_pair(bk_solution, bk_prob, 2)
171+
p3 = plot_bk_path_pair(bk_solution, bk_prob, 3)
172+
p4 = plot_bk_path_pair(bk_solution, bk_prob, 4)
173+
174+
# Plot them in a 2x2 grid
175+
path_plots = plot(p1, p2, p3, p4, layout=(2,2), size=(900, 700), title="Broadie-Kaya Antithetic Path Examples")
176+
display(path_plots)
177+
178+
# Also plot some in log scale
179+
p1_log = plot_bk_path_pair(bk_solution, bk_prob, 1, apply_exp=false)
180+
p2_log = plot_bk_path_pair(bk_solution, bk_prob, 2, apply_exp=false)
181+
182+
log_plots = plot(p1_log, p2_log, layout=(1,2), size=(900, 400), title="Broadie-Kaya Paths (Log Scale)")
183+
display(log_plots)
184+
185+
# Calculate statistics across all paths
186+
stats = calculate_bk_path_statistics(bk_solution, bk_prob)
187+
188+
# Print summary statistics
189+
println("=== Broadie-Kaya Antithetic Path Analysis ===")
190+
println("Number of paths: $num_paths")
191+
println("Average path correlation: $(round(stats.mean_path_correlation, digits=4))")
192+
println("Payoff correlation: $(round(stats.payoff_correlation, digits=4))")
193+
println("Standard estimator variance: $(round(stats.standard_variance, digits=8))")
194+
println("Antithetic estimator variance: $(round(stats.antithetic_variance, digits=8))")
195+
println("Variance reduction ratio: $(round(stats.var_reduction_ratio, digits=2))×")
196+
197+
# Plot distribution of path correlations
198+
histogram!(
199+
stats.path_correlations,
200+
bins=20,
201+
title="Distribution of Broadie-Kaya Path Correlations",
202+
xlabel="Correlation",
203+
ylabel="Frequency",
204+
legend=false
205+
)

‎examples/antithetic_analysis_plot_FINAL.jl

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
using Hedgehog2
2+
using Dates
3+
using Random
4+
using Statistics
5+
using Plots
6+
7+
function analyze_antithetic_paths(prob::PricingProblem, mc_method::MonteCarlo;
8+
num_paths=10, seed=nothing)
9+
# Set random seed if provided
10+
if seed !== nothing
11+
Random.seed!(seed)
12+
end
13+
14+
# Generate seeds for paths
15+
path_seeds = rand(1:10^9, num_paths ÷ 2)
16+
17+
# Create the antithetic method
18+
antithetic_method = @set mc_method.strategy.seeds = path_seeds
19+
antithetic_method = @set antithetic_method.strategy.kwargs = merge(antithetic_method.strategy.kwargs, (antithetic=true,))
20+
21+
# Solve the problem
22+
solution = solve(prob, antithetic_method)
23+
24+
# Return the ensemble solutions directly
25+
return solution.ensemble.solutions
26+
end
27+
28+
function plot_path_pair(paths, dynamics, strategy, method_name; apply_exp=true)
29+
original_path = paths[1]
30+
antithetic_path = paths[length(paths)÷2 + 1]
31+
32+
# Extract time points
33+
time_points = original_path.t
34+
35+
# Check if paths contain vector states (like Heston) or scalar states (like some BS)
36+
is_vector_state = eltype(original_path.u) <: AbstractVector
37+
38+
if is_vector_state
39+
# Handle paths with vector states (e.g., Heston model)
40+
if apply_exp
41+
original_prices = [exp(u[1]) for u in original_path.u]
42+
antithetic_prices = [exp(u[1]) for u in antithetic_path.u]
43+
else
44+
original_prices = [u[1] for u in original_path.u]
45+
antithetic_prices = [u[1] for u in antithetic_path.u]
46+
end
47+
48+
# Extract variance components if they exist
49+
has_variance = length(original_path.u[1]) >= 2
50+
51+
if has_variance
52+
original_vars = [u[2] for u in original_path.u]
53+
antithetic_vars = [u[2] for u in antithetic_path.u]
54+
var_corr = cor(original_vars, antithetic_vars)
55+
end
56+
else
57+
# Handle scalar state paths
58+
if apply_exp
59+
original_prices = exp.(original_path.u)
60+
antithetic_prices = exp.(antithetic_path.u)
61+
else
62+
original_prices = original_path.u
63+
antithetic_prices = antithetic_path.u
64+
end
65+
has_variance = false
66+
end
67+
68+
# Create plot
69+
p = plot(
70+
time_points, original_prices,
71+
label="Original Path",
72+
linewidth=2,
73+
title="$method_name: $(apply_exp ? "Stock Price" : "Log Price") Paths",
74+
xlabel="Time (years)",
75+
ylabel=apply_exp ? "Stock Price" : "Log Price",
76+
legend=:topleft
77+
)
78+
79+
plot!(
80+
p,
81+
time_points, antithetic_prices,
82+
label="Antithetic Path",
83+
linewidth=2,
84+
linestyle=:dash,
85+
color=:red
86+
)
87+
88+
# Calculate correlation
89+
price_corr = cor(original_prices, antithetic_prices)
90+
91+
# Add correlation annotation (only for price paths, ignoring variance)
92+
annotate!(p, [(0.5, minimum(original_prices),
93+
text("Path correlation: $(round(price_corr, digits=4))", 10, :left))])
94+
95+
return p
96+
end
97+
98+
# Set up common parameters
99+
reference_date = Date(2020, 1, 1)
100+
expiry = reference_date + Year(1)
101+
spot = 10.0
102+
strike = 10.0
103+
seed = 42
104+
105+
# Create payoff
106+
payoff = VanillaOption(strike, expiry, European(), Call(), Spot())
107+
108+
# Black-Scholes setup
109+
rate_bs = 0.05
110+
sigma_bs = 0.20
111+
bs_market = BlackScholesInputs(reference_date, rate_bs, spot, sigma_bs)
112+
bs_prob = PricingProblem(payoff, bs_market)
113+
114+
# Heston setup
115+
rate_heston = 0.03
116+
V0 = 0.04
117+
κ = 2.0
118+
θ = 0.04
119+
σ = 0.3
120+
ρ = -0.7
121+
heston_market = HestonInputs(reference_date, rate_heston, spot, V0, κ, θ, σ, ρ)
122+
heston_prob = PricingProblem(payoff, heston_market)
123+
124+
# Create Monte Carlo methods
125+
num_paths = 10
126+
steps_bs = 100
127+
steps_heston = 100
128+
129+
# 1. Black-Scholes with Exact Simulation
130+
bs_exact_strategy = BlackScholesExact(num_paths, steps=steps_bs)
131+
bs_exact_method = MonteCarlo(LognormalDynamics(), bs_exact_strategy)
132+
bs_exact_paths = analyze_antithetic_paths(bs_prob, bs_exact_method, seed=seed)
133+
134+
# 2. Black-Scholes with Euler-Maruyama
135+
bs_em_strategy = EulerMaruyama(num_paths, steps=steps_bs)
136+
bs_em_method = MonteCarlo(LognormalDynamics(), bs_em_strategy)
137+
bs_em_paths = analyze_antithetic_paths(bs_prob, bs_em_method, seed=seed)
138+
139+
# 3. Heston with Euler-Maruyama
140+
heston_em_strategy = EulerMaruyama(num_paths, steps=steps_heston)
141+
heston_em_method = MonteCarlo(HestonDynamics(), heston_em_strategy)
142+
heston_em_paths = analyze_antithetic_paths(heston_prob, heston_em_method, seed=seed)
143+
144+
# 4. Heston with Broadie-Kaya (for fewer paths due to computational complexity)
145+
heston_bk_strategy = HestonBroadieKaya(num_paths, steps=100)
146+
heston_bk_method = MonteCarlo(HestonDynamics(), heston_bk_strategy)
147+
heston_bk_paths = analyze_antithetic_paths(heston_prob, heston_bk_method, seed=seed)
148+
149+
# Generate plots with exponential transformation (natural price scale)
150+
p1 = plot_path_pair(bs_exact_paths, LognormalDynamics(), bs_exact_strategy, "Black-Scholes Exact", apply_exp=false)
151+
p2 = plot_path_pair(bs_em_paths, LognormalDynamics(), bs_em_strategy, "Black-Scholes Euler-Maruyama", apply_exp=true)
152+
p3 = plot_path_pair(heston_em_paths, HestonDynamics(), heston_em_strategy, "Heston Euler-Maruyama", apply_exp=false)
153+
p4 = plot_path_pair(heston_bk_paths, HestonDynamics(), heston_bk_strategy, "Heston Broadie-Kaya", apply_exp=true)
154+
155+
# If you want to view paths in log scale, use apply_exp=false
156+
# For example:
157+
# p1_log = plot_path_pair(bs_exact_paths, LognormalDynamics(), bs_exact_strategy, "Black-Scholes Exact (Log Scale)", apply_exp=false)
158+
159+
# Plot them in a 2x2 grid
160+
final_plot = plot(p1, p2, p3, p4, layout=(2,2), size=(900, 700))
161+
display(final_plot)

0 commit comments

Comments
 (0)