Boundary condition issues for Black-Scholes PDE using finite-differences

Question

I have been implementing an, in my opinion, interesting finite difference method (Runge-Kutta-Legendre of second order) to price American options in the standard Black-Scholes model (see "Pricing American Options with the Runge-Kutta-Legendre Finite Difference Scheme, F. Le Floc'h, International Journal of Theoretical and Applied Finance, vol 24, 2021"). This is an explicit method using so called super time-stepping.

My question is not related directly to the method itself though, but rather the choice of boundary conditions for the finite difference grid in the spatial direction. Many authors recommend assuming that the second derivative is zero on the boundary, $\frac{\partial^2 f}{\partial S^2}(t,S_{min}) = \frac{\partial^2 f}{\partial S^2}(t,S_{max}) = 0$.

However, this seems to lead to large errors in several cases, even if the upper boundary $S=S_{max}$ is chosen very large. The motivation given in articles and books seems to be that this condition is true for most payoffs.
While I agree that this is reasonable on maturity, the condition does not seem very appropriate when going backwards in time on the finite difference grid (Unless I made some trivial mistake, of course, which is certainly possible.).

To make things a bit more concrete, lets look at the (time-reversed) Black-Scholes equation: \begin{equation} \frac{\partial f}{\partial t} =\frac{1}{2} \sigma^2 S^2\frac{\partial^2f}{\partial S^2} + \mu S\frac{\partial f}{\partial S} - rf \end{equation} Here $r$ is the risk-free rate, $\mu = r - \delta$, where $\delta$ is the continuous dividend yield and $\sigma$ is the volatility.

Let us focus on the upper boundary. Assume we try to price a standard european put option with maturity in $5$ years, with strike $K=100$, where our upper spatial boundary on the grid is $S_{max} = 500$. For S = S_max, we then have \begin{equation} \frac{\partial f}{\partial t}(t,S_{max}) = \mu S\frac{\partial f}{\partial S}(t,S_{max}) - rf \end{equation} To make this clear, assume that $r=\mu=0$. Then our condition at the upper boundary becomes \begin{equation} \frac{\partial f}{\partial t}(t,S_{max}) = 0, \quad \forall t \end{equation} We have that \begin{equation} f(0, S_{max}) = max(K - S_{max}, 0) = 0. \end{equation} This together means that we will have $f(t, S_{max}) = 0$ for all t. The value gets stuck at 0. This results in large errors, especially for longer maturities and higher volatilities. In fact, we can see that the second order term we assumed to be zero is actually far from zero, especially compared to the other terms in the Black-Scholes PDE (which were actually zero in our case).
More precisely If we look at the second order term at the time we price the option (real non reversed time $t=0$) \begin{equation} \frac{1}{2} \sigma^2 S^2\frac{\partial^2f}{\partial S^2}, \end{equation} we have that this is $\frac{1}{2} \sigma^2 S_{max}^2 \cdot \Gamma$, where $\Gamma$ is the gamma of the option. We can calculate this explicitly for european options. We have \begin{equation} \frac{1}{2} \sigma^2 S_{max}^2 \cdot \Gamma = \frac{\sigma \cdot S_{max}}{2 \sqrt{T}} e^{-d_1^2/2}, \end{equation} where \begin{equation} d_1 = \frac{\log(\frac{S_{max}}{K}) + T \cdot \frac{\sigma^2}{2}}{\sigma \cdot \sqrt{T}} \end{equation} This is not close to zero unless $S_{max}$ is extremely large, so ignoring this term typically results in a substantial error, if I did not make a mistake.

So the question is if I am doing something wrong in my reasoning here. I can see the errors clearly on plots and also the theory seems to confirm that. But since the condition is so commonly recommended I wonder if I might have missed something here?

Is there any easy solution for this problem? One way I found around this is to make lots of points with distance between them growing exponentially outside the original upper boundary. But this makes calculations take more time as well.

Answer 1

Alright, so I solved my issues here. Since it might be useful to others stuck on the same thing, I mention it here. The problem I described about assuming the second order term was equal to zero is real, but as I mentioned, the problem goes away if the upper limit of the grid is huge.

My solution was to not use a uniform grid, but rather to create a uniform grid on some interval $[-C,C]$ for some large enough $C$, say $C=15$. We get grid points $x_i = -C + dx \cdot i, \quad i= 0, 1,2, \ldots$ for $dx = 2C/ (\textrm{number of spatial steps})$.
Then I create the grid points in the asset variable as $S_i = e^{x_i}$ for each gridpoint $x_i$ on the uniform grid. In this way one can reach the endpoint in the upper grid far way (where the second order term is negligible) without using too many steps in the spatial direction.

This is of course not some new magic trick, but rather imitates how the grid would transform after making a change of variable to the log space $x = log(S)$, which is one of the most common way to transform the Black-Scholes equation before doing any numerical calculations. However, I wanted to see if I could do everyting in the original formulation without using this transform. So it turned out fine, after using the "exponential" grid described. The problem was due to me using a uniform grid.