Numerical Solution of the Black Scholes equation for European Options

Under some assumptions, the Black-Scholes equation is satisfied by any derivative security whose price depends only on the current value of the stock price $S$ and time $t$. These assumptions and the derivation of the Black-Scholes equation can be found with ease, either online or by referring to a great reference such as The Mathematics of Financial Derivatives [1].

In this article we show that the Black-Scholes equation for European call options (more generally derivative securities) can be transformed into the heat equation, which makes a numerical solution very simple. In a later blog post we will see how to numerically solve problems in the Black-Scholes framework when American options are considered. If fact, the purpose of this post is to get used of the basics before moving on to some more exotic techniques. Also, it gives me a chance to develop test cases for my attempt at a modern C++ linear solver SolvAnt.

By following [1], we can write the Black-Scholes equation for a call option $C(S,t)$ as
$$ \frac{\partial C}{\partial t} + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 C}{\partial S^2} + rS \frac{\partial C}{\partial S} - rC = 0, $$
where $\sigma$ is the volatility and $r$ is the interest rate.
In order to obtain a unique solution we need an "end" condition, given by
$$C(S,T) = max(S-E, 0),  S > 0, $$
where T is the expiry date and $E$ is the exercise price. This essentially says that if the price of the underlying asset is greater than the exercise price then the value of the option is positive. However if the asset price is less than the exercise then the call option is worthless.

We also need "boundary" conditions, given by
$$C(0,t) = 0,  t > 0, $$
$$C(S,t) \rightarrow S - Ee^{-r(T-t)},  S \rightarrow \infty, t > 0$$
In order to understand the first of these, we just need to realise that the value of the underlying asset follows a lognormal random walk ($dS \sim S$). That is, if $S = 0$ the asset price does not change and so the "right to buy" such an asset is worthless. The second of these says that as the asset price gets large the value of the call option approaches that of the asset.

In this form we can most certainly apply the finite difference formulation and solve right away, and in some cases this is what we have to do. By using a change of variables,
$$ S = Ee^x, t = T - \frac{2\tau}{\sigma}, C = Ee^{-\frac{1}{2}(k-1)x-\frac{1}{4}(k+1)^2\tau}u(x,\tau),$$
we obtain the following non-dimensional problem
$$ \frac{\partial u}{\partial \tau} = \frac{\partial^2 u}{\partial x^2}, -\infty < x < \infty,  \tau > 0$$
$$u(x,0) = max(e^{ \frac{1}{2}(k+1)x} - e^{ \frac{1}{2}(k-1)x}, 0),$$
$$u(-\infty, \tau) = 0,$$
$$u(\infty, \tau) =  e^{\frac{1}{2}(k+1)x+\frac{1}{4}(k+1)^2\tau} -  e^{\frac{1}{2}(k-1)x+\frac{1}{4}(k-1)^2\tau}$$

This is similar to a physical problem of solving the heat equation on an infinite bar. What makes this problem interesting to solve numerically is that our domain is infinite in the $x$ direction. In order to overcome this problem we can restrict it to some small sufficiently large neighbourhood $[x_-, x_+]$. We can then evaluate the boundary conditions at the end points of the interval. Now, let $\mathbf{M}$ be a matrix of appropriate size with $-2$ on the diagonal and 1 on the super-diagonal and sub-diagonal, where $N$ is the number of spatial nodes. We can discretise the problem in space using central differences to obtain
$$\frac{d \mathbf{u}}{d t} = \frac{1}{2 \delta x}\mathbf{M} \mathbf{u} + \mathbf{b},$$
where $ \mathbf{b}$ contains the boundary information. We now need to discretise in time. If we use the Crank-Nicholson method we obtain the following matrix vector equation
$$(\mathbb{I} - \alpha \mathbf{M})\mathbf{u}_{n+1} = (\mathbb{I} + \alpha \mathbf{M})\mathbf{u}_{n} + \mathbf{c},$$
where $\alpha = \frac{\delta t}{2 \delta x}$.
We solve this each time step to obtain the option price for that given time. We can use SolvAnt to obtain the solution to the European call option using the following class

namespace solvant {
namespace finance {
template <std::size_t N, std::size_t M>
class EuroCallOptionSolver {
public:
    EuroCallOptionSolver();
    void init(float E, float T, float sig, float r,
              solvant::solver::DirectTriSolver<float, N - 2>* solver);
    void solve();
    void writeSolutionToFile(std::string filename);

private:
    bool m_init;  // has the solver been init'd?

    // Financial variables
    float m_E;    // Exercise price
    float m_T;    // Expiry date
    float m_sig;  // Volitility
    float m_r;    // interest rate

    // non-dimensional parameters
    float m_k;

    float m_dt;
    float m_dx;
    float m_alpha;

    solvant::solver::DirectTriSolver<float, N - 2>* p_solver;
    // pImpl -< can use differnt direct trisolver
    std::array<std::array<float, N - 2>, M> m_sols;
};
}  // namespace finance
}  // namespace solvant

The user of this class only has to provide the solver (making use of the bridge pattern) and then fire it off with the appropriate method. The solver in this case is provided by solvant, which makes use of move semantics to ensure that data is not unnecessarily duplicated.