diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index d22e82d..6253546 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -14,6 +14,9 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v2 + with: + version: '1' - uses: fredrikekre/runic-action@v1 with: version: '1' diff --git a/docs/pages.jl b/docs/pages.jl index 282b55b..e9a398b 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,3 +1,3 @@ # Put in a separate page so it can be used by SciMLDocs.jl -pages = ["Home" => "index.md", "diffeqfinancial.md"] +pages = ["Home" => "index.md", "Analytical Solutions" => "analytics.md", "API Reference" => "diffeqfinancial.md"] diff --git a/docs/src/analytics.md b/docs/src/analytics.md new file mode 100644 index 0000000..a27ecb3 --- /dev/null +++ b/docs/src/analytics.md @@ -0,0 +1,62 @@ +# Analytical Solutions + +DiffEqFinancial.jl provides analytical (closed-form) solutions for the mean, +variance, and standard deviation of the supported financial stochastic processes. +These functions provide exact results without requiring Monte Carlo simulation, +offering significant performance and accuracy benefits. + +## Geometric Brownian Motion + +The GBM process: ``dX = \mu X \, dt + \sigma X \, dW_t`` + +```@docs +gbm_mean +gbm_variance +gbm_std +``` + +## Ornstein-Uhlenbeck Process + +The OU process: ``dX = a(r - X) \, dt + \sigma \, dW_t`` + +```@docs +ou_mean +ou_variance +ou_std +ou_stationary_mean +ou_stationary_variance +``` + +## Cox-Ingersoll-Ross (CIR) Process + +The CIR process: ``dr = \kappa(\theta - r) \, dt + \sigma \sqrt{r} \, dW_t`` + +```@docs +cir_mean +cir_variance +cir_std +cir_stationary_mean +cir_stationary_variance +``` + +## Black-Scholes Log-Price + +The log-price process: ``d \ln S = (r - \sigma^2/2) \, dt + \sigma \, dW_t`` + +```@docs +bs_log_mean +bs_log_variance +bs_log_std +``` + +## Heston Model + +The Heston stochastic volatility model: +``dS = \mu S \, dt + \sqrt{v} S \, dW_1``, +``dv = \kappa(\Theta - v) \, dt + \sigma \sqrt{v} \, dW_2`` + +```@docs +heston_mean +heston_variance_mean +heston_variance_variance +``` diff --git a/src/DiffEqFinancial.jl b/src/DiffEqFinancial.jl index 0493d1d..d0eb4d9 100644 --- a/src/DiffEqFinancial.jl +++ b/src/DiffEqFinancial.jl @@ -8,6 +8,7 @@ using Markdown: @doc_str import RandomNumbers: Xorshifts include("problems.jl") +include("analytics.jl") export HestonProblem, BlackScholesProblem, GeneralizedBlackScholesProblem, ExtendedOrnsteinUhlenbeckProblem, OrnsteinUhlenbeckProblem, @@ -16,4 +17,11 @@ export HestonProblem, BlackScholesProblem, GeneralizedBlackScholesProblem, CIRProblem, CIRNoise +# Analytical solutions and moment functions +export gbm_mean, gbm_variance, gbm_std, + ou_mean, ou_variance, ou_std, ou_stationary_mean, ou_stationary_variance, + cir_mean, cir_variance, cir_std, cir_stationary_mean, cir_stationary_variance, + bs_log_mean, bs_log_variance, bs_log_std, + heston_mean, heston_variance_mean, heston_variance_variance + end # module diff --git a/src/analytics.jl b/src/analytics.jl new file mode 100644 index 0000000..d3c3c83 --- /dev/null +++ b/src/analytics.jl @@ -0,0 +1,378 @@ +# Analytical solutions, mean, and variance calculations for financial models + +# ============================================================================= +# Geometric Brownian Motion: dx = μx dt + σx dW_t +# ============================================================================= + +@doc doc""" +Analytical mean of Geometric Brownian Motion. + +```math +E[X(t)] = x_0 e^{\mu t} +``` + +### Arguments +- `μ`: Drift coefficient +- `u0`: Initial value ``x_0`` +- `t`: Time (scalar or array) + +### Returns +The expected value at time `t`. +""" +function gbm_mean(μ, u0, t) + return u0 * exp(μ * t) +end + +@doc doc""" +Analytical variance of Geometric Brownian Motion. + +```math +\text{Var}[X(t)] = x_0^2 e^{2\mu t} (e^{\sigma^2 t} - 1) +``` + +### Arguments +- `μ`: Drift coefficient +- `σ`: Volatility coefficient +- `u0`: Initial value ``x_0`` +- `t`: Time (scalar or array) + +### Returns +The variance at time `t`. +""" +function gbm_variance(μ, σ, u0, t) + return u0^2 * exp(2μ * t) * expm1(σ^2 * t) +end + +@doc doc""" +Analytical standard deviation of Geometric Brownian Motion. + +### Arguments +- `μ`: Drift coefficient +- `σ`: Volatility coefficient +- `u0`: Initial value ``x_0`` +- `t`: Time (scalar or array) + +### Returns +The standard deviation at time `t`. +""" +function gbm_std(μ, σ, u0, t) + return sqrt(gbm_variance(μ, σ, u0, t)) +end + +# ============================================================================= +# Ornstein-Uhlenbeck Process: dx = a(r - x)dt + σ dW_t +# ============================================================================= + +@doc doc""" +Analytical mean of the Ornstein-Uhlenbeck process. + +```math +E[X(t)] = r + (x_0 - r) e^{-at} +``` + +### Arguments +- `a`: Mean reversion speed +- `r`: Long-run mean level +- `u0`: Initial value ``x_0`` +- `t`: Time (scalar or array) + +### Returns +The expected value at time `t`. +""" +function ou_mean(a, r, u0, t) + return r + (u0 - r) * exp(-a * t) +end + +@doc doc""" +Analytical variance of the Ornstein-Uhlenbeck process. + +```math +\text{Var}[X(t)] = \frac{\sigma^2}{2a} (1 - e^{-2at}) +``` + +### Arguments +- `a`: Mean reversion speed +- `σ`: Volatility coefficient +- `t`: Time (scalar or array) + +### Returns +The variance at time `t`. +""" +function ou_variance(a, σ, t) + return (σ^2 / (2a)) * (-expm1(-2a * t)) +end + +@doc doc""" +Analytical standard deviation of the Ornstein-Uhlenbeck process. + +### Arguments +- `a`: Mean reversion speed +- `σ`: Volatility coefficient +- `t`: Time (scalar or array) + +### Returns +The standard deviation at time `t`. +""" +function ou_std(a, σ, t) + return sqrt(ou_variance(a, σ, t)) +end + +@doc doc""" +Stationary (long-run) mean of the Ornstein-Uhlenbeck process. + +```math +\lim_{t \to \infty} E[X(t)] = r +``` + +### Arguments +- `r`: Long-run mean level + +### Returns +The stationary mean. +""" +function ou_stationary_mean(r) + return r +end + +@doc doc""" +Stationary (long-run) variance of the Ornstein-Uhlenbeck process. + +```math +\lim_{t \to \infty} \text{Var}[X(t)] = \frac{\sigma^2}{2a} +``` + +### Arguments +- `a`: Mean reversion speed +- `σ`: Volatility coefficient + +### Returns +The stationary variance. +""" +function ou_stationary_variance(a, σ) + return σ^2 / (2a) +end + +# ============================================================================= +# Cox-Ingersoll-Ross (CIR) Process: dr = κ(θ - r)dt + σ√r dW_t +# ============================================================================= + +@doc doc""" +Analytical mean of the CIR process. + +```math +E[r(t)] = \theta + (r_0 - \theta) e^{-\kappa t} +``` + +### Arguments +- `κ`: Mean reversion speed +- `θ`: Long-run mean level +- `u0`: Initial value ``r_0`` +- `t`: Time (scalar or array) + +### Returns +The expected value at time `t`. +""" +function cir_mean(κ, θ, u0, t) + return θ + (u0 - θ) * exp(-κ * t) +end + +@doc doc""" +Analytical variance of the CIR process. + +```math +\text{Var}[r(t)] = r_0 \frac{\sigma^2}{\kappa} (e^{-\kappa t} - e^{-2\kappa t}) + \frac{\theta \sigma^2}{2\kappa} (1 - e^{-\kappa t})^2 +``` + +### Arguments +- `κ`: Mean reversion speed +- `θ`: Long-run mean level +- `σ`: Volatility coefficient +- `u0`: Initial value ``r_0`` +- `t`: Time (scalar or array) + +### Returns +The variance at time `t`. +""" +function cir_variance(κ, θ, σ, u0, t) + exp_κt = exp(-κ * t) + exp_2κt = exp(-2κ * t) + term1 = u0 * (σ^2 / κ) * (exp_κt - exp_2κt) + term2 = (θ * σ^2 / (2κ)) * (1 - exp_κt)^2 + return term1 + term2 +end + +@doc doc""" +Analytical standard deviation of the CIR process. + +### Arguments +- `κ`: Mean reversion speed +- `θ`: Long-run mean level +- `σ`: Volatility coefficient +- `u0`: Initial value ``r_0`` +- `t`: Time (scalar or array) + +### Returns +The standard deviation at time `t`. +""" +function cir_std(κ, θ, σ, u0, t) + return sqrt(cir_variance(κ, θ, σ, u0, t)) +end + +@doc doc""" +Stationary (long-run) mean of the CIR process. + +```math +\lim_{t \to \infty} E[r(t)] = \theta +``` + +### Arguments +- `θ`: Long-run mean level + +### Returns +The stationary mean. +""" +function cir_stationary_mean(θ) + return θ +end + +@doc doc""" +Stationary (long-run) variance of the CIR process. + +```math +\lim_{t \to \infty} \text{Var}[r(t)] = \frac{\theta \sigma^2}{2\kappa} +``` + +### Arguments +- `κ`: Mean reversion speed +- `θ`: Long-run mean level +- `σ`: Volatility coefficient + +### Returns +The stationary variance. +""" +function cir_stationary_variance(κ, θ, σ) + return θ * σ^2 / (2κ) +end + +# ============================================================================= +# Black-Scholes Log-Price: d ln S = (r - σ²/2)dt + σ dW_t +# ============================================================================= + +@doc doc""" +Analytical mean of the Black-Scholes log-price process. + +```math +E[\ln S(t)] = \ln S_0 + (r - \frac{\sigma^2}{2}) t +``` + +### Arguments +- `r`: Risk-free rate (constant) +- `σ`: Volatility coefficient +- `u0`: Initial log-price ``\ln S_0`` +- `t`: Time (scalar or array) + +### Returns +The expected log-price at time `t`. +""" +function bs_log_mean(r, σ, u0, t) + return u0 + (r - σ^2 / 2) * t +end + +@doc doc""" +Analytical variance of the Black-Scholes log-price process. + +```math +\text{Var}[\ln S(t)] = \sigma^2 t +``` + +### Arguments +- `σ`: Volatility coefficient +- `t`: Time (scalar or array) + +### Returns +The variance of log-price at time `t`. +""" +function bs_log_variance(σ, t) + return σ^2 * t +end + +@doc doc""" +Analytical standard deviation of the Black-Scholes log-price process. + +### Arguments +- `σ`: Volatility coefficient +- `t`: Time (scalar or array) + +### Returns +The standard deviation of log-price at time `t`. +""" +function bs_log_std(σ, t) + return σ * sqrt(t) +end + +# ============================================================================= +# Heston Model: dS = μS dt + √v S dW₁, dv = κ(Θ-v)dt + σ√v dW₂ +# ============================================================================= + +@doc doc""" +Analytical mean of the Heston model asset price. + +```math +E[S(t)] = S_0 e^{\mu t} +``` + +Note: The mean is independent of the stochastic volatility path because +the volatility term has zero expected increment. + +### Arguments +- `μ`: Drift coefficient +- `u0`: Initial asset price ``S_0`` (first component of initial state) +- `t`: Time (scalar or array) + +### Returns +The expected asset price at time `t`. +""" +function heston_mean(μ, u0, t) + return u0 * exp(μ * t) +end + +@doc doc""" +Analytical mean of the Heston model variance process. + +The variance process follows CIR dynamics: ``dv = \kappa(\Theta - v)dt + \sigma\sqrt{v}dW_2`` + +```math +E[v(t)] = \Theta + (v_0 - \Theta) e^{-\kappa t} +``` + +### Arguments +- `κ`: Mean reversion speed +- `Θ`: Long-run variance level +- `v0`: Initial variance ``v_0`` (second component of initial state) +- `t`: Time (scalar or array) + +### Returns +The expected variance at time `t`. +""" +function heston_variance_mean(κ, Θ, v0, t) + return cir_mean(κ, Θ, v0, t) +end + +@doc doc""" +Analytical variance of the Heston model variance process. + +The variance process follows CIR dynamics. + +### Arguments +- `κ`: Mean reversion speed +- `Θ`: Long-run variance level +- `σ`: Vol-of-vol coefficient +- `v0`: Initial variance ``v_0`` +- `t`: Time (scalar or array) + +### Returns +The variance of the variance process at time `t`. +""" +function heston_variance_variance(κ, Θ, σ, v0, t) + return cir_variance(κ, Θ, σ, v0, t) +end diff --git a/test/nopre/jet_tests.jl b/test/nopre/jet_tests.jl index 5936bee..68327ec 100644 --- a/test/nopre/jet_tests.jl +++ b/test/nopre/jet_tests.jl @@ -60,4 +60,28 @@ using Test 0.1, 0.5, 0.2, 0.0, 0.4 ) end + + @testset "Analytical functions type stability" begin + @test_opt target_modules = (DiffEqFinancial,) gbm_mean(0.05, 100.0, 1.0) + @test_opt target_modules = (DiffEqFinancial,) gbm_variance(0.05, 0.2, 100.0, 1.0) + @test_opt target_modules = (DiffEqFinancial,) gbm_std(0.05, 0.2, 100.0, 1.0) + @test_opt target_modules = (DiffEqFinancial,) ou_mean(0.5, 1.0, 0.5, 2.0) + @test_opt target_modules = (DiffEqFinancial,) ou_variance(0.5, 0.3, 2.0) + @test_opt target_modules = (DiffEqFinancial,) ou_std(0.5, 0.3, 2.0) + @test_opt target_modules = (DiffEqFinancial,) ou_stationary_mean(1.0) + @test_opt target_modules = (DiffEqFinancial,) ou_stationary_variance(0.5, 0.3) + @test_opt target_modules = (DiffEqFinancial,) cir_mean(0.5, 0.04, 0.03, 2.0) + @test_opt target_modules = (DiffEqFinancial,) cir_variance(0.5, 0.04, 0.1, 0.03, 2.0) + @test_opt target_modules = (DiffEqFinancial,) cir_std(0.5, 0.04, 0.1, 0.03, 2.0) + @test_opt target_modules = (DiffEqFinancial,) cir_stationary_mean(0.04) + @test_opt target_modules = (DiffEqFinancial,) cir_stationary_variance(0.5, 0.04, 0.1) + @test_opt target_modules = (DiffEqFinancial,) bs_log_mean(0.05, 0.2, 4.6, 1.0) + @test_opt target_modules = (DiffEqFinancial,) bs_log_variance(0.2, 1.0) + @test_opt target_modules = (DiffEqFinancial,) bs_log_std(0.2, 1.0) + @test_opt target_modules = (DiffEqFinancial,) heston_mean(0.05, 100.0, 1.0) + @test_opt target_modules = (DiffEqFinancial,) heston_variance_mean(2.0, 0.04, 0.04, 1.0) + @test_opt target_modules = (DiffEqFinancial,) heston_variance_variance( + 2.0, 0.04, 0.3, 0.04, 1.0 + ) + end end diff --git a/test/runtests.jl b/test/runtests.jl index 32bf2d6..d30e50b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,6 +50,176 @@ end @testset "Interface Compatibility" begin include("interface_tests.jl") end + + @testset "Analytical Solutions" begin + @testset "GBM analytical functions" begin + μ = 0.05 + σ = 0.2 + u0 = 100.0 + t = 1.0 + + # Test at t=0 (should return initial value for mean, 0 for variance) + @test gbm_mean(μ, u0, 0.0) ≈ u0 + @test gbm_variance(μ, σ, u0, 0.0) ≈ 0.0 atol = 1.0e-15 + + # Test analytical formulas + @test gbm_mean(μ, u0, t) ≈ u0 * exp(μ * t) + @test gbm_variance(μ, σ, u0, t) ≈ u0^2 * exp(2μ * t) * (exp(σ^2 * t) - 1) + @test gbm_std(μ, σ, u0, t) ≈ sqrt(gbm_variance(μ, σ, u0, t)) + + # Test broadcasting with array of times + times = [0.0, 0.5, 1.0, 2.0] + means = gbm_mean.(μ, u0, times) + @test length(means) == 4 + @test means[1] ≈ u0 + + # Monte Carlo validation + prob = GeometricBrownianMotionProblem(μ, σ, u0, (0.0, t)) + monte_prob = EnsembleProblem(prob) + sol = solve(monte_prob, EM(); dt = 0.01, trajectories = 10000) + final_values = [sol[i].u[end] for i in eachindex(sol)] + simulated_mean = mean(final_values) + simulated_var = var(final_values) + + analytical_mean = gbm_mean(μ, u0, t) + analytical_var = gbm_variance(μ, σ, u0, t) + + # Allow 5% relative error for Monte Carlo + @test abs(simulated_mean - analytical_mean) / analytical_mean < 0.05 + @test abs(simulated_var - analytical_var) / analytical_var < 0.15 + end + + @testset "OU analytical functions" begin + a = 0.5 + r = 1.0 + σ = 0.3 + u0 = 0.5 + t = 2.0 + + # Test at t=0 + @test ou_mean(a, r, u0, 0.0) ≈ u0 + @test ou_variance(a, σ, 0.0) ≈ 0.0 atol = 1.0e-15 + + # Test analytical formulas + @test ou_mean(a, r, u0, t) ≈ r + (u0 - r) * exp(-a * t) + @test ou_variance(a, σ, t) ≈ (σ^2 / (2a)) * (1 - exp(-2a * t)) + @test ou_std(a, σ, t) ≈ sqrt(ou_variance(a, σ, t)) + + # Test stationary values + @test ou_stationary_mean(r) == r + @test ou_stationary_variance(a, σ) ≈ σ^2 / (2a) + + # Test convergence to stationary distribution + large_t = 100.0 + @test ou_mean(a, r, u0, large_t) ≈ ou_stationary_mean(r) atol = 1.0e-10 + @test ou_variance(a, σ, large_t) ≈ ou_stationary_variance(a, σ) atol = 1.0e-10 + + # Monte Carlo validation + prob = OrnsteinUhlenbeckProblem(a, r, σ, u0, (0.0, t)) + monte_prob = EnsembleProblem(prob) + sol = solve(monte_prob, EM(); dt = 0.01, trajectories = 10000) + final_values = [sol[i].u[end] for i in eachindex(sol)] + simulated_mean = mean(final_values) + simulated_var = var(final_values) + + analytical_mean = ou_mean(a, r, u0, t) + analytical_var = ou_variance(a, σ, t) + + @test abs(simulated_mean - analytical_mean) < 0.05 + @test abs(simulated_var - analytical_var) / analytical_var < 0.15 + end + + @testset "CIR analytical functions" begin + κ = 0.5 + θ = 0.04 + σ = 0.1 + u0 = 0.03 + t = 2.0 + + # Verify Feller condition for test parameters + @test 2κ * θ >= σ^2 # Should satisfy Feller condition + + # Test at t=0 + @test cir_mean(κ, θ, u0, 0.0) ≈ u0 + @test cir_variance(κ, θ, σ, u0, 0.0) ≈ 0.0 atol = 1.0e-15 + + # Test analytical formulas + exp_κt = exp(-κ * t) + exp_2κt = exp(-2κ * t) + expected_mean = θ + (u0 - θ) * exp_κt + expected_var = u0 * (σ^2 / κ) * (exp_κt - exp_2κt) + + (θ * σ^2 / (2κ)) * (1 - exp_κt)^2 + + @test cir_mean(κ, θ, u0, t) ≈ expected_mean + @test cir_variance(κ, θ, σ, u0, t) ≈ expected_var + @test cir_std(κ, θ, σ, u0, t) ≈ sqrt(expected_var) + + # Test stationary values + @test cir_stationary_mean(θ) == θ + @test cir_stationary_variance(κ, θ, σ) ≈ θ * σ^2 / (2κ) + + # Test convergence to stationary distribution + large_t = 100.0 + @test cir_mean(κ, θ, u0, large_t) ≈ cir_stationary_mean(θ) atol = 1.0e-10 + + # Monte Carlo validation + prob = CIRProblem(κ, θ, σ, u0, (0.0, t)) + monte_prob = EnsembleProblem(prob) + sol = solve(monte_prob, EM(); dt = 0.001, trajectories = 10000) + final_values = [sol[i].u[end] for i in eachindex(sol)] + simulated_mean = mean(final_values) + simulated_var = var(final_values) + + analytical_mean = cir_mean(κ, θ, u0, t) + analytical_var = cir_variance(κ, θ, σ, u0, t) + + @test abs(simulated_mean - analytical_mean) / analytical_mean < 0.05 + @test abs(simulated_var - analytical_var) / analytical_var < 0.2 + end + + @testset "Black-Scholes log-price functions" begin + r = 0.05 + σ = 0.2 + u0 = log(100.0) # log of initial price + t = 1.0 + + # Test at t=0 + @test bs_log_mean(r, σ, u0, 0.0) ≈ u0 + @test bs_log_variance(σ, 0.0) ≈ 0.0 atol = 1.0e-15 + + # Test analytical formulas + @test bs_log_mean(r, σ, u0, t) ≈ u0 + (r - σ^2 / 2) * t + @test bs_log_variance(σ, t) ≈ σ^2 * t + @test bs_log_std(σ, t) ≈ σ * sqrt(t) + + # Test broadcasting + times = [0.25, 0.5, 1.0] + variances = bs_log_variance.(σ, times) + @test variances ≈ σ^2 .* times + end + + @testset "Heston analytical functions" begin + μ = 0.05 + κ = 2.0 + Θ = 0.04 + σ = 0.3 + u0_S = 100.0 + u0_v = 0.04 + t = 1.0 + + # Test Heston mean (should match GBM mean since volatility has zero mean increment) + @test heston_mean(μ, u0_S, t) ≈ u0_S * exp(μ * t) + @test heston_mean(μ, u0_S, 0.0) ≈ u0_S + + # Test variance process mean (should match CIR mean) + @test heston_variance_mean(κ, Θ, u0_v, t) ≈ cir_mean(κ, Θ, u0_v, t) + @test heston_variance_mean(κ, Θ, u0_v, 0.0) ≈ u0_v + + # Test variance of variance process (should match CIR variance) + @test heston_variance_variance(κ, Θ, σ, u0_v, t) ≈ + cir_variance(κ, Θ, σ, u0_v, t) + end + end end if GROUP == "all" || GROUP == "nopre"