Merged drif and diffusion functions into SDEFunction, combining SDEFunctions, also with analytical definitions, SDEProblems given a set of processes and a correlation matrix.

Author: aleCombi

examples/Project.toml

@@ -5,6 +5,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Hedgehog2 = "7f16798b-0e18-40de-98af-932948254698"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"

examples/correlated_gbms.jl

@@ -0,0 +1,25 @@
+using Revise, Hedgehog2, DifferentialEquations, Plots
+
+process = GBMProcess(0.05, 0.2)  # μ = 5%, σ = 20%
+sde_function_1 = get_sde_function(process)
+sde_function_2 = get_sde_function(process)
+combined_sde_function = combine_sde_functions((sde_function_1, sde_function_2))
+
+ρ = -0.8
+# Correlation matrix
+Γ = [1.0 ρ; 
+     ρ 1.0]
+
+u0 = [1.0, 1.0]  # Initial conditions
+tspan = (0.0, 2.0)  # Time span
+
+noise = CorrelatedWienerProcess(Γ, tspan[1], [0.0, 0.0])
+# Create multi-dimensional SDE problem
+prob = SDEProblem(combined_sde_function, u0, tspan, "p", noise=noise)
+prob2 = get_sde_problem((process, process), Γ, u0, tspan, "p")
+
+sol = solve(prob)
+plot(sol, plot_analytic=true)
+
+sol2 = solve(prob2)
+plot!(sol, plot_analytic=true)

examples/time_dependend_correlated_gbms.jl

@@ -0,0 +1,40 @@
+using Revise, Hedgehog2, DifferentialEquations, Plots
+
+# Define a constant volatility function and its integral
+σ_const = t -> 0.2  # Constant volatility (σ_0 = 0.2)
+function mean_σ²_const(t)
+    σ1² = 0.2^2  # σ_1^2
+    σ2² = 0.4^2  # σ_2^2
+
+    if t < 0.5
+        return σ1²
+    else
+        return (σ1² * 0.5 + σ2² * (t - 0.5)) / t
+    end
+end
+
+
+# Create the GBMTimeDependent process with constant volatility
+process_const = GBMTimeDependent(
+    0.05,          # Drift
+    σ_const,       # Constant volatility function
+    mean_σ²_const  # Integral of σ²(t)
+)
+
+# Get SDE function
+sde_func_const = get_sde_function(process_const)
+
+# Initial value and time span
+u0 = [1.0]
+tspan = (0.0, 2.0)
+
+
+sde_problem_td = SDEProblem(sde_func_const, u0, tspan, seed=12)
+sol_td = solve(sde_problem_td)
+plot(sol_td)
+
+process_gbm = GBMProcess(0.05, 0.2)
+sde_func = get_sde_function(process_gbm)
+sde_problem = SDEProblem(sde_func, u0, tspan, seed=12)
+sol = solve(sde_problem)
+plot!(sol)

examples/time_dependent_gbm.jl

@@ -1,6 +1,41 @@
 using DifferentialEquations
 
-export AbstractStochasticProcess, GBMProcess, HestonProcess, get_drift, get_diffusion, get_analytic_solution, get_sde_function
+export AbstractStochasticProcess, GBMProcess, HestonProcess, get_sde_function, combine_sde_functions, GBMTimeDependent, get_sde_problem
+
+function combine_sde_functions(sde_functions::Tuple{Vararg{SDEFunction}})
+    n = length(sde_functions)  # Number of 1D SDEs
+
+    # Construct multi-dimensional drift function
+    f = (u, p, t) -> begin
+        return map(i -> sde_functions[i].f([u[i]], p, t)[1], 1:n)
+    end
+
+    # Construct multi-dimensional diffusion function
+    g = (u, p, t) -> begin
+        return map(i -> sde_functions[i].g([u[i]], p, t)[1], 1:n)
+    end
+
+    # Check if all SDEFunctions have an analytic solution
+    if all(sf -> sf.analytic !== nothing, sde_functions)
+        # Construct a multi-dimensional analytical solution
+        analytic = (u0, p, t, W) -> begin
+            return map(i -> sde_functions[i].analytic([u0[i]], p, t, W[i])[1], 1:n)
+        end
+    else
+        analytic = nothing  # No combined analytic solution if any function lacks it
+    end
+
+    return SDEFunction(f, g, analytic=analytic)
+end
+
+function get_sde_problem(processes, Γ, u0, tspan, p; kwargs...)
+    sde_functions = get_sde_function.(processes)
+    combined_sde_function = combine_sde_functions(sde_functions)
+    noise = CorrelatedWienerProcess(Γ, tspan[1], zeros(length(processes)))
+
+    return SDEProblem(combined_sde_function, u0, tspan, p; noise=noise, kwargs...)
+end
+
 
 # Abstract type for stochastic processes
 """
@@ -21,6 +56,16 @@ struct GBMProcess <: AbstractStochasticProcess
     σ
 end
 
+function get_sde_function(process::GBMProcess)
+    drift(u, p, t) = process.μ .* u 
+
+    diffusion(u, p, t) = process.σ .* u
+
+    analytic_solution(u₀, p, t, W) = u₀ .* exp.((process.μ - (process.σ^2) / 2) .* t .+ process.σ * W)
+
+    return SDEFunction(drift, diffusion, analytic = analytic_solution)
+end
+
 # Define the Heston process
 """
     HestonProcess <: AbstractStochasticProcess
@@ -40,84 +85,44 @@ 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)
+function get_sde_function(process::HestonProcess)
+    function drift(u, p, t)
+        S, V = u
+        return [process.μ * S, process.κ * (process.θ - V)] 
+    end
+    
+    function diffusion(u, p, t)
+        S, V = u
+        return [process.σ * S * sqrt(V), process.σ * sqrt(V)]  
+    end
+
+    return SDEFunction(drift, diffusion)
 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)
+struct GBMTimeDependent
+    μ
+    σ
+    mean_σ²
 end
 
-"""
-    get_diffusion!(process::AbstractStochasticProcess)
+function get_sde_function(process::GBMTimeDependent)
+    drift(u, p, t) = process.μ .* u
 
-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
+    diffusion(u, p, t) = sqrt(process.mean_σ²(t)) .* u
 
-function get_analytic_solution(process::P) where P <: AbstractStochasticProcess
-    return Nothing()
-end
+    # observe that the noise here is integrated.
+    # That is, W is not such that dX = mu dt + sigma_t dW_t but it is W' such that W'_t = t (int_0^t sigma_s dW_t) / (int_0^t sigma_s^2)
+    function analytic_solution(u₀, p, t, W)
+        I_t = sqrt(process.mean_σ²(t)) ./ t ./ t
+        return u₀ .* exp.((process.μ .- 0.5 * I_t^2) .* t .+ I_t .* W)
+    end
 
-function get_analytic_solution(::GBMProcess)
-    return (u₀, p, t, W) -> u₀ * exp((0.05 - (0.2^2) / 2) * t + 0.2 * W)
+    return SDEFunction(drift, diffusion, analytic = analytic_solution)
 end
 
-function get_sde_function(process::P) where P<:AbstractStochasticProcess
-    return SDEFunction(get_drift(process), get_diffusion(process), analytic=get_analytic_solution(process))
+# dX_t = (θ_t - a t) dt + σ_t dW_t
+struct OUTimeDependent
+    θ
+    a
+    σ
 end
-