Skip to content

Commit f167047

Browse files
aleCombialeCombi
committed
Trying montecarlo derivatives with lenses
1 parent 593ccd0 commit f167047

File tree

2 files changed

+267
-0
lines changed

2 files changed

+267
-0
lines changed

‎examples/mc_black_scholes_lens_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+
using Accessors
7+
8+
# --------------------------
9+
# Custom Ensemble Framework
10+
# --------------------------
11+
12+
struct CustomEnsembleProblem{P, F}
13+
base_problem::P
14+
seeds::Vector{Int64}
15+
modify::F # (base_problem, seed, index) -> new_problem
16+
end
17+
18+
struct CustomEnsembleSolution{S}
19+
solutions::Vector{S}
20+
seeds::Vector{Int64}
21+
end
22+
23+
function solve_custom_ensemble(prob::CustomEnsembleProblem; dt, solver=EM(), trajectories=nothing)
24+
N = trajectories === nothing ? length(prob.seeds) : trajectories
25+
seeds = length(prob.seeds) == N ? prob.seeds : collect(1:N)
26+
27+
sols = Vector{Any}(undef, N)
28+
29+
Threads.@threads for i in 1:N
30+
seed = seeds[i]
31+
pmod = prob.modify(prob.base_problem, seed, i)
32+
sols[i] = DifferentialEquations.solve(pmod, solver; dt=dt)
33+
end
34+
35+
return CustomEnsembleSolution(sols, seeds)
36+
end
37+
38+
# --------------------------
39+
# GBM with NoiseProcess
40+
# --------------------------
41+
42+
function build_noise_gbm_ensemble(μ, σ, S0, tspan; seeds=nothing)
43+
T_ = typeof(σ)
44+
μ = convert(T_, μ)
45+
S0 = convert(T_, S0)
46+
t0, Tval = tspan
47+
tspan = (convert(T_, t0), convert(T_, Tval))
48+
49+
base_proc = GeometricBrownianMotionProcess(μ, σ, tspan[1], S0)
50+
base_prob = NoiseProblem(base_proc, tspan)
51+
52+
N = seeds === nothing ? 100 : length(seeds)
53+
seeds = seeds === nothing ? collect(1:N) : seeds
54+
modify = (prob, seed, i) -> remake(prob; seed=seed)
55+
return CustomEnsembleProblem(base_prob, collect(seeds), modify)
56+
end
57+
58+
# --------------------------
59+
# Param Wrapper & Lens
60+
# --------------------------
61+
62+
struct BSMonteCarloSetup{T}
63+
S0::T
64+
K::T
65+
r::T
66+
σ::T
67+
T::T
68+
N::Int
69+
dt::T
70+
end
71+
72+
struct EnsembleProblemLens end
73+
74+
function (lens::EnsembleProblemLens)(p::BSMonteCarloSetup)
75+
tspan = (zero(p.T), p.T)
76+
build_noise_gbm_ensemble(p.r, p.σ, p.S0, tspan; seeds=1:p.N)
77+
end
78+
79+
function Accessors.set(p::BSMonteCarloSetup, ::EnsembleProblemLens, new_σ)
80+
T = typeof(new_σ)
81+
return BSMonteCarloSetup{T}(T(p.S0), T(p.K), T(p.r), new_σ, T(p.T), p.N, T(p.dt))
82+
end
83+
84+
# --------------------------
85+
# Analytic Black-Scholes
86+
# --------------------------
87+
88+
function bs_call_price(S, K, r, σ, T)
89+
d1 = (log(S / K) + (r + 0.5 * σ^2) * T) / (σ * sqrt(T))
90+
d2 = d1 - σ * sqrt(T)
91+
return S * cdf(Normal(), d1) - K * exp(-r * T) * cdf(Normal(), d2)
92+
end
93+
94+
function bs_vega_analytic(S, K, r, σ, T)
95+
d1 = (log(S / K) + (r + 0.5 * σ^2) * T) / (σ * sqrt(T))
96+
return S * sqrt(T) * pdf(Normal(), d1)
97+
end
98+
99+
# --------------------------
100+
# AD Price & Vega via Lens
101+
# --------------------------
102+
103+
function price_from_params(p::BSMonteCarloSetup, lens::EnsembleProblemLens)
104+
prob = lens(p)
105+
sol = solve_custom_ensemble(prob; dt=p.dt)
106+
payoffs = map(s -> max(s.u[end] - p.K, 0.0), sol.solutions)
107+
return mean(payoffs) * exp(-p.r * p.T)
108+
end
109+
110+
function vega_from_ad(p::BSMonteCarloSetup, lens::EnsembleProblemLens)
111+
f = σ_val -> price_from_params(Accessors.set(p, lens, σ_val), lens)
112+
return ForwardDiff.derivative(f, p.σ)
113+
end
114+
115+
# --------------------------
116+
# Run Comparison
117+
# --------------------------
118+
119+
params = BSMonteCarloSetup(100.0, 100.0, 0.01, 0.2, 1.0, 10_000, 1/250)
120+
lens = EnsembleProblemLens()
121+
122+
price_mc = price_from_params(params, lens)
123+
vega_ad = vega_from_ad(params, lens)
124+
price_an = bs_call_price(params.S0, params.K, params.r, params.σ, params.T)
125+
vega_an = bs_vega_analytic(params.S0, params.K, params.r, params.σ, params.T)
126+
127+
println("Monte Carlo Price: ", price_mc)
128+
println("Analytic Price: ", price_an)
129+
println("Rel Error (Price): ", abs(price_mc - price_an) / price_an)
130+
131+
println("Monte Carlo Vega (AD): ", vega_ad)
132+
println("Analytic Vega: ", vega_an)
133+
println("Rel Error (Vega): ", abs(vega_ad - vega_an) / vega_an)

‎examples/mc_black_scholes_lens_wrapped.jl

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

0 commit comments

Comments
 (0)