Testing heston exact against heston

Author: aleCombi

examples/Project.toml

@@ -10,4 +10,5 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

examples/heston_noise.jl

@@ -1,11 +1,11 @@
 using Revise, Hedgehog2, Distributions, DifferentialEquations, Random, Plots
 
 # Define Heston model parameters
-S0 = 100.0    # Initial stock price
+S0 = 1.0    # Initial stock price
 V0 = 0.04     # Initial variance
 κ = 2.0       # Mean reversion speed
 θ = 0.04      # Long-run variance
-σ = 0.3       # Volatility of variance
+σ = 0.02     # Volatility of variance
 ρ = -0.7      # Correlation
 r = 0.05      # Risk-free rate
 T = 1.0       # Time to maturity
@@ -19,9 +19,19 @@ heston_noise = Hedgehog2.HestonNoise(0.0, heston_dist)
 # Define `NoiseProblem`
 problem = NoiseProblem(heston_noise, (0.0, T))
 
+trajectories = 10000
 # Solve with multiple trajectories
-solution = solve(EnsembleProblem(problem), dt=T, trajectories=10)
+solution_h = solve(EnsembleProblem(problem), dt=T, trajectories=trajectories)
 
-# Directly plot the solution
-plot(solution, title="Exact Heston Model Simulations",
-     xlabel="Time", ylabel="Log Stock Price")
+u0 = [S0, V0]  # Initial (Stock Price, Variance)
+tspan = (0.0, 1.0)  # Simulation for 1 year
+process = HestonProcess(r, κ, θ, σ, ρ)  # Heston parameters
+prob = SDEProblem(get_sde_function(process), u0, tspan, process)
+ensemble_problem = EnsembleProblem(prob)
+solution_h2 = solve(ensemble_problem, dt=T, trajectories=trajectories)
+
+final_prices_2 = [sol.u[end][1] for sol in solution_h2]  # Extract S_T for each trajectory
+final_prices = [sol.u[end] for sol in solution_h]  # Transform log-prices to prices
+
+println(mean(final_prices_2.^3))
+println(mean(final_prices.^3))

examples/montecarlo_bs.jl

@@ -15,7 +15,7 @@ strike = 1.2
 payoff = VanillaOption(strike, expiry, Hedgehog2.European(), Hedgehog2.Call(), Hedgehog2.Spot())
 
 # Define montecarlo methods
-trajectories = 100000000
+trajectories = 10000
 method = Hedgehog2.BSMontecarlo(trajectories)
 
 # Define pricer

src/distributions/heston.jl

@@ -30,21 +30,20 @@ end
 
 Samples `log(S_T)` using the exact Broadie-Kaya method.
 """
-function rand(rng::AbstractRNG, d::HestonDistribution)
+function rand(rng, d::HestonDistribution)
     κ, θ, σ, ρ, V0, T, S0, r = d.κ, d.θ, d.σ, d.ρ, d.V0, d.T, d.S0, d.r
-
     # 1. Sample V_T using the noncentral chi-squared distribution
-    ν = (4κθ) / σ^2
+    ν = (4κ * θ) / σ^2
     ψ = (4κ * exp(-κ * T) * V0) / (σ^2 * (1 - exp(-κ * T)))
-    V_T = (σ^2 * (1 - exp(-κ * T)) / (4κ)) * rand(rng, NoncentralChisq(ν, ψ))
+    V_T = (σ^2 * (1 - exp(-κ * T)) / (4κ)) * Distributions.rand(rng, NoncentralChisq(ν, ψ))
 
     # 2. Sample log(S_T) given V_T
     μ_XT = log(S0) + (r - 0.5 * V0) * T + (1 - exp(-κ * T)) * (θ - V0) / (2κ)
     var_XT = (ρ^2 * V0 * (1 - exp(-κ * T))^2) / (2κ) + (1 - ρ^2) * (V0 * T + θ * (T - (1 - exp(-κ * T)) / κ))
 
     log_S_T = μ_XT + sqrt(var_XT) * randn(rng)  # Sample from Normal(μ_XT, sqrt(var_XT))
 
-    return log_S_T
+    return exp(log_S_T)
 end
 
 """

src/pricing_methods/montecarlo.jl

@@ -22,5 +22,5 @@ function compute_price(payoff::VanillaOption{European, C, Spot}, market_inputs::
     final_payoffs = (payoff.(last.(solution.u)) + payoff.(last.(antithetic_solution.u))) * 0.5
     mean_payoff = mean(final_payoffs)
     var_payoff = var(final_payoffs)
-    return mean_payoff, var_payoff
+    return mean_payoff
 end

src/stochastic_processes/heston_noise.jl

@@ -17,7 +17,5 @@ function HestonNoise(t0, heston_dist::HestonDistribution, Z0=nothing)
         return rand(rng, heston_dist)  # Calls exact Heston sampler
     end
 
-    return NoiseProcess(t0, log_S0, Z0, Heston_dist, nothing)
+    return NoiseProcess{false}(t0, log_S0, Z0, Heston_dist, (dW, W, W0, Wh, q, h, u, p, t, rng) -> 1)
 end
-
-