Skip to content

Commit 593ccd0

Browse files
aleCombialeCombi
committed
Working on Mc pathwise derivative in examples
1 parent c5c097c commit 593ccd0

File tree

4 files changed

+349
-0
lines changed

4 files changed

+349
-0
lines changed

‎examples/mc_black_scholes_vega.jl

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
using DifferentialEquations
2+
using DiffEqNoiseProcess
3+
using ForwardDiff
4+
using Statistics
5+
using Distributions
6+
7+
# --------------------------
8+
# Custom Ensemble Framework
9+
# --------------------------
10+
11+
struct CustomEnsembleProblem{P, F}
12+
base_problem::P
13+
seeds::Vector{Int64}
14+
modify::F # (base_problem, seed, index) -> new_problem
15+
end
16+
17+
struct CustomEnsembleSolution{S}
18+
solutions::Vector{S}
19+
seeds::Vector{Int64}
20+
end
21+
22+
function solve_custom_ensemble(prob::CustomEnsembleProblem; dt, solver=EM(), trajectories=nothing)
23+
N = trajectories === nothing ? length(prob.seeds) : trajectories
24+
seeds = length(prob.seeds) == N ? prob.seeds : collect(1:N)
25+
26+
sols = Vector{Any}(undef, N)
27+
28+
Threads.@threads for i in 1:N
29+
seed = seeds[i]
30+
pmod = prob.modify(prob.base_problem, seed, i)
31+
sols[i] = DifferentialEquations.solve(pmod, solver; dt=dt)
32+
end
33+
34+
return CustomEnsembleSolution(sols, seeds)
35+
end
36+
37+
# --------------------------
38+
# GBM with NoiseProcess
39+
# --------------------------
40+
41+
function build_noise_gbm_ensemble(μ, σ, S0, tspan; seeds=nothing)
42+
base_proc = GeometricBrownianMotionProcess(μ, σ, tspan[1], S0)
43+
base_prob = NoiseProblem(base_proc, tspan)
44+
45+
N = seeds === nothing ? 100 : length(seeds)
46+
seeds = seeds === nothing ? collect(1:N) : seeds
47+
modify = (prob, seed, i) -> remake(prob; seed=seed)
48+
return CustomEnsembleProblem(base_prob, collect(seeds), modify)
49+
end
50+
51+
# --------------------------
52+
# Monte Carlo Vega
53+
# --------------------------
54+
55+
function mc_black_scholes_vega(; μ=0.0, σ=0.2, S0=100.0, K=100.0, r=0.01, T=1.0, ε=1e-4, N=10_000, dt=1/250)
56+
tspan = (0.0, T)
57+
58+
base_ens = build_noise_gbm_ensemble(μ, σ, S0, tspan; seeds=1:N)
59+
bump_ens = build_noise_gbm_ensemble(μ, σ + ε, S0, tspan; seeds=1:N)
60+
61+
sol_base = solve_custom_ensemble(base_ens; dt=dt)
62+
sol_bump = solve_custom_ensemble(bump_ens; dt=dt)
63+
64+
payoffs_base = max.(last.(map(s -> s.u, sol_base.solutions)) .- K, 0.0)
65+
payoffs_bump = max.(last.(map(s -> s.u, sol_bump.solutions)) .- K, 0.0)
66+
67+
df = exp(-r * T)
68+
vega_est = mean(payoffs_bump .- payoffs_base) / ε * df
69+
70+
return vega_est
71+
end
72+
73+
# --------------------------
74+
# Analytic Black-Scholes Vega
75+
# --------------------------
76+
77+
function bs_vega_analytic(S, K, r, σ, T)
78+
d1 = (log(S / K) + (r + 0.5 * σ^2) * T) / (σ * sqrt(T))
79+
return S * sqrt(T) * pdf(Normal(), d1)
80+
end
81+
82+
# --------------------------
83+
# Run Comparison
84+
# --------------------------
85+
86+
S0 = 100.0
87+
K = 100.0
88+
r = 0.01
89+
T = 1.0
90+
σ = 0.2
91+
N = 10_000
92+
93+
vega_mc = mc_black_scholes_vega(; σ=σ, S0=S0, K=K, r=r, T=T, N=N)
94+
vega_an = bs_vega_analytic(S0, K, r, σ, T)
95+
96+
println("Monte Carlo Vega: ", vega_mc)
97+
println("Analytic Vega: ", vega_an)
98+
println("Relative Error: ", abs(vega_mc - vega_an) / vega_an)

‎examples/mc_black_scholes_vega_ad.jl

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
using DifferentialEquations
2+
using DiffEqNoiseProcess
3+
using ForwardDiff
4+
using Statistics
5+
using Distributions
6+
7+
# --------------------------
8+
# Custom Ensemble Framework
9+
# --------------------------
10+
11+
struct CustomEnsembleProblem{P, F}
12+
base_problem::P
13+
seeds::Vector{Int64}
14+
modify::F # (base_problem, seed, index) -> new_problem
15+
end
16+
17+
struct CustomEnsembleSolution{S}
18+
solutions::Vector{S}
19+
seeds::Vector{Int64}
20+
end
21+
22+
function solve_custom_ensemble(prob::CustomEnsembleProblem; dt, solver=EM(), trajectories=nothing)
23+
N = trajectories === nothing ? length(prob.seeds) : trajectories
24+
seeds = length(prob.seeds) == N ? prob.seeds : collect(1:N)
25+
26+
sols = Vector{Any}(undef, N)
27+
28+
Threads.@threads for i in 1:N
29+
seed = seeds[i]
30+
pmod = prob.modify(prob.base_problem, seed, i)
31+
sols[i] = DifferentialEquations.solve(pmod, solver; dt=dt)
32+
end
33+
34+
return CustomEnsembleSolution(sols, seeds)
35+
end
36+
37+
# --------------------------
38+
# GBM with NoiseProcess
39+
# --------------------------
40+
41+
function build_noise_gbm_ensemble(μ, σ, S0, tspan; seeds=nothing)
42+
T_ = typeof(σ)
43+
μ = convert(T_, μ)
44+
S0 = convert(T_, S0)
45+
t0, Tval = tspan
46+
tspan = (convert(T_, t0), convert(T_, Tval))
47+
48+
base_proc = GeometricBrownianMotionProcess(μ, σ, tspan[1], S0)
49+
base_prob = NoiseProblem(base_proc, tspan)
50+
51+
N = seeds === nothing ? 100 : length(seeds)
52+
seeds = seeds === nothing ? collect(1:N) : seeds
53+
modify = (prob, seed, i) -> remake(prob; seed=seed)
54+
return CustomEnsembleProblem(base_prob, collect(seeds), modify)
55+
end
56+
57+
# --------------------------
58+
# Monte Carlo Vega (FD)
59+
# --------------------------
60+
61+
function mc_black_scholes_vega(; σ=0.2, S0=100.0, K=100.0, r=0.01, T=1.0, ε=1e-4, N=10_000, dt=1/250)
62+
tspan = (0.0, T)
63+
μ = r # risk-neutral drift
64+
65+
base_ens = build_noise_gbm_ensemble(μ, σ, S0, tspan; seeds=1:N)
66+
bump_ens = build_noise_gbm_ensemble(μ, σ + ε, S0, tspan; seeds=1:N)
67+
68+
sol_base = solve_custom_ensemble(base_ens; dt=dt)
69+
sol_bump = solve_custom_ensemble(bump_ens; dt=dt)
70+
71+
payoffs_base = max.(last.(map(s -> s.u, sol_base.solutions)) .- K, 0.0)
72+
payoffs_bump = max.(last.(map(s -> s.u, sol_bump.solutions)) .- K, 0.0)
73+
74+
df = exp(-r * T)
75+
vega_est = mean(payoffs_bump .- payoffs_base) / ε * df
76+
77+
return vega_est
78+
end
79+
80+
# --------------------------
81+
# Monte Carlo Vega (AD)
82+
# --------------------------
83+
84+
function mc_black_scholes_vega_ad(; σ=0.2, S0=100.0, K=100.0, r=0.01, T=1.0, N=10_000, dt=1/250)
85+
seeds = 1:N
86+
payoff = x -> max(x - K, 0.0)
87+
df = exp(-r * T)
88+
89+
f = σ_val -> begin
90+
T_ = typeof(σ_val)
91+
μ_ = convert(T_, r) # risk-neutral drift
92+
S0_ = convert(T_, S0)
93+
tspan = (convert(T_, 0.0), convert(T_, T))
94+
95+
base_ens = build_noise_gbm_ensemble(μ_, σ_val, S0_, tspan; seeds=seeds)
96+
sol = solve_custom_ensemble(base_ens; dt=dt)
97+
final_vals = map(s -> s.u[end], sol.solutions)
98+
mean(payoff.(final_vals)) * df
99+
end
100+
101+
vega = ForwardDiff.derivative(f, σ)
102+
return vega
103+
end
104+
105+
# --------------------------
106+
# Analytic Black-Scholes Vega
107+
# --------------------------
108+
109+
function bs_vega_analytic(S, K, r, σ, T)
110+
d1 = (log(S / K) + (r + 0.5 * σ^2) * T) / (σ * sqrt(T))
111+
return S * sqrt(T) * pdf(Normal(), d1)
112+
end
113+
114+
# --------------------------
115+
# Run Comparison
116+
# --------------------------
117+
118+
S0 = 100.0
119+
K = 100.0
120+
r = 0.01
121+
T = 1.0
122+
σ = 0.2
123+
N = 100000
124+
125+
vega_mc = mc_black_scholes_vega(; σ=σ, S0=S0, K=K, r=r, T=T, N=N)
126+
vega_ad = mc_black_scholes_vega_ad(; σ=σ, S0=S0, K=K, r=r, T=T, N=N)
127+
vega_an = bs_vega_analytic(S0, K, r, σ, T)
128+
129+
println("Monte Carlo Vega (FD): ", vega_mc)
130+
println("Monte Carlo Vega (AD): ", vega_ad)
131+
println("Analytic Vega: ", vega_an)
132+
println("Rel Error (FD): ", abs(vega_mc - vega_an) / vega_an)
133+
println("Rel Error (AD): ", abs(vega_ad - vega_an) / vega_an)

‎examples/mc_pathwise_derivative.jl

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
using DifferentialEquations
2+
using DiffEqNoiseProcess
3+
using ForwardDiff
4+
using Plots
5+
6+
# Parameters
7+
μ = 0.05
8+
σ = 0.2
9+
t0 = 0.0
10+
T = 1.0
11+
W0 = 100.0
12+
tspan = (t0, T)
13+
N = 250
14+
dt = (T - t0) / N
15+
seed = 1
16+
17+
# --- Terminal value of GBM with type stability
18+
function terminal_value(σ_val::Real)
19+
T_ = typeof(σ_val)
20+
μ_ = convert(T_, μ)
21+
W0_ = convert(T_, W0)
22+
t0_ = convert(T_, t0)
23+
24+
proc = GeometricBrownianMotionProcess(μ_, σ_val, t0_, W0_)
25+
prob = NoiseProblem(proc, tspan; seed=seed)
26+
sol = DifferentialEquations.solve(prob, EM(); dt=dt)
27+
return sol.u[end]
28+
end
29+
30+
# --- Finite difference
31+
function finite_difference(f, x; ε=1e-5)
32+
(f(x + ε) - f(x - ε)) / (2ε)
33+
end
34+
35+
# --- Compute values
36+
val = terminal_value(σ)
37+
grad_ad = ForwardDiff.derivative(terminal_value, σ)
38+
grad_fd = finite_difference(terminal_value, σ)
39+
40+
println("Terminal Value at σ = $σ: ", val)
41+
println("∂Value/∂σ (AD): ", grad_ad)
42+
println("∂Value/∂σ (FD): ", grad_fd)
43+
44+
# --- Plot
45+
σ_bumped = σ + 0.05
46+
sol = DifferentialEquations.solve(NoiseProblem(GeometricBrownianMotionProcess(μ, σ, t0, W0), tspan; seed=seed), dt=dt)
47+
sol_bumped = DifferentialEquations.solve(NoiseProblem(GeometricBrownianMotionProcess(μ, σ_bumped, t0, W0), tspan; seed=seed), dt=dt)
48+
49+
plot(sol.t, sol.u, label="σ = $σ", linewidth=2)
50+
plot!(sol_bumped.t, sol_bumped.u, label="σ = $σ_bumped", linestyle=:dash, linewidth=2)
51+
xlabel!("Time")
52+
ylabel!("GBM Value")
53+
title!("GBM: Vol Bump and Pathwise Sensitivity")

‎examples/mc_pathwise_derivative_discretized.jl

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
using DifferentialEquations
2+
using ForwardDiff
3+
using Plots
4+
5+
# Parameters
6+
μ = 0.05
7+
σ = 0.2
8+
t0 = 0.0
9+
T = 1.0
10+
W0 = 100.0
11+
tspan = (t0, T)
12+
N = 250
13+
dt = (T - t0) / N
14+
seed = 1
15+
16+
# --- Drift and diffusion functions
17+
function drift!(du, u, p, t)
18+
du[1] = p[1] * u[1] # p[1] = μ
19+
end
20+
21+
function diffusion!(du, u, p, t)
22+
du[1] = p[2] * u[1] # p[2] = σ
23+
end
24+
25+
# --- Terminal value of GBM using SDEProblem
26+
function terminal_value_sde(σ_val::Real)
27+
T_ = typeof(σ_val)
28+
μ_ = convert(T_, μ)
29+
W0_ = convert(T_, W0)
30+
31+
u0 = [W0_]
32+
p = [μ_, σ_val]
33+
34+
prob = SDEProblem(drift!, diffusion!, u0, tspan, p; seed=seed)
35+
sol = DifferentialEquations.solve(prob, EM(); dt=dt)
36+
return sol.u[end][1]
37+
end
38+
39+
# --- Finite difference
40+
function finite_difference(f, x; ε=1e-5)
41+
(f(x + ε) - f(x - ε)) / (2ε)
42+
end
43+
44+
# --- Compute values
45+
val = terminal_value_sde(σ)
46+
grad_ad = ForwardDiff.derivative(terminal_value_sde, σ)
47+
grad_fd = finite_difference(terminal_value_sde, σ)
48+
49+
println("Terminal Value at σ = $σ: ", val)
50+
println("∂Value/∂σ (AD): ", grad_ad)
51+
println("∂Value/∂σ (FD): ", grad_fd)
52+
53+
# --- Plot
54+
σ_bumped = σ + 0.05
55+
prob = SDEProblem(drift!, diffusion!, [W0], tspan, [μ, σ]; seed=seed)
56+
prob_bumped = SDEProblem(drift!, diffusion!, [W0], tspan, [μ, σ_bumped]; seed=seed)
57+
58+
sol = DifferentialEquations.solve(prob, EM(); dt=dt)
59+
sol_bumped = DifferentialEquations.solve(prob_bumped, EM(); dt=dt)
60+
61+
plot(sol.t, first.(sol.u), label="σ = $σ", linewidth=2)
62+
plot!(sol_bumped.t, first.(sol_bumped.u), label="σ = $σ_bumped", linestyle=:dash, linewidth=2)
63+
xlabel!("Time")
64+
ylabel!("GBM Value")
65+
title!("GBM (SDE): Vol Bump and Pathwise Sensitivity")

0 commit comments

Comments
 (0)