Finite Difference method in Matlab for SABR volatility model fails to provide correct option values

Question

Currently, I'm trying to implement a Finite Difference (FD) method in Matlab for my thesis (Quantitative Finance). I implemented the FD method for Black-Scholes already and got correct results. However, I want to extend it to work for the SABR volatility model. Although some information on this model can be found on the internet, this mainly regards Hagan's approximate formula for EU options. I am particularly interested in the numerical solution (FD).

I have taken the following steps already:

Substitute the derivatives in the SABR PDE with their finite difference approximations. I used the so-called implicit method (https://en.wikipedia.org/wiki/Finite_difference_method#Implicit_method) for this. Possibility for variable transformation from $F$ to $x$ and $\alpha$ to $y$ is incorporated in PDE discretization, but not applied yet (hence, $\dfrac{\partial x}{\partial F}=\dfrac{\partial y}{\partial \alpha}=1$ and $\dfrac{\partial^2 x}{\partial F^2}=\dfrac{\partial^2 y}{\partial \alpha^2}=0$ in formulas below). Substituting forward difference in time and central differences in space dimensions and rewriting, gives me the following equation: $$ V_{i,j,k} = V_{i-1,j-1,k+1}[-c] + V_{i,j-1,k+1}[-d+e+g] + V_{i+1,j-1,k+1}[c] + V_{i-1,j,k+1}[-a+b+f] + V_{i,j,k+1}[2a+2d+(1-h)] + V_{i+1,j,k+1}[-a-b-f] + V_{i-1,j+1,k+1}[c] + V_{i,j+1,k+1}[-d-e-g] + V_{i+1,j+1,k+1}[-c], (1) $$ where

$a = 0.5\sigma_x^2\dfrac{1}{dx^2}(\dfrac{\partial x}{\partial F})^2d\tau$

$b = 0.5\sigma_x^2\dfrac{1}{2dx}\dfrac{\partial^2 x}{\partial F^2}d\tau$

$c = \rho\sigma_x\sigma_{y}\dfrac{\partial x}{\partial F}\dfrac{\partial y}{\partial \alpha}\dfrac{1}{2dxdy}d\tau$

$d = 0.5\sigma_{y}^2\dfrac{1}{dy^2}(\dfrac{\partial y}{\partial \alpha})^2d\tau$

$e = 0.5\sigma_y^2\dfrac{1}{2dy}\dfrac{\partial^2 y}{\partial \alpha^2}d\tau$

$f = \mu_x\dfrac{\partial x}{\partial F}\dfrac{1}{2dx}d\tau$

$g = \mu_{y}\dfrac{\partial y}{\partial \alpha}\dfrac{1}{2dy}d\tau$

$h = - rdt$

Writing in matrix notation,

$V_{k+1} = A^{-1}( V_{k} - C_{k+1} )$

Note that vector $C$ contains the values that cannot be incorporated via the $A$ matrix, as they depend on boundary grid points.

Edit June 19: Since my other post is focused on the upper boundary in $F$ dimension, lets discuss upper bound in vol direction here, since @Yian_Pap provided an answer below regarding this. Note that I corrected the cross derivative, to be: $\dfrac{\partial^2 V}{\partial F \partial \alpha}=\dfrac{V_{i+1,j+1} - V_{i+1,j-1} - V_{i-1,j+1} + V_{i-1,j-1}}{2\Delta F\Delta \alpha}$, not containing $V_{i,j}$.

Now, as vol bound I set $\dfrac{\partial V}{\partial \alpha}=0$.

Substituting second order accurate backward FD approximation,

$\dfrac{1}{\Delta \alpha}(V_{i,M}-V_{i,M-1})=0$,

since the term in front is not zero, it should hold that,

$V_{i,M}-V_{i,M-1}=0$,

hence,

$V_{i,M}=V_{i,M-1}$ (2),

This can be implemented in the coefficient matrix $A$. Given (1),

$V_{i,j,k} = z_1 V_{i-1,j-1,k+1} + z_2 V_{i,j-1,k+1} + z_3 V_{i+1,j-1,k+1}+ z_4V_{i-1,j,k+1} + z_5V_{i,j,k+1} + z_6V_{i+1,j,k+1}+ z_7V_{i-1,j+1,k+1} + z_8V_{i,j+1,k+1} + z_9V_{i+1,j+1,k+1}$,

Impose the condition (2) as follows:

$V_{i,M,k} = z_1 V_{i-1,M-1,k+1} + z_2 V_{i,M-1,k+1} + z_3 V_{i+1,M-1,k+1}+ (z_4+z_7)V_{i-1,M,k+1} + (z_5+z_8)V_{i,M,k+1} + (z_6+z_9)V_{i+1,M,k+1}$,

Main question: Am I implementing this bound correctly?

Side question: When setting $\nu=0$ and $\beta=1$, $z_7$, $z_8$ and $z_9$ are equal to zero, correct? So, in this case, boundary condition for vol will not affect FD results?

Best,

Pim

Answer 1

The cross derivative FD formula is wrong to start with. If you take out the factor 4 from the denominator, it would become a valid (1st order) formula. As is, it's plain wrong and may well be the reason for your problems (though of course there may be more reasons). By the way, your other spatial FD approximations are second order, so not sure you intended this (a mix of 1st and 2nd order), probably not.

If that doesn't solve your problem and you are still interested, I'll edit and add more. For example, the choice of the cross derivative formula can indeed make the numerical solution more or less prone to negative values, depending on the sign of the correlation coefficient. Using central FD approximations for the convection (1st derivative) terms (as you are doing) can also lead to spurious oscillations and thus negative values, if your problem is convection dominated. That is if the coefficients of the convection terms are large compared to those of the diffusion terms.

EDIT Jun 15: I see that you still haven't corrected the cross derivative discretization, which you must do as it's plain wrong. Now, if you're using $ν=0$ it shouldn't matter, but if you want to get SABR to work in the general case and not only when it reduces to BS, you should correct it! Try the 2nd order formula that doesn't involve the $V_{i,j,k}$ point.

You now have a few separate posts around this problem, so I am a little confused now as to what boundary conditions you are using. Anyway, you don't mention the strike value in your sample calculation, so I have no idea as to how far your $F_{max}$ (250) is from the strike. When you're looking at such huge values of volatility (1000%), you need to truncate your underlying (F) grid very, very far from the strike in order for the solution to have the "space" to diffuse properly. So if your strike is say 50 and you are truncating at 250, then this would be the reason for your inaccurate results at such high vols. Try setting $F_{max}$ to say 50 times the strike and tell me if the situation improves. Of course if you do that, you will at the same time increase the discretization error if you keep using 60 grid points in the F-direction (because you are using a uniform grid and your $dF$ will become much larger). So you'd need to increase the number of F-grid points too possibly. Lots of things to consider:)

Alternatively, a better boundary condition imo for the upper vol boundary is to set $\dfrac{\partial V}{\partial α}=0$ (Neumann). Also, since you are using $ν=0$ at this stage, there is no need for your volatility upper boundary to be that far out. You can truncate at a value just a little higher than the vol you are interested in finding the option price for. In practice the vol will never be 1000% and you would only need to truncate the grid that far out if your $ν$ was non zero and quite large.