Fast special cases for MonteCarlo

[JADekker]

In the current implementation, I believe every MonteCarlo problem gets mapped to a NoiseProblem/EnsembleProblem. While the benefit of this framework is that it works in great generality, it is much slower and requires much more allocations than a more direct approach for problems that admit this. The following is 300x faster on my device and cuts allocations from 800.000ish to 18 in a simple example.

If such special-cases are desirable, what would be a good infrastructure (what I did now is just the first thing that came to mind); and would it be good to have a tracking issue to identify (and solve) such special cases?

[JADekker]

using Revise
using Hedgehog
using Dates
using Printf
using BenchmarkTools
using Random 
struct BSExactProblem{T,U,V,W,X}
    r::T
    σ::U
    t₀::V
    S₀::W
    T::X
end

struct BSExactProblemSimulated{T}
    S::T
end

function sde_problem2(
    problem::PricingProblem{Payoff, Inputs},
    ::LognormalDynamics,
    ::BlackScholesExact,
    ) where {Payoff, Inputs <: BlackScholesInputs}
    market = problem.market_inputs
    T = yearfrac(market.referenceDate, problem.payoff.expiry)

    r = zero_rate(market.rate, 0.0)
    σ = get_vol(market.sigma, nothing, nothing)
    S₀ = market.spot
    t₀ = zero(S₀)

    r, σ, t₀, S₀ = promote(r, σ, t₀, S₀)

    return BSExactProblem(r, σ, t₀, S₀, T)
end

function simulate_paths(
    sde_prob::BSExactProblem,
    method::MonteCarlo,
    ::Hedgehog.NoVarianceReduction
)
    Δt = sde_prob.T - sde_prob.t₀
    σ = sde_prob.σ
    S₀ = sde_prob.S₀
    r = sde_prob.r
    N = method.config.trajectories

    drift_part = ((r - 0.5 * σ^2) * Δt)
    diffusion_part = σ * sqrt(Δt)
    S_out = S₀ .* exp.(drift_part .+ diffusion_part .* randn(N))
    return BSExactProblemSimulated(S_out)
end

function simulate_paths(
    sde_prob::BSExactProblem,
    method::MonteCarlo,
    ::Hedgehog.Antithetic
)
    Δt = sde_prob.T - sde_prob.t₀
    σ = sde_prob.σ
    S₀ = sde_prob.S₀
    r = sde_prob.r
    N = method.config.trajectories

    drift_part = ((r - 0.5 * σ^2) * Δt)
    diffusion_part = σ * sqrt(Δt)

    # We ensure that normal_sol is of the right type
    normal_sol = Vector{typeof(drift_part*diffusion_part)}(undef, N)
    randn!(normal_sol)

    antithetic_sol = S₀ .* exp.(drift_part .- diffusion_part .* normal_sol)
    normal_sol .= S₀ .* exp.(drift_part .+ diffusion_part .* normal_sol)
    return (BSExactProblemSimulated(normal_sol), BSExactProblemSimulated(antithetic_sol))
end

function reduce_payoffs(
    result::Tuple{BSExactProblemSimulated, BSExactProblemSimulated},
    payoff::F,
    ::Hedgehog.Antithetic,
    dynamics::Hedgehog.PriceDynamics,
    strategy::Hedgehog.SimulationStrategy,
) where {F}
    S1, S2 = result
    return (payoff.(S1.S) + payoff.(S2.S))/2
end

function solve2(
    prob::PricingProblem{VanillaOption{TS, TE, European, C, Spot}, I},
    method::MonteCarlo{D, S},
) where {TS, TE, C, I, D, S}
    strategy = method.strategy
    dynamics = method.dynamics
    config = method.config

    sde_prob = sde_problem2(prob, method.dynamics, method.strategy)
    ens = simulate_paths(sde_prob, method, config.variance_reduction)
    payoffs = reduce_payoffs(ens, prob.payoff, config.variance_reduction, dynamics, strategy)
    discount = df(prob.market_inputs.rate, prob.payoff.expiry)
    price = discount * mean(payoffs)

    return Hedgehog.MonteCarloSolution(prob, method, price, ens)
end

function test()
    spot = 100.0
    strike = 100.0
    rate = 0.05
    sigma = 0.20
    reference_date = Date(2023, 1, 1)
    expiry = reference_date + Year(1)

    # Create the payoff (European call option)
    payoff = VanillaOption(strike, expiry, European(), Call(), Spot())

    # Create market inputs
    market_inputs = BlackScholesInputs(reference_date, rate, spot, sigma)

    # Create pricing problem
    prob = PricingProblem(payoff, market_inputs)

    trajectories = 5_000
    mc_exact_method =
        MonteCarlo(LognormalDynamics(), BlackScholesExact(), SimulationConfig(trajectories))
    mc_exact_solution = solve(prob, mc_exact_method)
    display(mc_exact_solution.price)
    mc_exact_solution = solve2(prob, mc_exact_method)
    display(mc_exact_solution.price)
    display(@benchmark solve($prob, $mc_exact_method))
    display(@benchmark solve2($prob, $mc_exact_method))
end

test()

gives

10.44471025405915
10.230865999334211
BenchmarkTools.Trial: 272 samples with 1 evaluation per sample.
 Range (min … max):  15.796 ms … 26.483 ms  ┊ GC (min … max):  0.00% … 15.83%
 Time  (median):     18.651 ms              ┊ GC (median):    13.39%
 Time  (mean ± σ):   18.368 ms ±  2.243 ms  ┊ GC (mean ± σ):  10.25% ±  9.00%

  ▆█▄             ▃▃▂▄▂▂▂▁▃▁▂ ▂                                
  █████▁▄▄▁▁▁▁▁▁▆▆███████████████▇▄▄▆▄▆▁▆▆▄▆▁▄▁▆▄▁▁▁▄▁▁▁▄▁▄▁▄ ▆
  15.8 ms      Histogram: log(frequency) by time      25.2 ms <

 Memory estimate: 39.72 MiB, allocs estimate: 890149.
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  43.834 μs …   7.981 ms  ┊ GC (min … max):  0.00% … 99.06%
 Time  (median):     57.438 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   73.604 μs ± 166.850 μs  ┊ GC (mean ± σ):  19.74% ± 11.18%

  ▅█▅▃▁                                                        ▁
  ██████▇▆▅▅▅▄▄▅▅▅▅▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▄▅▄▅▁▄▅▃▄▅▅▅▅ █
  43.8 μs       Histogram: log(frequency) by time       598 μs <

 Memory estimate: 384.38 KiB, allocs estimate: 18.

on julia 1.11.5

[JADekker]

See #13 for a PR with an example implementation!

[aleCombi]

My understanding is that some copying of the noise process is the cause of this slowness of solve(NoiseProblem, ...). I checked your code against NoiseProblem code with ProfileView to see this.

I want to open an issue on https://github.com/SciML/DiffEqNoiseProcess.jl as well.

One issue though is that the machinery of DifferentialEquations.jl is used to reuse seeds (for pathwise montecarlo greeks and for antithetic variates).
Perhaps the cleanest way is to define a dispatch of StochasticDiffEq.solve(ensemble_prob, EM(); dt = dt, trajectories=config.trajectories), something like solve(ensemble_prob, HedgehogExact(); dt = dt, trajectories=config.trajectories).
If we can return an EnsembleSolution this way, we might be able to get the performance benefits without losing the infrastructure.
The simulation scheme (EM(), HedgehogExact()) could be another field of the Montecarlo struct (of part of SimulationConfig).

[JADekker]

My understanding is that some copying of the noise process is the cause of this slowness of solve(NoiseProblem, ...). I checked your code against NoiseProblem code with ProfileView to see this.

I want to open an issue on https://github.com/SciML/DiffEqNoiseProcess.jl as well.

One issue though is that the machinery of DifferentialEquations.jl is used to reuse seeds (for pathwise montecarlo greeks and for antithetic variates). Perhaps the cleanest way is to define a dispatch of StochasticDiffEq.solve(ensemble_prob, EM(); dt = dt, trajectories=config.trajectories), something like solve(ensemble_prob, HedgehogExact(); dt = dt, trajectories=config.trajectories). If we can return an EnsembleSolution this way, we might be able to get the performance benefits without losing the infrastructure. The simulation scheme (EM(), HedgehogExact()) could be another field of the Montecarlo struct (of part of SimulationConfig).

Alternatively, you could define an own layer of infrastructure, which both EM and HedgehogExact would be tasked to work with: for European MC pricing for example, it is only necessary to store the final values, the rest of the things tracked by EnsembleSolution are not necessary there (and a similar case can be made for other applications).

[JADekker]

One issue though is that the machinery of DifferentialEquations.jl is used to reuse seeds (for pathwise montecarlo greeks and for antithetic variates).

I think my implementation in #13 also reuses the seed and handles antithetic variates. It should not be too hard to allow HedgeHogExact to handle such features analogously.

[aleCombi]

I tried to get the best of both worlds in #17. Still kind of a draft though.