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".
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).
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.
In [9]:
from IPython.display import YouTubeVideo
YouTubeVideo('EfmTWu2yn5Q')
Out[9]:
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 we 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.
Buyers and sellers use financial derivatives (derivative here means a contract of some underlying assets) to offset risk in their portfolios. The key idea behind the model is to hedge the option by buying and selling the underlying asset in just the right way and, as a consequence, to eliminate risk. This type of hedging is the basis of more complicated hedging strategies such as those engaged in by investment banks and hedge funds. If investors want to trade their financial derivatives, this equation gives a view of the pricing strategy.
Now let's solve a practical problem using our numerical skills!
Consider an European call style option is made for a security currently trading at 40 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?
The time interval $[0,T]$ is divided into $N$ equally sized subintervals of length $\Delta t$. The price of the underlying asset will take values in the unbounded interval $[0,\infty)$. However, an artifical limit $S_{max}$ is introduced. The interval $[0, S_{max}]$ is divided into $M$ equally sized subintervals of length $\Delta S$.
In [10]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
T = 0.25 # expiration time
N = 2000
dt = T/N
S_max = 40 # price of the underlying asset
M = 200
s = np.linspace(0, S_max, M+1)
r = 0.1 # no-risk interest rate
sigma = 0.4 # volatility of underlying asset
E = 10. # the exercise price
In the following code, we are calculating an European call option. "European" means options can only be exercised at time $t=T$. So we have the final condition such as $C(S,T) = C_T(S)$. However, in a numerical code, we prefer a initial condition rather than final condition. So we apply the change of variable: $t_{new} = T -t$. Thus, the new equation we are solving should look like this (note the sign changes):
$$\begin{equation} \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 \end{equation}$$For the explicit method, we use a forward difference scheme in time $t$, and a central difference scheme in value of underlying asset $S$. From the discretized form of equation, we can rearrange the equation as following:
$$\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}$$We can find that the price of option at $(n+1)$ time step $C_i^{n+1}$ does not depend on $\Delta S$. And since there is $i$ in the equation, we can no longer use vectorized stencil operation instead of double loops.
The boundary conditions are:
$$\begin{equation} C(0,t) = 0\\ C(S_{max},t) = S_{max} - E \end{equation}$$The initial conditions are: $$\begin{equation} C(S,0) = max (S-E,0) \end{equation}$$
In [11]:
# initial condition
C = np.zeros(M+1)
C = s - E
for i in range(M+1):
if C[i] < 0:
C[i] = 0
Cn = np.zeros(M+1)
In [12]:
for n in range(N):
Cn = C.copy()
for i in range(1,M):
C[i] = 0.5 * (sigma**2 * i**2 * dt - r*i*dt) * Cn[i-1] \
+ (1 - sigma**2* i**2 *dt - r*dt) * Cn[i] \
+ 0.5 * (sigma**2 * i**2 * dt + r*i*dt) * Cn[i+1]
print 'the price of the call option should be around {}, \
if the price of stock is 20 dollar.'.format(C[M/2])
Here we get the solution of the European call option problem. Then we compared the numerical result with the result given by a statistical method:
$$\begin{equation} C(S,t) = SN(d_1) -Ee^{-rT}N(d_2) \end{equation}$$N(x) is the cumulative distribution function for a standardised normal random variable (we use norm.cdf to calculate 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} d_1 = \frac{ln(S/E)+(r + \frac{1}{2}\sigma^2)T}{\sigma \sqrt{T}}\\ d_2 = d_1 - \sigma \sqrt{T} \end{equation}$$
In [13]:
from scipy.stats import norm
import math
C_stat = np.zeros(M+1)
for i in range(M+1):
d1 = (math.log1p(s[i]/E) + (r+0.5*sigma**2)*T) / (sigma * math.sqrt(T))
d2 = d1 - (sigma * math.sqrt(T))
C_stat[i] = s[i] * norm.cdf(d1) - E*math.exp(-r*T) * norm.cdf(d2)
In [14]:
plt.figure(figsize=(8,5), dpi=100)
plt.plot(s,C,color='#003366', ls='--', lw=3, label='Computational');
plt.plot(s,C_stat, label='Statistical');
plt.xlabel('S')
plt.ylabel('C')
plt.legend();
We compare the results between numerical method and statistical method. They match well after $S>10$. When $S<10$, the statisical model is not accurate (below zero is meaningless).
In [15]:
from IPython.core.display import HTML
css_file = '../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[15]:
In [15]: