Errors on Finite Differences + Implicit Scheme + Black & Scholes

Question

I'm solving the classical Black & Scholes (BS) PDE for a European option using finite difference and the implicit scheme. In other words, I'm trying to solve

$\displaystyle\frac{\partial V}{\partial t} + \frac 12\sigma^2\frac{\partial^2V}{\partial x^2} + \mu\frac{\partial V}{\partial x} - rV = 0$

with boundary condition $V(T,x) = \Phi(e^x) = \max\{e^x - K, 0\}$. The drift parameter is as usual $\mu = r - y - \frac 12 \sigma^2$.

The problem is when my $x$ grid is too narrow, my option prices diverge from the BS. Specifically, I'm solving the ATM case with S0 = K = 2775 and if I use [np.log(2675), np.log(2875)] as the boundary points, the solution has a huge error (~400). However, when I use a wider grid [np.log(1775), np.log(3775)], it works like a charm.

I thought the implicit scheme would be numerically stable enough to not exhibit this kind of behavior, so should I assume there is something wrong with my Python code?

My other parameters are r, y, T, sigma = (0.0278, 0.0189, 1, 0.15).

def buildDs(nx):
    D1 = np.diag(np.ones(nx), 1) - np.diag(np.ones(nx), -1)
    D2 = np.diag(np.ones(nx), 1) + np.diag(np.ones(nx), -1) - 2 * np.eye(nx + 1)

    D1[0, 2] = -1
    D1[0, 1] = 4
    D1[0, 0] = -3

    D1[-1, -1] = 3
    D1[-1, -2] = -4
    D1[-1, -3] = 1


    D2[0, 2] = 1
    D2[0, 1] = -2
    D2[0, 0] = 1

    D2[-1, -1] = 1
    D2[-1, -2] = -2
    D2[-1, -3] = 1

    return csc_matrix(D1), csc_matrix(D2)

def pxFD(Phi, r, mu, sigma,
         xs, ts,
         nx = 2000, nt = 10000):

    dx = (xs[1] - xs[0]) / nx
    dt = (ts[1] - ts[0]) / nt

    xs = np.linspace(xs[0], xs[1], num = nx + 1, endpoint = True)
    ts = np.linspace(ts[0], ts[1], num = nt + 1, endpoint = True)

    V = np.zeros((nt + 1, nx + 1))
    I = identity(nx + 1)

    D1, D2 = buildDs(nx)
    L = 1 / 2 * (sigma / dx) ** 2 * D2 + mu / (2 * dx) * D1 - r * I
    P = I - dt * L

    V[-1] = Phi(xs)
    for j in reversed(range(nt)):
        V[j] = spsolve(P, V[j + 1])    
    return V, xs, ts

def pxFDGBM(Phi, r, y, sigma,
            Ss, ts,
            nx = 2000, nt = 10000):

    mu = r - y - 1 / 2 * sigma ** 2
    xs = np.log(Ss)
    Philog = lambda x: Phi(np.exp(x))

    Vs, xs, ts = pxFD(Phi = Philog,
                      r = r,
                      mu = mu,
                      sigma = sigma,
                      xs = xs,
                      ts = ts,
                      nx = nx,
                      nt = nt)
    return Vs, np.exp(xs), ts

Answer 1

The PDE is defined for $x \in ]-\infty, +\infty[$ but the finite difference scheme requires a truncated domain $[x_{\min}, x_{\max}]$, and the choice of $x_{\min}$ and $x_{\max}$ will affect the quality of the result, regardless of the scheme being explicit, implicit, or mixed.

A good rule of thumb is to choose the truncation $[x_{\min}, x_{\max}]$ such that the probability of $X_t < x_{\min}$ or $X_t > x_{\max}$ is very low under the original SDE, because these are states that won't have much impact on the option value.

In the BS case a truncation at $\pm 5$ standard deviation is usually a good choice. In you example that would be $x_{\min} = \ln(S_0) - 5 \sigma \sqrt{T} \approx \ln(1311)$ and $x_{\max} = \ln(S_0) + 5 \sigma \sqrt{T} \approx \ln(5875)$.

Another thing that can improve accuracy is the choice of the boundary condition at $x_{\min}$ and $x_{\max}$.