Steven R. Dunbar
Department of Mathematics
203 Avery Hall
Lincoln, NE 68588-0130
http://www.math.unl.edu
Voice: 402-472-3731
Fax: 402-472-8466

Stochastic Processes and

__________________________________________________________________________

Solution of the Black-Scholes Equation

_______________________________________________________________________

Note: These pages are prepared with MathJax. MathJax is an open source JavaScript display engine for mathematics that works in all browsers. See http://mathjax.org for details on supported browsers, accessibility, copy-and-paste, and other features.

_______________________________________________________________________________________________

### Rating

Mathematically Mature: may contain mathematics beyond calculus with proofs.

_______________________________________________________________________________________________

### Section Starter Question

What is the solution method for the Cauchy-Euler type of ordinary diﬀerential equation:

${x}^{2}\frac{{d}^{2}v}{d{x}^{2}}+ax\frac{dv}{dx}+bv=0\phantom{\rule{0.3em}{0ex}}?$

_______________________________________________________________________________________________

### Key Concepts

1. We solve the Black-Scholes equation for the value of a European call option on a security by judicious changes of variables that reduce the equation to the heat equation. The heat equation has a solution formula. Using the solution formula with the changes of variables gives the solution to the Black-Scholes equation.
2. Solving the Black-Scholes equation is an example of how to choose and execute changes of variables to solve a partial diﬀerential equation.

__________________________________________________________________________

### Vocabulary

1. A diﬀerential equation with auxiliary initial conditions and boundary conditions, that is an initial value problem, is said to be well-posed if the solution exists, is unique, and small changes in the equation parameters, the initial conditions or the boundary conditions produce only small changes in the solutions.

__________________________________________________________________________

### Mathematical Ideas

#### Conditions for Solution of the Black-Scholes Equation

We have to start somewhere, and to avoid the problem of deriving everything back to calculus, we will assert that the initial value problem for the heat equation on the real line is well-posed. That is, consider the solution to the partial diﬀerential equation

$\frac{\partial u}{\partial \tau }=\frac{{\partial }^{2}u}{\partial {x}^{2}}\phantom{\rule{2em}{0ex}}-\infty 0$

with the initial condition

$u\left(x,0\right)={u}_{0}\left(x\right).$

Assume the initial condition and the solution satisfy the following technical requirements:

1. ${u}_{0}\left(x\right)$ has at most a ﬁnite number of discontinuities of the jump kind.
2. $\underset{|x|\to \infty }{lim}{u}_{0}\left(x\right){e}^{-a{x}^{2}}=0$ for any $a>0$.
3. $\underset{|x|\to \infty }{lim}u\left(x,\tau \right){e}^{-a{x}^{2}}=0$ for any $a>0$.

Under these mild assumptions, the solution exists for all time and is unique. Most importantly, the solution is represented as

$u\left(x,\tau \right)=\frac{1}{2\sqrt{\pi \tau }}{\int }_{-\infty }^{\infty }{u}_{0}\left(s\right){e}^{\frac{-{\left(x-s\right)}^{2}}{4\tau }}\phantom{\rule{0.3em}{0ex}}ds.$

Remark. This solution can be derived in several diﬀerent ways, the easiest way is to use Fourier transforms. The derivation of this solution representation is standard in any course or book on partial diﬀerential equations.

Remark. Mathematically, the conditions above are unnecessarily restrictive, and can be considerably weakened. However, they will be more than suﬃcient for all practical situations we encounter in mathematical ﬁnance.

Remark. The use of $\tau$ for the time variable (instead of the more natural $t$) is to avoid a conﬂict of notation in the several changes of variables we will soon have to make.

The Black-Scholes terminal value problem for the value $V\left(S,t\right)$ of a European call option on a security with price $S$ at time $t$ is

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

with $V\left(0,t\right)=0$, $V\left(S,t\right)\sim S$ as $S\to \infty$ and

$V\left(S,T\right)=max\left(S-K,0\right).$

Note that this looks a little like the heat equation on the inﬁnite interval in that it has a ﬁrst derivative of the unknown with respect to time and the second derivative of the unknown with respect to the other (space) variable. On the other hand, notice:

1. Each time the unknown is diﬀerentiated with respect to $S$, it also multiplied by the independent variable $S$, so the equation is not a constant coeﬃcient equation.
2. There is a ﬁrst derivative of $V$ with respect to $S$ in the equation.
3. There is a zero-th order derivative term $V$ in the equation.
4. The sign on the second derivative is the opposite of the heat equation form, so the equation is of backward parabolic form.
5. The data of the problem is given at the ﬁnal time $T$ instead of the initial time $0$, consistent with the backward parabolic form of the equation.
6. There is a boundary condition $V\left(0,t\right)=0$ specifying the value of the solution at one sensible boundary of the problem. The boundary is sensible since security values must only be zero or positive. This boundary condition says that any time the security value is 0, then the call value (with strike price $K$) is also worth 0.
7. There is another boundary condition $V\left(S,t\right)\sim S$, as $S\to \infty$, but although this is ﬁnancially sensible, (it says that for very large security prices, the call value with strike price $K$ is approximately $S$) it is more in the nature of a technical condition, and we will ignore it without consequence.

We eliminate each objection with a suitable change of variables. The plan is to change variables to reduce the Black-Scholes terminal value problem to the heat equation, then to use the known solution of the heat equation to represent the solution, and ﬁnally change variables back. This is a standard solution technique in partial diﬀerential equations. All the transformations are standard and well-motivated.

#### Solution of the Black-Scholes Equation

First we take $t=T-\frac{\tau }{\left(1∕2\right){\sigma }^{2}}$ and $S=K{e}^{x}$, and we set

$V\left(S,t\right)=Kv\left(x,\tau \right).$

Remember, $\sigma$ is the volatility, $r$ is the interest rate on a risk-free bond, and $K$ is the strike price. In the changes of variables above, the choice for $t$ reverses the sense of time, changing the problem from backward parabolic to forward parabolic. The choice for $S$ is a well-known transformation based on experience with the Euler equidimensional equation in diﬀerential equations. In addition, the variables have been carefully scaled so as to make the transformed equation expressed in dimensionless quantities. All of these techniques are standard and are covered in most courses and books on partial diﬀerential equations and applied mathematics.

Some extremely wise advice adapted from Stochastic Calculus and Financial Applications by J. Michael Steele, [1, page 186], is appropriate here.

“There is nothing particularly diﬃcult about changing variables and transforming one equation to another, but there is an element of tedium and complexity that slows us down. There is no universal remedy for this molasses eﬀect, but the calculations do seem to go more quickly if one follows a well-deﬁned plan. If we know that $V\left(S,t\right)$ satisﬁes an equation (like the Black-Scholes equation) we are guaranteed that we can make good use of the equation in the derivation of the equation for a new function $v\left(x,\tau \right)$ deﬁned in terms of the old if we write the old $V$ as a function of the new $v$ and write the new $\tau$ and $x$ as functions of the old $t$ and $S$. This order of things puts everything in the direct line of ﬁre of the chain rule; the partial derivatives ${V}_{t}$, ${V}_{S}$ and ${V}_{SS}$ are easy to compute and at the end, the original equation stands ready for immediate use.”

$\tau =\frac{{\sigma }^{2}}{2}\cdot \left(T-t\right)$

and

$x=log\left(\frac{S}{K}\right).$

The ﬁrst derivatives are

$\frac{\partial V}{\partial t}=K\frac{\partial v}{\partial \tau }\cdot \frac{d\tau }{dt}=K\frac{\partial v}{\partial \tau }\cdot \frac{-{\sigma }^{2}}{2}$

and

$\frac{\partial V}{\partial S}=K\frac{\partial v}{\partial x}\cdot \frac{dx}{dS}=K\frac{\partial v}{\partial x}\cdot \frac{1}{S}.$

The second derivative is

$\begin{array}{llll}\hfill \frac{{\partial }^{2}V}{\partial {S}^{2}}& =\frac{\partial }{\partial S}\left(\frac{\partial V}{\partial S}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =\frac{\partial }{\partial S}\left(K\frac{\partial v}{\partial x}\frac{1}{S}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =K\frac{\partial v}{\partial x}\cdot \frac{-1}{{S}^{2}}+K\frac{\partial }{\partial S}\left(\frac{\partial v}{\partial x}\right)\cdot \frac{1}{S}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =K\frac{\partial v}{\partial x}\cdot \frac{-1}{{S}^{2}}+K\frac{\partial }{\partial x}\left(\frac{\partial v}{\partial x}\right)\cdot \frac{dx}{dS}\cdot \frac{1}{S}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =K\frac{\partial v}{\partial x}\cdot \frac{-1}{{S}^{2}}+K\frac{{\partial }^{2}v}{\partial {x}^{2}}\cdot \frac{1}{{S}^{2}}.\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

The terminal condition is

$V\left(S,T\right)=max\left(S-K,0\right)=max\left(K{e}^{x}-K,0\right)$

but $V\left(S,T\right)=Kv\left(x,0\right)$ so $v\left(x,0\right)=max\left({e}^{x}-1,0\right)$.

Now substitute all of the derivatives into the Black-Scholes equation to obtain:

$K\frac{\partial v}{\partial \tau }\cdot \frac{-{\sigma }^{2}}{2}+\frac{{\sigma }^{2}}{2}{S}^{2}\left(K\frac{\partial v}{\partial x}\cdot \frac{-1}{{S}^{2}}+K\frac{{\partial }^{2}v}{\partial {x}^{2}}\cdot \frac{1}{{S}^{2}}\right)+rS\left(K\frac{\partial v}{\partial x}\cdot \frac{1}{S}\right)-rKv=0.$

Now begin the simpliﬁcation:

1. Divide out the common factor $K$.
2. Transpose the $\tau$-derivative to the other side, and divide through by $\frac{{\sigma }^{2}}{2}$.
3. Rename the remaining constant $\frac{r}{{\sigma }^{2}∕2}$ as $k$ which measures the ratio between the risk-free interest rate and the volatility.
4. Cancel the ${S}^{2}$ terms in the second derivative.
5. Cancel the $S$ terms in the ﬁrst derivative.
6. Gather up like order terms.

What remains is the rescaled, constant coeﬃcient equation:

$\frac{\partial v}{\partial \tau }=\frac{{\partial }^{2}v}{\partial {x}^{2}}+\left(k-1\right)\frac{\partial v}{\partial x}-kv.$

1. Now there is only one dimensionless parameter $k$ measuring the risk-free interest rate as a multiple of the volatility and a rescaled time to expiry $\frac{{\sigma }^{2}}{2}T$, not the original $4$ dimensioned quantities $K$, $T$, ${\sigma }^{2}$ and $r$.
2. The equation is deﬁned on the interval $-\infty , since this $x$-interval deﬁnes $0 through the change of variables $S=K{e}^{x}$.
3. The equation now has constant coeﬃcients.

In principle, we could now solve the equation directly.

Instead, we will simplify further by changing the dependent variable scale yet again, by

$v={e}^{\alpha x+\beta \tau }u\left(x,\tau \right)$

where $\alpha$ and $\beta$ are yet to be determined. Using the product rule:

${v}_{\tau }=\beta {e}^{\alpha x+\beta \tau }u+{e}^{\alpha x+\beta \tau }{u}_{\tau }$

and

${v}_{x}=\alpha {e}^{\alpha x+\beta \tau }u+{e}^{\alpha x+\beta \tau }{u}_{x}$

and

${v}_{xx}={\alpha }^{2}{e}^{\alpha x+\beta \tau }u+2\alpha {e}^{\alpha x+\beta \tau }{u}_{x}+{e}^{\alpha x+\beta \tau }{u}_{xx}.$

Put these into our constant coeﬃcient partial diﬀerential equation, divide the common factor of ${e}^{\alpha x+\beta \tau }$ throughout and obtain:

$\beta u+{u}_{\tau }={\alpha }^{2}u+2\alpha {u}_{x}+{u}_{xx}+\left(k-1\right)\left(\alpha u+{u}_{x}\right)-ku.$

Gather like terms:

${u}_{\tau }={u}_{xx}+\left[2\alpha +\left(k-1\right)\right]{u}_{x}+\left[{\alpha }^{2}+\left(k-1\right)\alpha -k-\beta \right]u.$

Choose $\alpha =-\frac{k-1}{2}$ so that the ${u}_{x}$ coeﬃcient is $0$, and then choose $\beta ={\alpha }^{2}+\left(k-1\right)\alpha -k=-\frac{{\left(k+1\right)}^{2}}{4}$ so the $u$ coeﬃcient is likewise $0$. With this choice, the equation is reduced to

${u}_{\tau }={u}_{xx}.$

We need to transform the initial condition too. This transformation is

$\begin{array}{llll}\hfill u\left(x,0\right)& ={e}^{-\left(-\frac{\left(k-1\right)}{2}\right)x-\left(-\frac{{\left(k+1\right)}^{2}}{4}\right)\cdot 0}v\left(x,0\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & ={e}^{\left(\frac{\left(k-1\right)}{2}\right)x}max\left({e}^{x}-1,0\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =max\left({e}^{\left(\frac{\left(k+1\right)}{2}\right)x}-{e}^{\left(\frac{\left(k-1\right)}{2}\right)x},0\right).\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

For future reference, we notice that this function is strictly positive when the argument $x$ is strictly positive, that is ${u}_{0}\left(x\right)>0$ when $x>0$, otherwise, ${u}_{0}\left(x\right)=0$ for $x\le 0$.

We are in the ﬁnal stage since we are ready to apply the heat-equation solution representation formula:

$u\left(x,\tau \right)=\frac{1}{2\sqrt{\pi \tau }}{\int }_{-\infty }^{\infty }{u}_{0}\left(s\right){e}^{-\frac{{\left(x-s\right)}^{2}}{4\tau }}\phantom{\rule{0.3em}{0ex}}ds.$

However, ﬁrst we want to make a change of variable in the integration, by taking $z=\frac{\left(s-x\right)}{\sqrt{2\tau }}$, (and thereby $\phantom{\rule{0.3em}{0ex}}dz=\left(-\frac{1}{\sqrt{2\tau }}\right)\phantom{\rule{0.3em}{0ex}}dx$) so that the integration becomes:

$u\left(x,\tau \right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }{u}_{0}\left(z\sqrt{2\tau }+x\right){e}^{-\frac{{z}^{2}}{2}}\phantom{\rule{0.3em}{0ex}}dz.$

We may as well only integrate over the domain where ${u}_{0}>0$, that is for $z>-\frac{x}{\sqrt{2\tau }}$. On that domain, ${u}_{0}={e}^{\frac{k+1}{2}\left(x+z\sqrt{2\tau }\right)}-{e}^{\frac{k-1}{2}\left(x+z\sqrt{2\tau }\right)}$ so we are down to:

$\frac{1}{\sqrt{2\pi }}{\int }_{-x∕\sqrt{2\tau }}^{\infty }{e}^{\frac{k+1}{2}\left(x+z\sqrt{2\tau }\right)}{e}^{-\frac{{z}^{2}}{2}}\phantom{\rule{0.3em}{0ex}}dz-\frac{1}{\sqrt{2\pi }}{\int }_{-x∕\sqrt{2\tau }}^{\infty }{e}^{\frac{k-1}{2}\left(x+z\sqrt{2\tau }\right)}{e}^{-\frac{{z}^{2}}{2}}\phantom{\rule{0.3em}{0ex}}dz.$

Call the two integrals ${I}_{1}$ and ${I}_{2}$ respectively.

We will evaluate ${I}_{1}$ (the integral with the $k+1$ term) ﬁrst. This is easy, completing the square in the exponent yields a standard, tabulated integral. The exponent is

$\begin{array}{llll}\hfill & \frac{k+1}{2}\left(x+z\sqrt{2\tau }\right)-\frac{{z}^{2}}{2}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & \phantom{\rule{2em}{0ex}}=\left(\frac{-1}{2}\right)\left({z}^{2}-\sqrt{2\tau }\left(k+1\right)z\right)+\left(\frac{k+1}{2}\right)x\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & \phantom{\rule{2em}{0ex}}=\left(\frac{-1}{2}\right)\left({z}^{2}-\sqrt{2\tau }\left(k+1\right)z+\tau \frac{{\left(k+1\right)}^{2}}{2}\right)+\left(\frac{\left(k+1\right)}{2}\right)x+\tau \frac{{\left(k+1\right)}^{2}}{4}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & \phantom{\rule{2em}{0ex}}=\left(\frac{-1}{2}\right){\left(z-\sqrt{\tau ∕2}\left(k+1\right)\right)}^{2}+\frac{\left(k+1\right)x}{2}+\frac{\tau {\left(k+1\right)}^{2}}{4}.\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

Therefore

$\begin{array}{c}\frac{1}{\sqrt{2\pi }}{\int }_{-x∕\sqrt{2\tau }}^{\infty }{e}^{\frac{k+1}{2}\left(x+z\sqrt{2\tau }\right)}{e}^{-\frac{{z}^{2}}{2}}\phantom{\rule{0.3em}{0ex}}dz\\ =\frac{{e}^{\frac{\left(k+1\right)x}{2}+\tau \frac{{\left(k+1\right)}^{2}}{4}}}{\sqrt{2\pi }}{\int }_{-x∕\sqrt{2\tau }}^{\infty }{e}^{-\frac{1}{2}{\left(z-\sqrt{\tau ∕2}\left(k+1\right)\right)}^{2}}\phantom{\rule{0.3em}{0ex}}dz.\end{array}$

Now, change variables again on the integral, choosing $y=z-\sqrt{\tau ∕2}\left(k+1\right)$ so $\phantom{\rule{0.3em}{0ex}}dy=\phantom{\rule{0.3em}{0ex}}dz$, and all we need to change are the limits of integration:

$\frac{{e}^{\frac{\left(k+1\right)x}{2}+\tau \frac{{\left(k+1\right)}^{2}}{4}}}{\sqrt{2\pi }}{\int }_{-x∕\sqrt{2\tau }-\sqrt{\tau ∕2}\left(k+1\right)}^{\infty }{e}^{\left(-\frac{{y}^{2}}{2}\right)}\phantom{\rule{0.3em}{0ex}}dz.$

The integral can be represented in terms of the cumulative distribution function of a normal random variable, usually denoted $\Phi$. That is,

$\Phi \left(d\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{d}{e}^{-\frac{{y}^{2}}{2}}\phantom{\rule{0.3em}{0ex}}dy$

so

${I}_{1}={e}^{\frac{\left(k+1\right)x}{2}+\tau \frac{{\left(k+1\right)}^{2}}{4}}\Phi \left({d}_{1}\right)$

where ${d}_{1}=\frac{x}{\sqrt{2\tau }}+\sqrt{\frac{\tau }{2}}\left(k+1\right)$. Note the use of the symmetry of the integral! The calculation of ${I}_{2}$ is identical, except that $\left(k+1\right)$ is replaced by $\left(k-1\right)$ throughout.

The solution of the transformed heat equation initial value problem is

$u\left(x,\tau \right)={e}^{\frac{\left(k+1\right)x}{2}+\frac{\tau {\left(k+1\right)}^{2}}{4}}\Phi \left({d}_{1}\right)-{e}^{\frac{\left(k-1\right)x}{2}+\frac{\tau {\left(k-1\right)}^{2}}{4}}\Phi \left({d}_{2}\right)$

where ${d}_{1}=\frac{x}{\sqrt{2\tau }}+\sqrt{\frac{\tau }{2}}\left(k+1\right)$ and ${d}_{2}=\frac{x}{\sqrt{2\tau }}+\sqrt{\frac{\tau }{2}}\left(k-1\right)$.

Now we must systematically unwind each of the changes of variables, starting from $u$. First, $v\left(x,\tau \right)={e}^{-\frac{\left(k-1\right)x}{2}-\frac{{\left(k+1\right)}^{2}\tau }{4}}u\left(x,\tau \right)$. Notice how many of the exponentials neatly combine and cancel! Next put $x=log\left(S∕K\right)$, $\tau =\left(\frac{1}{2}\right){\sigma }^{2}\left(T-t\right)$ and $V\left(S,t\right)=Kv\left(x,\tau \right)$.

The ultimate result is the Black-Scholes formula for the value of a European call option at time $T$ with strike price $K$, if the current time is $t$ and the underlying security price is $S$, the risk-free interest rate is $r$ and the volatility is $\sigma$:

$\begin{array}{llll}\hfill V\left(S,t\right)& =S\Phi \left(\frac{log\left(S∕K\right)+\left(r+{\sigma }^{2}∕2\right)\left(T-t\right)}{\sigma \sqrt{T-t}}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & \phantom{\rule{2em}{0ex}}-K{e}^{-r\left(T-t\right)}\Phi \left(\frac{log\left(S∕K\right)+\left(r-{\sigma }^{2}∕2\right)\left(T-t\right)}{\sigma \sqrt{T-t}}\right).\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

Usually one doesn’t see the solution as this full closed form solution. Most versions of the solution write intermediate steps in small pieces, and then present the solution as an algorithm putting the pieces together to obtain the ﬁnal answer. Speciﬁcally, let

$\begin{array}{llll}\hfill {d}_{1}& =\frac{log\left(S∕K\right)+\left(r+{\sigma }^{2}∕2\right)\left(T-t\right)}{\sigma \sqrt{T-t}}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {d}_{2}& =\frac{log\left(S∕K\right)+\left(r-{\sigma }^{2}∕2\right)\left(T-t\right)}{\sigma \sqrt{T-t}}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

and writing ${V}_{C}\left(S,t\right)$ to remind ourselves this is the value of a call option,

${V}_{C}\left(S,t\right)=S\cdot \Phi \left({d}_{1}\right)-K{e}^{-r\left(T-t\right)}\cdot \Phi \left({d}_{2}\right).$

#### Solution of the Black-Scholes Equation Graphically

Consider for purposes of graphical illustration the value of a call option with strike price $K=100$. The risk-free interest rate per year, continuously compounded is 12%, so $r=0.12$, the time to expiration is $T=1$ measured in years, and the standard deviation per year on the return of the stock, or the volatility is $\sigma =0.10$. The value of the call option at maturity plotted over a range of stock prices $70\le S\le 130$ surrounding the strike price is illustrated in 1

We use the Black-Scholes formula above to compute the value of the option before expiration. With the same parameters as above the value of the call option is plotted over a range of stock prices $70\le S\le 130$ at time remaining to expiration $t=1$ (red), $t=0.8$, (orange), $t=0.6$ (yellow), $t=0.4$ (green), $t=0.2$ (blue) and at expiration $t=0$ (black).

Using this graph, notice two trends in the option value:

1. For a ﬁxed time, as the stock price increases the option value increases.
2. As the time to expiration decreases, for a ﬁxed stock value the value of the option decreases to the value at expiration.

We predicted both trends from our intuitive analysis of options, see Options.. The Black-Scholes option pricing formula makes the intuition precise.

We can also plot the solution of the Black-Scholes equation as a function of security price and the time to expiration as value surface.

This value surface shows both trends.

#### Sources

This discussion is drawn from Section 4.2, pages 59–63; Section 4.3, pages 66–69; Section 5.3, pages 75–76; and Section 5.4, pages 77–81 of The Mathematics of Financial Derivatives: A Student Introduction by P. Wilmott, S. Howison, J. Dewynne, Cambridge University Press, Cambridge, 1995. Some ideas are also taken from Chapter 11 of Stochastic Calculus and Financial Applications by J. Michael Steele, Springer, New York, 2001.

_______________________________________________________________________________________________

### Algorithms, Scripts, Simulations

#### Algorithm

For given parameter values, the Black-Scholes-Merton solution formula is sampled at a speciﬁed $m×1$ array of times and at a speciﬁed $1×n$ array of security prices using vectorization and broadcasting. The result can be plotted as functions of the security price as done in the text. The calculation is vectorized for an array of $S$ values and an array of $t$ values, but it is not vectorized for arrays in the parameters $K$, $r$, $T$, and $\sigma$. This approach is taken to illustrate the use of vectorization and broadcasting for eﬃcient evaluation of an array of solution values from a complicated formula.

In particular, the calculation of ${d}_{1}$ and ${d}_{2}$ then uses broadcasting, also called binary singleton expansion, recycling, single-instruction multiple data, threading or replication.

The calculation relies on using the rules for calcuation and handling of inﬁnity and NaN (Not a Number) which come from divisions by $0$, taking logarithms of $0$, and negative numbers and calculating the normal cdf at inﬁnity and negative inﬁnity. The plotting routines will not plot a NaN which accounts for the gap at $S=0$ in the graph line for $t=1$.

#### Scripts

Geogebra

R
1m <- 6;
2n <- 61;
3S0 <- 70;
4S1 <- 130;
5K <- 100;
6r <- 0.12;
7T <- 1.0;
8sigma <- 0.10;
9
10time <- seq(T,0, length=m);
11S <- seq(S0,S1, length=n);
12
13numerd1 <- outer( ( (r + sigma^2/2)*(T-time)), log(S/K), "+");
14numerd2 <- outer( ( (r - sigma^2/2)*(T-time)), log(S/K), "+");
15
16d1 <- numerd1/(sigma*sqrt(T-time));
17d2 <- numerd2/(sigma*sqrt(T-time));
18part1 <- t(t(pnorm(d1))*S);
19part2 <- K*exp(-r*(T-time))*pnorm(d2);
20VC <- part1 - part2;
21
22matplot(S, t(VC), type = "l");
23
Octave
1m = 6;
2n = 61;
3S0 = 70;
4S1 = 130;
5K = 100;
6r = 0.12;
7T = 1.0;
8sigma = 0.10;
9
10time = transpose(linspace(T, 0, m) );
11S = linspace(S0,S1,n);
12
13d1 = (log(S/K) + ( (r + sigma^2/2)*(T-time)))./(sigma*sqrt(T-time));
14d2 = (log(S/K) + ( (r - sigma^2/2)*(T-time)))./(sigma*sqrt(T-time));
15part1 = bsxfun(@times, normcdf(d1), S);
16part2 = bsxfun( @times, K*exp(-r*(T-time)), normcdf(d2));
17VC = part1 - part2;
18
19plot(S,VC)
20
Perl
1sub pnorm {
2    my ( $x,$sigma, $mu ) = @_; 3$sigma = 1 unless defined($sigma); 4$mu    = 0 unless defined($mu); 5 6 return 0.5 * ( 1 + erf( ($x - $mu ) / ( sqrt(2) *$sigma ) ) );
7}
8$m = 6; 9$n = 61;
10$S0 = 70; 11$S1 = 130;
12$K = 100; 13$r = 0.12;
14$T = 1.0; 15$sigma = 0.10;
16
17$time = zeroes($m)->xlinvals($T,0.0); 18$S = zeroes($n)->xlinvals($S0,$S1); 19 20$logSoverK = log($S/$K);
21$n12 = (($r + $sigma**2/2)*($T-$time)); 22$n22 = (($r -$sigma**2/2)*($T-$time));
23$numerd1 =$logSoverK + $n12->transpose; 24$numerd2 = $logSoverK +$n22->transpose;
25$d1 =$numerd1/( ($sigma*sqrt($T-$time))->transpose); 26$d2 = $numerd2/( ($sigma*sqrt($T-$time))->transpose);
27
28$part1 =$S * pnorm($d1); 29$part2 = pnorm($d2) * ($K*exp(-$r*($T-$time)))->transpose; 30$VC = $part1 -$part2;
31
32# file output to use with external plotting programming
33# such as gnuplot, R, octave, etc.
34# Start gnuplot, then from gnuplot prompt
35#    plot for [n=2:7] solution.dat using 1:(column(n)) with lines
36
37
38open( F, ">solution.dat" ) || die "cannot write: $! "; 39wcols$S, $VC, *F; 40close(F); 41 SciPy 1import scipy 2 3m = 6 4n = 61 5S0 = 70 6S1 = 130 7K = 100 8r = 0.12 9T = 1.0 10sigma = 0.10 11 12time = scipy.linspace(T, 0.0, m) 13S = scipy.linspace(S0,S1, n) 14 15logSoverK = scipy.log(S/K) 16n12 = ((r + sigma**2/2)*(T-time)) 17n22 = ((r - sigma**2/2)*(T-time)) 18numerd1 = logSoverK[scipy.newaxis,:] + n12[:,scipy.newaxis] 19numerd2 = logSoverK[scipy.newaxis,:] + n22[:,scipy.newaxis] 20d1 = numerd1/(sigma*scipy.sqrt(T-time)[:,scipy.newaxis]) 21d2 = numerd2/(sigma*scipy.sqrt(T-time)[:,scipy.newaxis]) 22 23from scipy.stats import norm 24part1 = S[scipy.newaxis] * norm.cdf(d1) 25part2 = norm.cdf(d2) * K*scipy.exp(-r*(T-time))[:,scipy.newaxis] 26VC = part1 - part2 27 28# optional file output to use with external plotting programming 29# such as gnuplot, R, octave, etc. 30# Start gnuplot, then from gnuplot prompt 31# plot for [n=2:7] solution.dat using 1:(column(n)) with lines 32 33scipy.savetxt(solution.dat, 34 scipy.column_stack((scipy.transpose(S),scipy.transpose(VC))), 35 fmt=(%4.3f)) __________________________________________________________________________ ### Problems to Work for Understanding 1. Explicitly evaluate the integral ${I}_{2}$ in terms of the c.d.f. $\Phi$ and other elementary functions as was done for the integral ${I}_{1}$. 2. What is the price of a European call option on a non-dividend-paying stock when the stock price is$52, the strike price is $50, the risk-free interest rate is 12% per annum (compounded continuously), the volatility is 30% per annum, and the time to maturity is 3 months? 3. What is the price of a European call option on a non-dividend paying stock when the stock price is$30, the exercise price is $29, the risk-free interest rate is 5%, the volatility is 25% per annum, and the time to maturity is 4 months? 4. Show that the Black-Scholes formula for the price of a call option tends to $max\left(S-K,0\right)$ as $t\to T$. 5. For a particular scripting language of your choice, modify the scripts to create a function within that language that will evaluate the Black-Scholes formula for a call option at a time and security value for given parameters. 6. For a particular scripting language of your choice, modify the scripts to create a script within that language that will evaluate the Black-Scholes formula for a call option at a speciﬁed time for given parameters, and return a function of the security prices that can plotted by calling the function over an interval of security prices. 7. Plot the price of a European call option on a non-dividend paying stock over the stock prices$20 to $40, given that the exercise price is$29, the risk-free interest rate is 5%, the volatility is 25% per annum, and the time to maturity is 4 months.
8. For a particular scripting language of your choice, modify the scripts to create a script within that language that will plot the Black-Scholes solution for ${V}_{C}\left(S,t\right)$ as a surface over the two variables $S$ and $t$.

__________________________________________________________________________

### References

[1]   J. Michael Steele. Stochastic Calculus and Financial Applications. Springer-Verlag, 2001. QA 274.2 S 74.

[2]   Paul Wilmott, S. Howison, and J. Dewynne. The Mathematics of Financial Derivatives. Cambridge University Press, 1995.

__________________________________________________________________________

1. Cornell University, Department of Computer Science, Prof. T. Coleman Rhodes and Prof. R. Jarrow. Numerical Solution of Black-Scholes Equation, Submitted by Chun Fan, Nov. 12, 2002.
2. An applet for calculating the option value. based on the Black-Scholes model. Also contains tips on options, business news and literature on options. Submitted by Yogesh Makkar, November 19, 2003.

__________________________________________________________________________

I check all the information on each page for correctness and typographical errors. Nevertheless, some errors may occur and I would be grateful if you would alert me to such errors. I make every reasonable eﬀort to present current and accurate information for public use, however I do not guarantee the accuracy or timeliness of information on this website. Your use of the information from this website is strictly voluntary and at your risk.

I have checked the links to external sites for usefulness. Links to external websites are provided as a convenience. I do not endorse, control, monitor, or guarantee the information contained in any external website. I don’t guarantee that the links are active at all times. Use the links here with the same caution as you would all information on the Internet. This website reﬂects the thoughts, interests and opinions of its author. They do not explicitly represent oﬃcial positions or policies of my employer.

Information on this website is subject to change without notice.