support for Symbolics

src/FinancialToolbox.jl

@@ -13,6 +13,7 @@ function __init__()
     @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" include("financial_reverse_diff.jl")
     @require HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" include("financial_hyper.jl")
     @require TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" include("financial_taylor.jl")
+    @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("financial_symbolics.jl")
 end
 export normcdf,
     normpdf,

src/financial.jl

@@ -357,7 +357,20 @@ function blslambda(S0, K, r, T, σ, d, FlagIsCall::Bool = true)
 end
 
 #Check input for Black Scholes Formula
-function blcheck(S0::num1, K::num2, T::num4, σ::num5) where {num1, num2, num4, num5}
+function blcheck_impl(::num0, S0::num1, K::num2, T::num4, σ::num5) where {num0, num1, num2, num4, num5}
+    lesseq(x, y) = x <= y
+    if (lesseq(S0, zero(num1)))
+        throw(DomainError(S0, "Spot Price Cannot Be Negative"))
+    elseif (lesseq(K, zero(num2)))
+        throw(DomainError(K, "Strike Price Cannot Be Negative"))
+    elseif (lesseq(T, zero(num4)))
+        throw(DomainError(T, "Time to Maturity Cannot Be Negative"))
+    elseif (lesseq(σ, zero(num5)))
+        throw(DomainError(σ, "Volatility Cannot Be Negative"))
+    end
+    return
+end
+function blcheck_impl(::num0, S0::num1, K::num2, T::num4, σ::num5) where {num0 <: Complex, num1, num2, num4, num5}
     lesseq(x::Complex, y::Complex) = real(x) <= real(y)
     lesseq(x, y) = x <= y
     if (lesseq(S0, zero(num1)))
@@ -371,6 +384,12 @@ function blcheck(S0::num1, K::num2, T::num4, σ::num5) where {num1, num2, num4,
     end
     return
 end
+
+function blcheck(S0::num1, K::num2, T::num4, σ::num5) where {num1, num2, num4, num5}
+    zero_typed = S0 * K * T * σ
+    blcheck_impl(zero_typed, S0, K, T, σ)
+    return
+end
 ##	ADDITIONAL Functions
 
 """

src/financial_symbolics.jl

@@ -0,0 +1,72 @@
+using .Symbolics
+
+function blcheck_impl(::num0, S0::num1, K::num2, T::num4, σ::num5) where {num0 <: Num, num1, num2, num4, num5}
+    return
+end
+
+@register_symbolic blimpv(S0, K, T, Price, FlagIsCall::Bool, xtol::Float64, n_iter_max::Int64)
+
+function blvega_impl_registered(S0, K, T, σ)
+    return blvega_impl(S0, K, T, σ)
+end
+
+function Symbolics.derivative(::typeof(blvega_impl_registered), ::NTuple{4, Any}, ::Any)
+    return 0
+end
+function finalize_derivative_fwd(S0, K, T, σ, der_fwd)
+    vega = blvega_impl_registered(S0, K, T, σ)
+    return der_fwd / vega
+end
+# function finalize_derivative_fwd(vega, price_diff, target_var, original_var)
+#     der_fwd = Symbolics.derivative(price_diff, target_var)
+#     der_fwd = substitute(der_fwd, Dict(target_var => original_var))
+#     return der_fwd / vega
+# end
+function bldelta(S0, K, T, σ, FlagIsCall)
+    sigma_sqrtT = σ * sqrt(T)
+    S0_K = S0 / K
+    d1 = log(S0_K) / sigma_sqrtT + sigma_sqrtT / FinancialToolbox.two_internal
+    iscall = ifelse(FlagIsCall, FinancialToolbox.one_internal, FinancialToolbox.minus_one_internal)
+    Δ = normcdf(iscall * d1) * iscall
+    return Δ
+end
+#TODO: implement derivatives
+function Symbolics.derivative(::typeof(blimpv), args::NTuple{7, Any}, ::Val{1})
+    S0, K, T, price_d, FlagIsCall, xtol, n_iter_max = args
+    σ = blimpv(S0, K, T, price_d, FlagIsCall, xtol, n_iter_max)
+    @variables new_S0
+    price_diff = -blprice_impl(new_S0, K, T, σ, FlagIsCall)
+    der_fwd = Symbolics.derivative(price_diff, new_S0)
+    der_fwd = substitute(der_fwd, Dict(new_S0 => S0))
+    #der_fwd = bldelta(S0, K, r, T, σ, d, FlagIsCall)
+    return finalize_derivative_fwd(S0, K, T, σ, der_fwd)
+end
+function Symbolics.derivative(::typeof(blimpv), args::NTuple{7, Any}, ::Val{2})
+    S0, K, T, price_d, FlagIsCall, xtol, n_iter_max = args
+    σ = blimpv(S0, K, T, price_d, FlagIsCall, xtol, n_iter_max)
+    @variables new_K
+    price_diff = -blprice_impl(S0, new_K, T, σ, FlagIsCall)
+    der_fwd = Symbolics.derivative(price_diff, new_K)
+    der_fwd = substitute(der_fwd, Dict(new_K => K))
+    #der_fwd = bldelta(S0, K, r, T, σ, d, FlagIsCall)
+    return finalize_derivative_fwd(S0, K, T, σ, der_fwd)
+end
+function Symbolics.derivative(::typeof(blimpv), args::NTuple{7, Any}, ::Val{3})
+    S0, K, T, price_d, FlagIsCall, xtol, n_iter_max = args
+    σ = blimpv(S0, K, T, price_d, FlagIsCall, xtol, n_iter_max)
+    @variables new_T
+    price_diff = -blprice_impl(S0, K, new_T, σ, FlagIsCall)
+    der_fwd = Symbolics.derivative(price_diff, new_T)
+    der_fwd = substitute(der_fwd, Dict(new_T => T))
+    #der_fwd = bldelta(S0, K, r, T, σ, d, FlagIsCall)
+    return finalize_derivative_fwd(S0, K, T, σ, der_fwd)
+end
+function Symbolics.derivative(::typeof(blimpv), args::NTuple{7, Any}, ::Val{4})
+    S0, K, T, price_d, FlagIsCall, xtol, n_iter_max = args
+    σ = blimpv(S0, K, T, price_d, FlagIsCall, xtol, n_iter_max)
+    vega = blvega_impl_registered(S0, K, T, σ)
+    return inv(vega)
+end
+function Symbolics.derivative(::typeof(blimpv), ::NTuple{7, Any}, ::Any)
+    return 0
+end

test/Project.toml

@@ -6,6 +6,7 @@ DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

test/runtests.jl

@@ -6,7 +6,7 @@ function print_colored(in::String, color1)
     end
 end
 
-test_list = ["testRealNumbers.jl", "test_implied_volatility.jl", "test_implied_volatility_black.jl", "testComplexNumbers.jl", "testForwardDiff.jl", "testReverseDiff.jl", "testDual_.jl", "testHyperDualNumbers.jl", "testDates.jl", "testTaylor.jl", "testZygote.jl", "testDiffractor.jl"]
+test_list = ["testRealNumbers.jl", "test_implied_volatility.jl", "test_implied_volatility_black.jl", "testComplexNumbers.jl", "testForwardDiff.jl", "testReverseDiff.jl", "testDual_.jl", "testHyperDualNumbers.jl", "testDates.jl", "testTaylor.jl", "testZygote.jl", "testDiffractor.jl", "testSymbolics.jl"]
 
 println("Running tests:\n")
 for (current_test, i) in zip(test_list, 1:length(test_list))

test/testSymbolics.jl

@@ -0,0 +1,56 @@
+using Test
+using FinancialToolbox, Symbolics
+toll = 1e-7
+@variables S0_s, r_s, d_s, T_s, sigma_s, K_s, price_s
+S0 = 100.0;
+K = 100.0;
+r = 0.02;
+T = 1.2;
+sigma = 0.2;
+d = 0.01;
+price = blsprice(S0, K, r, T, sigma, d)
+dict_vals = Dict(S0_s => S0, r_s => r, d_s => d, K_s => K, T_s => T, sigma_s => sigma, price_s => price)
+price_s_c = blsprice(S0_s, K_s, r_s, T_s, sigma_s, d_s)
+price_v = substitute(price_s_c, dict_vals)
+@test abs(price - price_v) < toll
+
+vol_s = blsimpv(S0_s, K_s, r_s, T_s, price_s, d_s)
+vol_v = substitute(vol_s, dict_vals)
+vol = blsimpv(S0, K, r, T, price, d)
+@test abs(vol - vol_v) < toll
+
+#Test multiple order derivatives of blsimpv
+using HyperDualNumbers
+#S0
+der_vol_S0 = Symbolics.derivative(vol_s, S0_s, simplify = true);
+der_vol_S02 = Symbolics.derivative(der_vol_S0, S0_s, simplify = true);
+vol_h_S0 = blsimpv(hyper(S0, 1.0, 1.0, 0.0), K, r, T, price, d)
+delta = substitute(der_vol_S0, dict_vals)
+gamma = substitute(der_vol_S02, dict_vals)
+@test abs(delta - vol_h_S0.epsilon1) < toll
+@test abs(gamma - vol_h_S0.epsilon12) < toll
+
+#K
+der_vol_K = Symbolics.derivative(vol_s, K_s, simplify = true);
+der_vol_K2 = Symbolics.derivative(der_vol_K, K_s, simplify = true);
+vol_h_K = blsimpv(S0, hyper(K, 1.0, 1.0, 0.0), r, T, price, d)
+delta_k = substitute(der_vol_K, dict_vals)
+gamma_k = substitute(der_vol_K2, dict_vals)
+@test abs(delta_k - vol_h_K.epsilon1) < toll
+@test abs(gamma_k - vol_h_K.epsilon12) < toll
+#r
+der_vol_r = Symbolics.derivative(vol_s, r_s, simplify = true);
+der_vol_r2 = Symbolics.derivative(der_vol_r, r_s, simplify = true);
+vol_h_r = blsimpv(S0, K, hyper(r, 1.0, 1.0, 0.0), T, price, d)
+delta_r = substitute(der_vol_r, dict_vals)
+gamma_r = substitute(der_vol_r2, dict_vals)
+@test abs(delta_r - vol_h_r.epsilon1) < toll
+@test abs(gamma_r - vol_h_r.epsilon12) < toll
+#T
+der_vol_T = Symbolics.derivative(vol_s, T_s, simplify = true);
+der_vol_T2 = Symbolics.derivative(der_vol_T, T_s, simplify = true);
+vol_h_T = blsimpv(S0, K, r, hyper(T, 1.0, 1.0, 0.0), price, d)
+delta_T = substitute(der_vol_T, dict_vals)
+gamma_T = substitute(der_vol_T2, dict_vals)
+@test abs(delta_T - vol_h_T.epsilon1) < toll
+@test abs(gamma_T - vol_h_T.epsilon12) < toll