0. Using 2nd order divided differences is error-prone: round-off errors, catastrophic cancellation, spurious oscillation. Complicated and tedious mathematics is needed in order to verify the Carr/Madan conditions for no-arbitrage.

1. When calculating sensitivities we choose for differentiation-then-approximation approach (1) to the BS PDE. In particular, each sensitivity satisfies a BS-type equation with a possible inhomogeneous forcing term + BC + IC (2).

2. Instead of (trying to) compute the solution of the PDE analytically we use the maximum principle (positive input implies positive output/solution (not many quants use/know this). It is crucial and it allows us to give sufficient conditions for no-arbitrage without much effort at all. The result is the continuous analogue of the Carr/Madan conditions.

3. The next step is to approximate the solution of the PDE by monotone finite difference schemes that satisfy the discrete maximum principle for any mesh size and set of input parameters (not all schemes are monotonic and L-stable which is the reason for all that hullabaloo around Crank-Nicolson). CN gives oscillations around the strike, making it unsuitable as it can cause discontinuities in LV function. Fully implicit is L-stable but not perfect. It is only 1

^{st}-order.

4. For convection-dominant problems or tricky diffusion terms, we can use exponential fitting to preserve monotonicity and avoid spatial instability.

5. Example; gamma and strike gamma solve Fokker-Planck PDE and you can see that the approximate integral adds up to 1. The initial condition is a Dirac delta function. The maximum principle can be applied to prove positivity of this PDE.

6. For sensitivity ‘theta’ it is non-decreasing in maturity as can be seen by looking directly at the Dupire PDE. No need for messing with integrals which I find to be a bit grungy TBH. It soon becomes a Pandora’s box (what’s the Laplace-Carson transform?)

7. I am fairly confident that the same techniques can be applied to compute duration and convexity sensitivities etc. for IR models. I suspect it should also be possible to compute swaption and CVA greeks assuming you have a PDE for the underlying.

8. I don’t know yet how this approach fits into the LV calibration algorithm. At this stage, we can at least say that the volatility in Dupire’s formula will always be positive no matter what. Not sure if piecewise constant vol or cubic splines are necessarily superior.

(1) This is in essence the Continuous Sensitivity Equation (CSE) approach and is used in shape design and optimization.

(2) Many current approaches use the approximation-then- differentiation approach, for example discretization of a SDE by Euler and then differentiating the discrete equation wrt vol using AD or the complex step method.