options - SABR PDE spot/forward upper boundary condition implementation

When running my Finite Difference code, I observe something odd.

Although implementing a classical (non-reverting) SABR model, I initialized the variables such that it should be equal to Black-Scholes.

Boundary conditions for

• 
  
  

  Lower bound on price

  
  


• 
  

  Lower bound on volatility

  
  


• 
  

  Upper bound on volatility

  
  

are Dirichlet-style (just set equal to a value).

For the upper bound on forward price ($F$), I want to set the Neumann condition of $\Gamma=0$, since I believe this is true for all (put and call) options for large underlying price. This condition yields:

$\dfrac{\partial^2 V}{\partial x^2}(\dfrac{\partial x}{\partial F})^2+\dfrac{\partial V}{\partial x}\dfrac{\partial^2 x}{\partial F^2}=\Gamma=0$

Note that the extra terms with $x$ and $F$ are due to my variable transformation from $F$ to $x$.

Edit June 19: @Yian_Pap I agree, lets forget about smoothing condition and focus on getting Neumann condition working. Let me be specific regarding my steps applying your approach :) This is how I would currently implement it, but once again I'm not very familiar with FD (this is my first use case):

• 
  

  First of all, it is not immediately clear to me what PDE to start with (pricing equation? or just the $\Gamma=0$ equation?). If we start from $\Gamma=0$, why would we bother about cross derivative term being 0 or not?

  
  


• 
  

  If I start with $\Gamma=0$, no transformations in F, then

  
  
  

$\dfrac{\partial^2 V}{\partial F^2} = 0$

discretize using 2nd order central FD approx,

$\dfrac{1}{\Delta F^2} (V_{N-1,j} - 2V_{N,j} + V_{N+1,j}) = 0$,

$V_{N,j}=\dfrac{1}{2}(V_{N-1,j}+V_{N+1,j})$ (1),

which contains the outer grid point $V_{N+1,j}$, which I assume you mean in your answer?

Then, as I understand your comment regarding that the first derivative at $F \to \infty$ is known and constant, you mean that the delta is 1 (0) for a call (put)?

$\dfrac{\partial V}{\partial F} = \Delta = 1$,

$V_{N,j} = \Delta F + V_{N-1,j}$,

use this to substitute $V_{N+1,j}$ in (1),

$V_{N,j} = \Delta F + V_{N-1,j}$,

with $\Delta F$ step size in F. Then I can incorporate this in the coefficient matrix as follows. (Side question: why not use $\Delta=1$ condition directly?)

Say I have implicit scheme $V_{i,j,k} = z_1V_{i-1,j,k+1} + z_2V_{i,j,k+1} + z_3V_{i+1,j,k+1}$. Then for some i (i.e. at some places in matrix) $V_{i+1,j}$ is a bound, hence matrix entry is set $0$. Incorporate bound by adding $z_3$ to $z_2$ at every place in matrix where $i=N-1$. Additionally add $z_3\Delta F$. Mathematically,

$V_{N-1,j,k} = z_1V_{N-2,j,k+1} + (z_2+z_3)V_{N-1,j,k+1} + z_3\Delta F$.

Would this be a correct understanding of how the upper bound F Neumann condition should be implemented?

Answer

Wherever your discretized PDE references $V_{N,j,k+1}$, you will use your expression for $V_{N,j,k+1}$ in terms of $V_{N−1,j,k+1}$ and $V_{N−2,j,k+1}$ to eliminate $V_{N,j,k+1}$ from the discretized equation. This means that the algebraic system you will solve will not involve any $V_{N,j,k+1}$ values directly (so no $C_{k+1}$ contribution from the F-discretization). Rather you will be augmenting the system coefficients of the $V_{N−1,j,k+1}$ and $V_{N−2,j,k+1}$ points accordingly.

Edit Jun 18. OK, let's start again and assume you are using the pricing PDE in the forward F (no transformations). If you use a smoothing condition you'd need to use one-sided finite difference approximations so that you can close your algebraic system at the last line $F=F_{max}$. That said, the reason why you opted for this kind of (smoothing) condition is because when you tried with Neumann, you didn't get correct results. But this was most probably because you made an error somewhere. I would try again as follows for the last line of points of your grid: the first derivative in F (let's call it $D_{F\to\infty}$) you already know (you know its value from the Neumann condition), so this doesn't need to be discretized on the boundary, you just use its known value. The cross derivative also becomes zero there since the first derivative in the F-direction is constant (even in the general case where $\nu$ is not zero). The only derivative in the F-direction that you'll need to discretize is the second derivative and that will reference a point outside your grid (with the standard 2nd order central FD approximation). This "fictitious" point in your discretized equations you will replace by its extrapolated value as $V_{N+1,j,k+1} = V_{N,j,k+1} + D_{F\to\infty} d_F $. This will again close your system and you're good to go. Try and let me know if something is still not clear.