Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c) 2014 Tingyu Wang.

Black-Scholes Equation

In the class, we discussed how to use numerical methods to solve wave propagation problems and heat equation, which both represent a physical phenomenon. In fact, this powerful tool is not only widely applied by engineers and physicists, but also by traders in a stock market and economists working for the government. Economists would build their own financial models in PDE to estimate future economic growth for each year, and PDE solvers can give traders a sketch about the price of a stock or an option.

Black-Scholes equation is one of the most important models concerning the option pricing. It was first published by Fischer Black and Myron Scholes in 1973. I first read this equation through a book called "17 equations that changed the world".

What is Black-Scholes Equation?

Before we introduce the equation, let us start with some important terminologies behind the equation. Options are used on markets and exchanges. From the definition in Wikipedia, an option is a contract which gives the buyer (the owner) the right, but not the obligation, to buy or sell an underlying asset at a specified price on or before a specified date. We call the price determined in the contract "exercise price (or strike price) ($E$), and call the date "the expiry time ($T$)". The underlying asset has its own value, and we denote the current value of underlying asset as ($S$). If the asset is a car, then $S$ may be equal to the present value of the car, say 20,000 dollars. If the asset is a stock, then $S$ is equal to the current price of the stock (find it in Nasdaq, Dow or S&P).

From the definition, we know that the owner of the option can buy or sell some amount of assets on or before the expiry time. We call it an European option if the option can only be executed on the expiry time ($t_{exercise} = T$). If it is an American option, than means it can be exercised at any time before the expiry time($t_{exercise} \leq T$).

We can understand the option as a right or a priviledge, so the option itself should also have a price, "the price of an option ($V$)", which depends on the current value of underlying asset ($S$) and time ($t$): $V = V(S,t)$. Notice that there is a "buy or sell" in the definition of option above. If the owner of the option has the right to buy some shares of Apple's stock at 50 dollars per share from the buyer before Feb 1, 2015, then this type of option is called call option. We use $C(S,t)$ instead of $V(S,t)$. If the owner of option has the right to "sell", then this type of option is called put option. We denote it as $P(S,t)$ instead of $V(S,t)$. Click the short video below for real-world examples of call options and put options.

Finally, we introduce the Black-Scholes equation as follows.

$$\begin{equation} \frac{\partial V}{\partial t}(S, t) + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}(S, t) + rS \frac{\partial V}{\partial S}(S, t) - rV(S, t) = 0 \end{equation}$$

$V$ stands for the price of the option, $S$ denotes current value of underlying assets, $t$ is time, $r$ denotes the no-risk interest rate (ex.saving rate in bank), and $\sigma$ denotes the volatility of underlying asset (ex.the standard deviation of the stock's returns).

As you can see, the Black-Scholes equation is aimed at calculating the price of an option, so it is the most popular option-pricing method. The only unknown here is the price of option ($V$), depending on $S$ and $t$. In the following part, we will use finite difference methods (both explicit and implicit schemes) to solve this nonlinear second-order parabolic equation.

Initial Conditions

Now, consider we are holding an European, call option $C(S,t)$ of a stock and we are at the expiry time: $t=T$.

• If the current price of the stock ($S$) > exercise price ($E$), then we definitely want to execute our right to buy the stock at the price $E$, and the profit we gain from the option is $S-E$.

• If the current price of the stock ($S$) < exercise price ($E$), then there is no need to exercise the call option, therefore, in this case the value of this option is zero.

Thus, for a call option,

$$C(S,T) = max(S-E,0)$$

This is the final condition for the Black-Scholes equation, however, usually we prefer to deal with an inital condition with a form of $C(S,0) = f(S)$. So here we apply change of variables:

$$ t^*(t) = T - t$$

then,

$$\frac{\partial V(S,t)}{\partial t} = - \frac{\partial V(S,t^*)}{\partial t^*}$$

substitute in the Black-Scholes equation:

$$-\frac{\partial V}{\partial t^*}(S, t^*) + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}(S, t^*) + rS \frac{\partial V}{\partial S}(S, t^*) - rV(S, t^*) = 0 $$

changing the sign, and use $C$ instead of $V$ since it is a call option:

$$\frac{\partial C}{\partial t^*}(S, t^*) - \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 C}{\partial S^2}(S, t^*) - rS \frac{\partial C}{\partial S}(S, t^*) + rC(S, t^*) = 0 $$

with an initial condition: $C(S,0) = max(S-E,0)$. Be careful that this initial condition reflects the expiry moment $t = T$, and when we are stepping forward to calculate the last time step in $t^*$, we are actually calculating the price of the call option at $t =0$ (that is when we are buying an option).

Boundary Conditions

In order to make it a well-posed partial differential equation, we need to define the boundary conditions on both sides. The current price of the stock $S$ can take value in the interval $[0,\infty)$, which is unbounded.

• left boundary: $S=0$, then $C(0,t^*) = 0$
• right boundary: $S\rightarrow \infty$, then $C(S,t^*) = (S-E) $

However, in the numerical method, we need to set an artificial bound $S_{max}$ concerning the space discretization, otherwise we will have inifinity number of grids. Based on Dura's work [ref 1], it is recommended to set $S_{max}$ around four times the exercise price $E$. Here we assume that

$$S_{max} = 4E $$

thus, the boundary conditions can be written as follows:

$$\begin{cases} C(0,t^*) = 0\\ C(S_{max},t^*) = S_{max} - E\\ \end{cases}$$

Problem Statement

Now let's solve a practical problem using our numerical skills!

Consider an European call style option is made for a stock currently trading at 20 dollar per share with volatility 0.4. The term is 3 months (0.25 year) and the strike price is 10 dollar. The prevailing no-risk interest rate is 10 %. What should the price per share be for the option?

Explicit Scheme

For the explicit method, we use a forward difference scheme in time $t$, and a central difference scheme in the current price of underlying asset $S$. The truncation error here is second-order in space and first-order in time $O(\Delta t, \Delta S^2)$. Write out discretized form for each term:

$$\begin{equation} \frac{\partial C}{\partial t^*} = \frac{C_i^{n+1} - C_i^{n}}{\Delta t}\\ \frac{\partial C}{\partial S} = \frac{C_{i+1}^{n} - C_{i-1}^{n}}{2\Delta S}\\ \frac{\partial^2 C}{\partial S^2} = \frac{C_{i+1}^{n} - 2C_i^{n} + C_{i-1}^{n}}{(\Delta S)^2}\\ \end{equation}$$

the discretized Black-Scholes equation: $$\frac{C_i^{n+1} - C_i^{n}}{\Delta t} - \frac{1}{2} \sigma^2 S^2 \frac{C_{i+1}^{n} - 2C_i^{n} + C_{i-1}^{n}}{(\Delta S)^2} - rS \frac{C_{i+1}^{n} - C_{i-1}^{n}}{2\Delta S} + rC_i^{n} = 0 $$

since $S=i \Delta S$, we can rearrange the above equation as:

$$\frac{C_i^{n+1} - C_i^{n}}{\Delta t} - \frac{1}{2} \sigma^2 i^2 (C_{i+1}^{n} - 2C_i^{n} + C_{i-1}^{n}) - \frac{1}{2} ri (C_{i+1}^{n} - C_{i-1}^{n}) + rC_i^{n} = 0 $$

the option price at the next time step can be expressed as: $$\begin{equation} C_i^{n+1} = \frac{1}{2}(\sigma^2 i^2 \Delta t-ri\Delta t)C_{i-1}^n + (1 - \sigma^2 i^2 \Delta t -r \Delta t) C_i^n + \frac{1}{2}(\sigma^2 i^2 \Delta t+ri\Delta t)C_{i+1}^n \end{equation}$$

Keep in mind that we are using an explicit scheme, so start with a small time step to satisfy the stability requirement.

Run the code and calculate the price of the call option when the stock is currently trading at 20 dollars. (See problem statement)

Let us plot the call option price against the current stock price based on our result.

Implicit Scheme

As we learned in class that the explicit schemes contraint the size of time step, and they will become unstable if we use a larger $\Delta t$. That is why we want to solve it in an implicit way that is unconditionally stable. However, the fully implicit scheme only gives us the fisrt-order accurary in time. Remember in the fourth module, we introduced the Crank-Nicolson scheme to solve the heat equation. The Black-Scholes equation and heat equation are both parabolic equations. So why not try to apply Crank-Nicolson scheme in order to achieve a better accurary?

$$\frac{\partial C}{\partial t^*}(S, t^*) - \frac{1}{2} \sigma^2 S^2 \frac{\partial^2 C}{\partial S^2}(S, t^*) - rS \frac{\partial C}{\partial S}(S, t^*) + rC(S, t^*) = 0 $$

Recall that Crank-Nicolson method is an average of implicit scheme and explicit scheme. Based on the discretized form of Black-Scholes equation, we rewrite each partial derivative in space as an average of the schemes in time step $n$ and $n+1$.

$$\frac{C_i^{n+1} - C_i^{n}}{\Delta t} - \frac{1}{4} \sigma^2 S^2 (\frac{C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}}{(\Delta S)^2} + \frac{C_{i+1}^{n} - 2C_i^{n} + C_{i-1}^{n}}{(\Delta S)^2}) -\\ \frac{1}{2} rS (\frac{C_{i+1}^{n+1} - C_{i-1}^{n+1}}{2\Delta S} + \frac{C_{i+1}^{n} - C_{i-1}^{n}}{2\Delta S}) + \frac{1}{2} r (C_i^{n+1}+C_i^{n} )= 0 $$

replace $S$ with $i\Delta S$, and rearrange terms so that all the unknowns (values in $n+1$ time step) are on the LHS:

$$\frac{\Delta t}{4} (ri-\sigma^2i^2) C_{i-1}^{n+1} + (1+\frac{\Delta t}{2}(r+\sigma^2i^2))C_i^{n+1} + (-\frac{\Delta t}{4}(ri+\sigma^2i^2))C_{i+1}^{n+1}\\= \frac{\Delta t}{4} (\sigma^2i^2-ri) C_{i-1}^{n} + (1-\frac{\Delta t}{2}(r+\sigma^2i^2))C_i^{n} + (\frac{\Delta t}{4}(ri+\sigma^2i^2))C_{i+1}^{n}$$

for notation convenience, we define:

$$\alpha_i = \frac{\Delta t}{4} (ri-\sigma^2i^2) \\ \beta_i = \frac{\Delta t}{2}(r+\sigma^2i^2)\\ \gamma_i = -\frac{\Delta t}{4}(ri+\sigma^2i^2)$$

for $i = 1,2,3, ..., M-1$

so we can simplify the discretized equation as follows:

$$\alpha_i C_{i-1}^{n+1} + (1+\beta_i)C_i^{n+1} + \gamma_i C_{i+1}^{n+1} = -\alpha_i C_{i-1}^{n} + (1-\beta_i)C_i^{n} - \gamma_iC_{i+1}^{n}$$

This linear system of equations can be written as:

$$ AC^{n+1} = A_*C^{n} + B $$

A is the LHS coefficient matrix: \begin{align}\left[ \begin{array}{cccccc} \left(1+\beta_1 \right) & \gamma_1 & 0 & \cdots & 0 \\ \alpha_2 & \left(1+\beta_2 \right) & \gamma_2 & & \vdots \\ 0 & \alpha_3 & \ddots& \ddots & 0 \\ \vdots & & \ddots & \left(1 + \beta_{M-2}\right) & \gamma_{M-2} \\ 0 & \cdots & 0 & \alpha_{M-1} & \left(1 + \beta_{M-1}\right) \end{array} \right] \end{align}

At the RHS of the equation, $A_*$ is the coefficient matrix of $C^{n}$ (the call option price array at the previous time step $n$). From above we know matrix $A_*$ is very similar to matrix $A$:

\begin{align}\left[ \begin{array}{cccccc} \left(1-\beta_1 \right) & -\gamma_1 & 0 & \cdots & 0 \\ -\alpha_2 & \left(1-\beta_2 \right) & -\gamma_2 & & \vdots \\ 0 & -\alpha_3 & \ddots& \ddots & 0 \\ \vdots & & \ddots & \left(1 - \beta_{M-2}\right) & -\gamma_{M-2} \\ 0 & \cdots & 0 & -\alpha_{M-1} & \left(1 - \beta_{M-1}\right) \end{array} \right] \end{align}

Vector $B$ is the boundary condition vector based on the Dirichlet boundary condition: $C(0,t) = 0, C(S_{max},t) = S_{max} -E$. Consider the total number of spatial grids $M = 200$ (the 201 grid points' indices are from 0 to 200), then the last equation in the linear system should be:

$$\alpha_{199}C_{198}^{n+1} + (1+\beta_{199})C_{199}^{n+1} + \gamma_{199}C_{200}^{n+1} \\ = -\alpha_{199}C_{198}^{n} + (1-\beta_{199})C_{199}^{n} - \gamma_{199}C_{200}^{n}$$

and we know from the BC that $C_{200}^{n} = C_{200}^{n+1} = (S_{max}-E)$, so rearrange the equation as： $$\alpha_{199}C_{198}^{n+1} + (1+\beta_{199})C_{199}^{n+1} \\ = -\alpha_{199}C_{198}^{n} + (1-\beta_{199})C_{199}^{n} - 2\cdot \gamma_{199}(S_{max}-E)$$

and since the BC at the lower bound is $C(0,t) = 0$, we don't need add any term in the first equation. Thus, the boundary condition vector $B$ should look like:

\begin{align}\left[ \begin{array}{cccccc} 0 \\ \vdots \\ 0 \\ -2\gamma_{199}(S_{max}-E) \end{array} \right] \end{align}

Knowing the value of matrix $A$ and vector $b$, we can solve the linear system using the numpy.linalg.solve function.

We first set the initial condition and boundary condition, and calculate the $\alpha$, $\beta$, $\gamma$ coefficient arrays. Since the Crank-Nicolson scheme is unconditionally stable, we can choose a larger value for time step $dt$.

Given the coefficient arrays $\alpha$, $\beta$, $\gamma$, we can produce the LHS matrix $A$, update the RHS vector $b$ at each time step, and solve the linear system again. Finally we have the solution by Crank-Nicolson scheme.

Plot the result and compare with the solution of explicit scheme (FTCS).

These two solutions are very close to each other. Are they both correct? Let us compare them with an exact solution.

An Exact Solution of Black-Scholes Equation

As for an European call option problem, mathematicians provide us with an analytical solution (here is the derivation of the solution):

$$\begin{equation} C(S,t) = SN(d_1) -Ee^{-r(T-t)}N(d_2) \end{equation}$$

where $N(x)$ is the cumulative distribution function for a standardised normal random variable (we use norm.cdf to calculate it in python), given by

$$\begin{equation} N(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} {e^{-\frac{1}{2}y^2}dy} \end{equation}$$

and

$$\begin{equation}\begin{split} d_1 & = \frac{ln(S/E)+(r + \frac{1}{2}\sigma^2)T}{\sigma \sqrt{(T-t)}}\\ d_2 & = d_1 - \sigma \sqrt{(T-t)} \end{split}\end{equation}$$

Recall that we change variable $t^* = T-t$ when we are dealing with the initial condition, so the analytical solution should be:

$$\begin{equation} C(S,t^*) = SN(d_1) -Ee^{-rT}N(d_2) \end{equation}$$

and

$$\begin{equation}\begin{split} d_1 & = \frac{ln(S/E)+(r + \frac{1}{2}\sigma^2)T}{\sigma \sqrt{t*}}\\ d_2 & = d_1 - \sigma \sqrt{t^*} \end{split}\end{equation}$$

Don't forget there is still an underlying condition that when $S<E$, we won't exercise the call option. Therefore, the price of call option should be zero rather than a negative value. This explains the if condition inside the loop.

Let plot the exact solution together with our numerical results.

Provided the fact that the stock is currently trading at 20 dollars, our pricing strategy of the call option is shown below. The result of Crank-Nicolson method is slightly closer to analytical result.

Using the parameters: $r = 0.1$, $\sigma = 0.4$, $E = 10$, we can generate a pricing strategy by changing current the stock price (S) and expiry time (T) as the table below (given by Crank-Nicolson scheme):

Reference

1. Dura, Gina, and Ana-Maria Moşneagu. "Numerical approximation of Black-Scholes equation." Annals of the Alexandru Ioan Cuza University-Mathematics 56.1 (2010): 39-64.
2. Stewart, Ian. In Pursuit of the Unknown: 17 Equations That Changed the World Paperback, Basic Books, 2013. Print
3. "Option (finance)." Wikipedia. Wikimedia Foundation, n.d. Web. 13 Jan. 2015.
4. Fischer Black and Myron Scholes (1973). The Pricing of Options and Corporate Liabilities, Journal of Political Economy, Vol. 81, No. 3 (May - Jun., 1973), pp. 637-654
5. Francois Coppex. Solving the Black-Scholes equation: a demystification. 2009. http://www.francoiscoppex.com/blackscholes.pdf
6. Jiang, L. (2005). Mathematical modeling and methods of option pricing. Singapore: World Scientific.