Finite Difference Method in Greeks (Options)

Question

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?

Answer 1

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.