Added initial SDE DifferentialEquations.jl compatibility to be used in the future with Montecarlo and FFT methods, plus an example where Heston and BS are plotted.

Author: aleCombi

Project.toml

@@ -7,13 +7,15 @@ version = "1.0.0-DEV"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 
 [compat]
 Accessors = "0.1.42"
 BenchmarkTools = "1.6.0"
 Dates = "1.11.0"
+DifferentialEquations = "7.16.0"
 Distributions = "0.25.117"
 ForwardDiff = "0.10.38"
 julia = "1.11"

examples/Project.toml

@@ -2,8 +2,10 @@
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Hedgehog2 = "7f16798b-0e18-40de-98af-932948254698"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

examples/gbm_simulation.jl

@@ -0,0 +1,52 @@
+
+# Simulate GBM process
+"""
+    simulate_gbm()
+
+Simulates a Geometric Brownian Motion (GBM) process using the Euler-Maruyama method.
+Returns the numerical solution.
+"""
+function simulate_gbm()
+    u0 = 100.0  # Initial stock price
+    tspan = (0.0, 1.0)  # Simulation for 1 year
+    process = GBMProcess(0.05, 0.2)  # μ = 5%, σ = 20%
+    sde_function = get_sde_function(process)
+    prob = SDEProblem(sde_function, u0, tspan)
+    sol = solve(prob, EM(), dt=0.01)  # Euler-Maruyama solver
+
+    return sol
+end
+
+# Simulate Heston process
+"""
+    simulate_heston()
+
+Simulates the Heston stochastic volatility model using the Euler-Maruyama method.
+Returns the numerical solution.
+"""
+function simulate_heston()
+    u0 = [100.0, 0.04]  # Initial (Stock Price, Variance)
+    tspan = (0.0, 1.0)  # Simulation for 1 year
+    process = HestonProcess(0.05, 2.0, 0.04, 0.3, -0.5)  # Heston parameters
+    prob = SDEProblem(get_sde_function(process), u0, tspan, process)
+    sol = solve(prob, EM(), dt=0.01)  # Euler-Maruyama solver
+
+    return sol
+end
+
+# Running the simulations
+using Revise, Plots, Hedgehog2, DifferentialEquations
+
+sol_gbm = simulate_gbm()
+sol_heston = simulate_heston()
+
+# Plot GBM
+gbm_plot = plot(sol_gbm, plot_analytic = true, label="GBM Path", xlabel="Time", ylabel="Stock Price", lw=2)
+
+# Plot Heston
+heston_plot = plot!(sol_heston.t, [u[1] for u in sol_heston.u], label="Heston Stock Price", xlabel="Time", ylabel="Stock Price", lw=2)
+heston_plot_var = plot!(sol_heston.t, [u[2] for u in sol_heston.u], label="Variance", lw=2)
+
+display(gbm_plot)
+display(heston_plot)
+display(heston_plot_var)

src/Hedgehog2.jl

@@ -12,5 +12,6 @@ include("pricing_methods/pricing_methods.jl")
 include("pricing_methods/black_scholes.jl")
 include("pricing_methods/cox_ross_rubinstein.jl")
 include("sensitivity_methods/delta_methods.jl")
+include("stochastic_processes/stochastic_processes.jl")
 
 end

src/stochastic_processes/stochastic_processes.jl

@@ -0,0 +1,123 @@
+using DifferentialEquations
+
+export AbstractStochasticProcess, GBMProcess, HestonProcess, get_drift, get_diffusion, get_analytic_solution, get_sde_function
+
+# Abstract type for stochastic processes
+"""
+    AbstractStochasticProcess
+
+An abstract type representing a generic stochastic process.
+"""
+abstract type AbstractStochasticProcess end
+
+# Define the GBM process
+"""
+    GBMProcess <: AbstractStochasticProcess
+
+Represents a Geometric Brownian Motion (GBM) process with drift `μ` and volatility `σ`.
+"""
+struct GBMProcess <: AbstractStochasticProcess
+    μ
+    σ
+end
+
+# Define the Heston process
+"""
+    HestonProcess <: AbstractStochasticProcess
+
+Represents the Heston stochastic volatility model with parameters:
+- `μ`: Drift of the asset price
+- `κ`: Mean reversion speed of variance
+- `θ`: Long-run variance
+- `σ`: Volatility of variance
+- `ρ`: Correlation between asset and variance processes
+"""
+struct HestonProcess <: AbstractStochasticProcess
+    μ
+    κ
+    θ
+    σ
+    ρ
+end
+
+# Drift function for GBM
+"""
+    drift(process::GBMProcess, u, t)
+
+Computes the drift term of the GBM process at time `t` for state `u`.
+Drift equation: `du/dt = μ * u`
+"""
+function drift(process::GBMProcess, u, t)
+    return process.μ .* u  # Element-wise broadcasting for array compatibility
+end
+
+# Drift function for Heston
+"""
+    drift(process::HestonProcess, u, t)
+
+Computes the drift term of the Heston process at time `t` for state `u = (S, V)`.
+Drift equations:
+- `dS/dt = μ * S`
+- `dV/dt = κ * (θ - V)`
+"""
+function drift(process::HestonProcess, u, t)
+    S, V = u
+    return [process.μ * S, process.κ * (process.θ - V)]  # Drift for (Stock price, Variance)
+end
+
+# Diffusion function for GBM
+"""
+    diffusion(process::GBMProcess, u, t)
+
+Computes the diffusion term of the GBM process at time `t` for state `u`.
+Diffusion equation: `dW_t = σ * u * dB_t`
+"""
+function diffusion(process::GBMProcess, u, t)
+    return process.σ .* u  # Element-wise broadcasting for array compatibility
+end
+
+# Diffusion function for Heston
+"""
+    diffusion(process::HestonProcess, u, t)
+
+Computes the diffusion term of the Heston process at time `t` for state `u = (S, V)`.
+Diffusion equations:
+- `dS_t = σ * sqrt(V) * dB_t`
+- `dV_t = ξ * sqrt(V) * dW_t`, with correlation `ρ`.
+"""
+function diffusion(process::HestonProcess, u, t)
+    S, V = u
+    return [process.σ * S * sqrt(V), process.σ * sqrt(V)]  # Diffusion for (Stock price, Variance)
+end
+
+# Higher-order functions to return drift! and diffusion!
+"""
+    get_drift!(process::AbstractStochasticProcess)
+
+Returns the in-place drift function `drift!` for the given process.
+"""
+function get_drift(process::P) where P<:AbstractStochasticProcess
+    return (u, p, t) -> drift(process, u, t)
+end
+
+"""
+    get_diffusion!(process::AbstractStochasticProcess)
+
+Returns the in-place diffusion function `diffusion!` for the given process.
+"""
+function get_diffusion(process::P) where P<:AbstractStochasticProcess
+    return (u, p, t) -> diffusion(process, u, t)
+end
+
+function get_analytic_solution(process::P) where P <: AbstractStochasticProcess
+    return Nothing()
+end
+
+function get_analytic_solution(::GBMProcess)
+    return (u₀, p, t, W) -> u₀ * exp((0.05 - (0.2^2) / 2) * t + 0.2 * W)
+end
+
+function get_sde_function(process::P) where P<:AbstractStochasticProcess
+    return SDEFunction(get_drift(process), get_diffusion(process), analytic=get_analytic_solution(process))
+end
+