Merge pull request #4 from aleCombi/TestSuiteRedesign

Test suite redesign

Author: aleCombi

Project.toml

@@ -5,18 +5,13 @@ version = "1.0.0-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
-DiffEqGPU = "071ae1c0-96b5-11e9-1965-c90190d839ea"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
-JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
@@ -26,23 +21,16 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Accessors = "0.1.42"
-BenchmarkTools = "1.6.0"
-Bessels = "0.2.8"
 DataInterpolations = "7.2.0"
 Dates = "1.11.0"
-DiffEqGPU = "3.4.1"
 DiffEqNoiseProcess = "5.24.1"
 DifferentialEquations = "7.16.0"
 Distributions = "0.25"
-FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.38"
 Integrals = "4.5.0"
-JuliaFormatter = "1.0.62"
 NonlinearSolve = "4.5.0"
 Optimization = "4.1.2"
 Polynomials = "4.0.19"
@@ -53,7 +41,6 @@ SpecialFunctions = "2.5.0"
 StaticArrays = "1.9.13"
 Statistics = "1.11.1"
 Test = "1.11.0"
-Zygote = "0.6.76"
 julia = "1.11"
 
 [extras]

examples/calib.jl

@@ -0,0 +1,22 @@
+using Revise, Dates, Hedgehog2
+
+reference_date = Date(2020, 1, 1)
+r, S0, sigma = 0.05, 100.0, 0.25
+market_inputs = BlackScholesInputs(reference_date, r, S0, sigma)
+
+strikes = collect(60.0:5.0:140.0)
+expiry = reference_date + Day(365)
+payoffs = [VanillaOption(K, expiry, European(), Call(), Spot()) for K in strikes]
+
+quotes = [
+    solve(PricingProblem(p, market_inputs), BlackScholesAnalytic()).price for
+    p in payoffs
+]
+
+accessors = [VolLens(1,1)]
+initial_guess = [0.15]
+
+basket = BasketPricingProblem(payoffs, market_inputs)
+calib = CalibrationProblem(basket, BlackScholesAnalytic(), accessors, quotes, 0.7*ones(length(payoffs)))
+result = solve(calib, OptimizerAlgo())
+@show result.objective

examples/calibration_heston.jl

@@ -48,11 +48,11 @@ initial_guess = [0.02, 3.0, 0.03, 0.4, -0.3]  # [v0, κ, θ, σ, ρ]
 
 # Accessors for each parameter in HestonInputs
 accessors = [
-    @optic(_.market.V0),
-    @optic(_.market.κ),
-    @optic(_.market.θ),
-    @optic(_.market.σ),
-    @optic(_.market.ρ)
+    @optic(_.market_inputs.V0),
+    @optic(_.market_inputs.κ),
+    @optic(_.market_inputs.θ),
+    @optic(_.market_inputs.σ),
+    @optic(_.market_inputs.ρ)
 ]
 
 

examples/montecarlo_blackscholes.jl

@@ -1,3 +1,4 @@
+using Revise
 using Hedgehog2
 using Dates
 using Printf
@@ -27,14 +28,15 @@ bs_solution = solve(prob, bs_method)
 # Using 100,000 paths for more accurate results
 trajectories = 5_000
 mc_exact_method =
-    MonteCarlo(LognormalDynamics(), BlackScholesExact(trajectories, antithetic = true))
+    MonteCarlo(LognormalDynamics(), BlackScholesExact(), SimulationConfig(trajectories))
 mc_exact_solution = solve(prob, mc_exact_method)
 
 # Method 3: Monte Carlo with Euler-Maruyama discretization
 # Using 100,000 paths and 100 time steps
+trajectories = 10_000
 steps = 100
 mc_euler_method =
-    MonteCarlo(LognormalDynamics(), EulerMaruyama(trajectories, steps, antithetic = true))
+    MonteCarlo(LognormalDynamics(), EulerMaruyama(), SimulationConfig(trajectories, steps=steps, variance_reduction=Hedgehog2.Antithetic()))
 mc_euler_solution = solve(prob, mc_euler_method)
 
 # Print the results

src/Hedgehog2.jl

@@ -19,11 +19,11 @@ end
 
 if false
     include("../test/runtests.jl")
+    include("../test/calibration.jl")
 end
 
 # utilities
 include("date_functions.jl")
-include("solutions/pricing_solutions.jl")
 
 # payoffs
 include("payoffs/payoffs.jl")
@@ -35,6 +35,7 @@ include("market_inputs/market_inputs.jl")
 
 # pricing methods
 include("pricing_methods/pricing_methods.jl")
+include("solutions/pricing_solutions.jl")
 include("pricing_methods/black_scholes.jl")
 include("pricing_methods/cox_ross_rubinstein.jl")
 include("pricing_methods/montecarlo.jl")

src/date_functions.jl

@@ -1,21 +1,27 @@
 const SECONDS_IN_YEAR_365 = 365 * 86400
 const MILLISECONDS_IN_YEAR_365 = SECONDS_IN_YEAR_365 * 1000
+const MILLISECONDS_IN_DAY = 86400000
 
 # --- Time Conversions ---
 
 """
     to_ticks(x::Date)
 
-Convert a `Date` to milliseconds since the Unix epoch.
+Convert a `Date` to milliseconds since the Julia `Dates` module epoch (0000-01-01).
+
+Note: This is calculated by converting the `Date` to days since epoch 
+(`Dates.date2epochdays`) and multiplying by `MILLISECONDS_IN_DAY`.
 """
 function to_ticks(x::Date)
-    return Dates.datetime2epochms(DateTime(x))
+    return Dates.date2epochdays(x) * MILLISECONDS_IN_DAY
 end
 
 """
     to_ticks(x::DateTime)
 
-Convert a `DateTime` to milliseconds since the Unix epoch.
+Convert a `DateTime` to milliseconds since the Julia `Dates` module epoch (0000-01-01T00:00:00).
+
+Uses `Dates.datetime2epochms`.
 """
 function to_ticks(x::DateTime)
     return Dates.datetime2epochms(x)
@@ -24,8 +30,11 @@ end
 """
     to_ticks(x::Real)
 
-Assume `x` is already a timestamp in milliseconds (e.g., `Float64` or `Int`).
-Used to normalize mixed inputs.
+Assume `x` is already a timestamp in milliseconds since the Julia `Dates` module epoch
+(0000-01-01T00:00:00) and return it unchanged.
+
+Used to normalize mixed inputs (e.g., `Date`, `DateTime`, `Real`) to a common
+tick representation for calculations like `yearfrac` or `add_yearfrac`.
 """
 function to_ticks(x::Real)
     return x
@@ -38,7 +47,9 @@ end
 
 Compute the ACT/365 year fraction between two time points.
 
-Supports `Date`, `DateTime`, or ticks (`Int` or `Float64`).
+Inputs `start` and `stop` can be `Date`, `DateTime`, or ticks (`Int` or `Float64`).
+If ticks are provided, they are assumed to be milliseconds since the Julia `Dates`
+module epoch (0000-01-01T00:00:00), consistent with the output of `to_ticks`.
 """
 function yearfrac(start, stop)
     ms_start = to_ticks(start)
@@ -47,12 +58,17 @@ function yearfrac(start, stop)
 end
 
 """
-    yearfrac(p::Period)
+    yearfrac(p::AbstractTime)
 
-Compute the ACT/365 year fraction from a `Period` object.
+Compute the ACT/365 year fraction from a `Period` object (e.g., `Year(1)`, `Day(180)`).
+It allows also for `CompoundPeriod`, and even Dates, interpreting them as time periods from year 0 of the Gregorian calendar.
+Calculates the duration represented by the period in milliseconds and divides by
+`MILLISECONDS_IN_YEAR_365`.
 """
-function yearfrac(p::Period)
-    ref = DateTime(1970, 1, 1)
+function yearfrac(p::Dates.AbstractTime)
+    # The choice of reference date here is arbitrary and does not affect the result,
+    # as only the difference (duration) matters.
+    ref = DateTime(1970, 1, 1) 
     return yearfrac(ref, ref + p)
 end
 
@@ -61,8 +77,10 @@ end
 """
     add_yearfrac(t::Real, yf::Real) -> Real
 
-Add a fractional number of years (ACT/365) to a timestamp in milliseconds since epoch.
-Returns the updated timestamp in milliseconds as `Float64`.
+Add a fractional number of years (`yf`, computed as ACT/365) to a timestamp `t`.
+
+The timestamp `t` is assumed to be in milliseconds since the Julia `Dates` module
+epoch (0000-01-01T00:00:00). Returns the updated timestamp in milliseconds as `Float64`.
 
 This version is AD-compatible.
 """
@@ -73,9 +91,14 @@ end
 """
     add_yearfrac(t::TimeType, yf::Real) -> DateTime
 
-Add a fractional number of years (ACT/365) to a `Date` or `DateTime`.
-Returns a `DateTime` object.
+Add a fractional number of years (`yf`, computed as ACT/365) to a `Date` or `DateTime` `t`.
+
+Converts `t` to milliseconds since the Julia `Dates` module epoch, adds the
+duration corresponding to `yf`, and converts the result back to a `DateTime` object.
 """
 function add_yearfrac(t::TimeType, yf::Real)
+    # `to_ticks(t)` uses the Julia Dates epoch (Year 0000)
+    # `add_yearfrac(Real,Real)` adds the correct millisecond duration
+    # `epochms2datetime` converts back from the Julia Dates epoch (Year 0000)
     return Dates.epochms2datetime(add_yearfrac(to_ticks(t), yf))
-end
+end

src/greeks/greeks_problem.jl

@@ -63,7 +63,6 @@ end
 
 function compute_fd_derivative(::FDCentral, prob, lens, ε, pricing_method)
     x₀ = lens(prob)
-
     prob_up = set(prob, lens, x₀ * (1 + ε))
     prob_down = set(prob, lens, x₀ * (1 - ε))
     v_up = solve(prob_up, pricing_method).price
@@ -173,16 +172,14 @@ function solve(
         # Delta = ∂V/∂S = ∂V/∂F * ∂F/∂S = (cp * N(cp·d1)) * (1/D)
         cp * Φ(cp * d1)
 
-    elseif lens === @optic _.market_inputs.sigma
+    elseif lens === VolLens(1,1) #TODO: add logic
         # Vega = ∂V/∂σ = D · F · φ(d1) · √T
         D * F * ϕ(d1) * √T
 
     elseif lens === @optic _.payoff.expiry
-        @assert is_flat(prob.market_inputs.rate)
-
         # Assume flat rate: z(T) = r ⇒ D(T) = exp(-rT), F(T) = S / D(T)
         r = zero_rate_yf(prob.market_inputs.rate, T)
-        (r * prob.payoff.strike * D * Φ(d2) + F * D * prob.market_inputs.sigma * ϕ(d1) / (2√T)) / (MILLISECONDS_IN_YEAR_365) #against ticks, to match AD and FD. Observe that the sign is counterintuitive as it is a derivative against expiry in tticks, not against time-to-maturity in yearfrac
+        (r * prob.payoff.strike * D * Φ(d2) + F * D * σ * ϕ(d1) / (2√T)) / (MILLISECONDS_IN_YEAR_365) #against ticks, to match AD and FD. Observe that the sign is counterintuitive as it is a derivative against expiry in tticks, not against time-to-maturity in yearfrac
 
     else
         error("Unsupported lens for analytic Greek")
@@ -198,9 +195,9 @@ function solve(gprob::SecondOrderGreekProblem, ::AnalyticGreek, ::BlackScholesAn
     inputs = prob.market_inputs
 
     S = inputs.spot
-    σ = inputs.sigma
     T = yearfrac(inputs.referenceDate, prob.payoff.expiry)
     K = prob.payoff.strike
+    σ = get_vol_yf(inputs.sigma, T, K)
 
     D = df(inputs.rate, prob.payoff.expiry)
     F = S / D
@@ -214,7 +211,7 @@ function solve(gprob::SecondOrderGreekProblem, ::AnalyticGreek, ::BlackScholesAn
         # Gamma = ∂²V/∂S² = φ(d1) / (Sσ√T)
         ϕ(d1) / (S * σ * √T)
 
-    elseif (lens1 === @optic _.market_inputs.sigma) && (lens2 === @optic _.market_inputs.sigma)
+    elseif (lens1 === VolLens(1,1)) && (lens2 === VolLens(1,1)) #TODO: introduce logic for sigma
         # Volga = Vega * d1 * d2 / σ
         vega = D * F * ϕ(d1) * √T
         vega * d1 * d2 / σ
@@ -223,7 +220,7 @@ function solve(gprob::SecondOrderGreekProblem, ::AnalyticGreek, ::BlackScholesAn
         error("Unsupported second-order analytic Greek")
     end
 
-    return GreekResult(deriv)
+    return GreekResult(greek)
 end
 
 struct BatchGreekProblem{P,L}

src/market_inputs/rate_curve.jl

@@ -64,7 +64,7 @@ zero_rate(curve::RateCurve, ticks::T) where T <: Number =
 zero_rate(curve::FlatRateCurve, ticks::T) where T <: Number =
     curve.rate
 
-zero_rate(curve::R, t::Date) where R <: AbstractRateCurve = zero_rate(curve, to_ticks(t))
+zero_rate(curve::R, t::D) where {R <: AbstractRateCurve, D<:TimeType} = zero_rate(curve, to_ticks(t))
 
 zero_rate_yf(curve::RateCurve, yf::R) where R <: Number = curve.interpolator(yf)
 zero_rate_yf(curve::FlatRateCurve, yf::R) where R <: Number = curve.rate

src/market_inputs/vol_surface.jl

@@ -211,7 +211,7 @@ function RectVolSurface(
     end
 
     vols = Matrix{Float64}(undef, nrows, ncols)
-    accessor = @optic _.market_inputs.sigma
+    accessor = VolLens(1,1)
     for i = 1:nrows, j = 1:ncols
         expiry = reference_date + tenors[i]
         strike = strikes[j]

src/pricing_methods/black_scholes.jl

@@ -37,7 +37,7 @@ Computes the price of a European vanilla option under the Black-Scholes model.
 """
 function solve(
     prob::PricingProblem{VanillaOption{TS,TE,European,B,C}, I},
-    ::BlackScholesAnalytic,
+    method::BlackScholesAnalytic,
 ) where {TS,TE,B,C, I <: BlackScholesInputs}
 
     payoff = prob.payoff
@@ -60,5 +60,5 @@ function solve(
         D * cp * (F * cdf(N, cp * d1) - K * cdf(N, cp * d2))
     end
 
-    return AnalyticSolution(price)
+    return AnalyticSolution(prob, method, price)
 end