Finite Difference Method in Greeks (Options)

I need a way to approximate the analytical formula of Greeks of a generic call option using the Finite Difference Method.

For example, the FD method for Delta/Gamma is the following one:

Now, I am in trouble with respect to the denominator "DeltaS"; how can I find the optimal value that minimize the distance between the analytical formula of Delta/Gamma obtained with Black-Scholes?

• If your C() is standard black scholes, the smaller the closer by definition of the derivative Jul 24, 2021 at 14:11
• Yes, you are right for the first derivative, but we get a worse second derivative with smaller value of deltaS @AKdemy Jul 24, 2021 at 15:19
• I suspect this is because either (1) you are using values of $C(\dots)$ from a tree or grid, giving rise to discretization errors, or (2) you are setting $\Delta S$ so small that you are running into numerical problems with your approximation to the gaussian CDF. Jul 24, 2021 at 17:20

1 Answer

Agree with @Brian B. With BS, you cannot have the issue in (1). Tree, grid, Monte Carlo could all result in errors though.

(2) is a likely reason. I just tried in Julia for ATM, 0 div and rates plus 0.2 vol and 1 year tenor. Shifts smaller than ~ 0.00008 result in an error for Gamma. Delta seems to be less sensitive for this, and it is fine for at least 1e-7 and deviates for 1e-8. So anywhere in between.

I don't think there is a serious issue that requires much thought to be put into the exact point when it deviates. At least not if you ask about Black Scholes. In case of Monte Carlo pricing for example, that will be quite important as you don't want to end up shifting in an area inside your standard error which will just be noise.

Personally, I prefer shifting up and down (central difference) for most greeks. Based on my experience, this seems to also be the consensus (or at least most frequently used implementation). Below is an intuitive explanation why I think it is better (compared to your forward difference -> only shifting up).

Consider this toy example where BSM is a custom function for generic Black Scholes where first $$[1]$$ index provides the call option value and $$[2]$$ the delta:

K=10 # strike
t = 1 # 1 year
d = 0 # zero dividends
rf = 0  # zero rates
σ = 0.2 # 20% IVOL

function deltaBumpReprice(S,bump)
up = BSM(S+bump/2,K, t, rf, d, σ)[1]
down = BSM(S-bump/2,K, t, rf, d, σ)[1]
delta = BSM(S,K, t, rf, d, σ)[2]
approx = (up-down)/bump
difference = delta-approx
return approx, delta, difference
end


vs your single shift up

function deltaBumpRepriceqse(S,bump)
up = BSM(S+bump,K, t, rf, d, σ)[1]
down = BSM(S,K, t, rf, d, σ)[1]
delta = BSM(S,K, t, rf, d, σ)[2]
approx = (up-down)/bump
difference = delta-approx
return approx, delta, difference
end


Now assume we are ATMS (S=K=10) and shift in integers (1,2,3,..., 20) which is obviously extreme. The dataframe shows approximate delta with up/down, delta, the difference between the two, a single upshift approximation and the difference to delta.

You can see that even with crazy shifts, delta with shifting up and down is still sort of "reasonable". How come?

This bump corresponds to a spot of 4.7 and 14.7 respectively, as opposed to a spot of 9.7 where analytical delta is computed. Ignore the shift number, that is a simplification as I used the DataFrame index directly. Yet, the approximation is not too bad. The chart also shows what happens for such a large shift in your implementation.

Obviously this is unrealistic, but bump up and down will always be better than shifting in one direction. Below is an example for an OTM call.

Lastly, redo the same exercise as above for very small shifts from 0.0000001 up to 0.001, in 0.00005 shifts.

• Thank you very much!! Jul 26, 2021 at 7:25
• On a side note, actual shifts are usually a lot bigger than the small shifts shown above (e.g. 1% of spot: $shift = spot * 1\%/ 2$). P. Glasserman. Monte Carlo Methods in Financial Engineering, 2010 as well as [M. Henrard. Sensitivity computation. OpenGamma (July 2014)](developers. opengamma.com/quantitative-research/sensitivity-computation.pdf) are frequently referenced sources on determining bump sizes to avoid numerical instabilities associated with shift sizes that are too small or too large. Jul 26, 2021 at 11:53