Black–Scholes equation

Suppose that at time $t=0$ you buy a stock whose share price is $S(t)$. At a later time, if $S(t)>S(0)$, you can sell the stock and make money. But if $S(t)<S(0)$, you stand to lose money—potentially, your entire investment. You may prefer to mitigate this risk.

One way to do so is to buy a call option instead of the stock. This is a contract with a fixed strike time $T$ and strike price $K$ that gives you the right to buy the stock from the contract issuer at cost $K$ at time $T$.^[1]

Suppose $S(T)>K$. Then you could buy the stock at price $K$ and instantly resell it at price $S(T)$, so you make profit $S(T)-K$. On the other hand, if $S(T)\le K$, there is no advantage to exercising the option because the stock is less valuable than the guaranteed purchase price. However, you have lost only what you paid for the option. These observations are summarized by the payoff function

What is the fair market value of the option contract? This question is answered to close approximation by the famous Black–Scholes equation,

where $r$ is the risk-free interest rate (i.e., what you could earn with a very safe investment), and σ is the volatility of the stock (essentially, the standard deviation in the rate of return for the stock). In the Black–Scholes model, both $r$ and σ are assumed to be known and constant.

The value $v(S,t)$ of the option depends on the time $t$ and on the stock price $S$, which may be varied independently. At the strike time, the payoff condition $v(S,T)=H(S)$ is imposed, and the goal is to solve the equation backward in time to find $v(S,0)$. Equation (11.1.2) is a time-dependent partial differential equation, or evolutionary PDE.

11.1.1Initial and boundary conditions¶

Equation (11.1.2) has independent variables $t$ and $S$. Because of the first-order derivative in time, we require an initial value for the dependent variable.

In order to follow the usual convention of having time flow forward instead of backward, we define a new variable $\eta=T-t$, in which case (11.1.2) becomes

now defined for $0\le \eta \le T$. Here we have adopted the common notation of using subscripts for partial derivatives. In the new time variable we have the initial condition

and the goal is to find $v(S,\eta)$ for $\eta>0$. Henceforth we will simply rename the η variable as $t$ again, with the understanding that it runs in the opposite direction to actual time in this application.

We have not yet specified the domain of the independent variable $S$. Stock prices cannot be negative, so $S= 0$ is the left end. In principle, $S$ is not bounded from above. There are different ways to handle that fact computationally, but here we choose to truncate the domain at some positive value $S_\text{max}$.

Because (11.1.2) has two derivatives in $S$, we also require a condition at each end of its domain. At $S=0$, both the stock and the option are worthless, so we impose $v=0$ there. In light of the payoff function $H$, initially the value $v$ has slope equal to 1 at $S=S_\text{max}$, and we choose to enforce that condition for all time. We summarize the boundary conditions as

At the left we have a homogeneous Dirichlet condition, and at the right we have a nonhomogeneous Neumann condition.

The complete problem consists of (11.1.3), (11.1.4), and (11.1.5). An evolutionary PDE has both initial and boundary conditions and is sometimes called an initial-boundary-value problem or IBVP as a result. Solving it numerically requires aspects of methods for both IVPs and BVPs.

11.1.2The heat equation¶

The Black–Scholes equation can be transformed by a change of variables into a simpler canonical PDE.^[2]

The heat equation is the archetype differential equation for the class known as parabolic PDEs. A diffusive process is one in which velocity is proportional to the gradient of the solution. Thus, rapid changes in the solution flatten out quickly.

In certain cases there are explicit formulas describing the solution of (11.1.6) and hence (11.1.2). These provide a lot of insight and can lead to practical solutions. But the exact solutions are fragile in the sense that minor changes to the model can make them impractical or invalid. By contrast, numerical methods are often more robust to such changes.

11.1.3A naive numerical solution¶

Let’s return to the Black–Scholes equation (11.1.3), rewritten here as

subject to the initial condition (11.1.4) and the boundary conditions (11.1.5). We now try to solve it numerically, without any transformation tricks or other insight. Define

Observe that we have used the more conventional $x$ and $t$ in place of $S$ and η. The result is a grid function $\mathbf{V}$ whose entries are $V_{ij} \approx v(x_i,t_j)$. By the initial condition (11.1.4), we set $V_{i,0} = H(x_i)$ for all $i$.

For the moment, let us pretend that $i$ is unbounded in both directions. Replacing the derivatives in (11.1.9) with some simple finite-difference formulas, we get

We can rearrange (11.1.11) into

where

If we put $j=0$ into (11.1.12), then everything on the right-hand side is known for each value of $i$, and thus we can get values for $V_{i,1}$. We can then use $j=1$ and get all of $V_{i,2}$, and so on.

Now we take the boundaries on $x$ into account. The value $V_{0,j+1}$ is zero, so we don’t even need to compute the solution using (11.1.12) at $i=0$. At $i=m$, or $x_i=S_\text{max}$, things are trickier. We don’t know $V_{m,j+1}$ explicitly and would like to solve for it, yet formula (11.1.12) refers to the fictitious value $V_{m+1,j}$ when $i=m$. This is where the Neumann condition must be applied. If the value $V_{m+1,j}$ did exist, then we could discretize that condition as

We can therefore solve (11.1.14) for the fictitious $V_{m+1,j}$ and use it where called for in the right-hand side of (11.1.12).

Everything in Demo 11.1.2 seems to go smoothly. However, trouble lurks just around the corner.

The explosive growth of error in Demo 11.1.3 suggests that there is instability at work. Understanding the source of that instability comes later in this chapter. First, though, we consider a general and robust strategy for solving evolutionary PDEs.

11.1.4Exercises¶

1. ✍ Show that $u(x,t) = e^{-4 \pi^2 t} \cos (2 \pi x)$ is a solution to the IBVP

   $$\begin{align*} \text{PDE:} \quad & \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2},\quad 0<x<1,\; t > 0, \\ \text{BC:} \quad & \frac{\partial u}{\partial x}(0,t) = \frac{\partial u}{\partial x}(1,t) = 0,\quad t > 0, \\ \text{IC:} \quad & u(x,0) = \cos (2 \pi x),\quad 0 \le x \le 1. \end{align*}$$

   (11.1.15)

2. ✍ Show that $u(x,t) = t^{-1/2} \exp(-x^2/4t)$ solves the heat equation (11.1.6) at any value of $t>0$.

3. ✍ Equation (11.1.11) results from applying finite differences to the derivatives in (11.1.9), including a forward difference for the term $v_t$.

   (a) Write out the method that results if a backward difference is used for $v_t$ instead.

   (b) Explain why modifying the code from Demo 11.1.2 to implement this formula requires the use of matrix algebra.

4. ⌨ In this problem you are asked to revisit Demo 11.1.2 in order to examine the instability phenomenon more closely.

   (a) Leaving other parameters alone, let $m=100$. To the nearest ten, find the minimum value of $n$ that leads to a stable (i.e., not exponentially growing) solution.

   (b) Repeat (a) for $m=120,140,\ldots,200$. Make a table of the minimum stable $n$ for each $m$. Is the relationship $n=O(m)$, or something else?