This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0340. It is primarily for students who have some experience using Mathematica. If you have never used Mathematica before and would like to learn more of the basics for this computer algebra system, it is strongly recommended looking at the APMA 0330 tutorial. As a friendly reminder, don't forget to clear variables in use and/or the kernel. The Mathematica commands in this tutorial are all written in bold black font, while Mathematica output is in normal font.

Finally, you can copy and paste all commands into your Mathematica notebook, change the parameters, and run them because the tutorial is under the terms of the GNU General Public License (GPL). You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. The tutorial accompanies the textbook Applied Differential Equations. The Primary Course by Vladimir Dobrushkin, CRC Press, 2015;

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica

Black Scholes model

Fischer Black Myron Scholes Robert C. Merton

An options contract is a financial instrument which gives the owner the right, but not the requirement, to buy or sell stock on a particular date (the expiration date) for a specified price (the strike price). A call contract is the contract that provides the owner the right to purchase stock at the strike price on the expiration date, and a put contract allows the owner to sell stock at the strike price on the expiration date. These contracts are commonly bought and sold on public exchanges, and options markets exist for most publicly traded stocks.

Given that the contract has value to the owner and places a liability on the seller, it stands to reason that each contract has a fair value at any given point in time. In other words, if we assume that a buyer and a seller of an options contract have the same imperfect information about changes in the stock price, there should be some price at which the buyer does not stand to have a positive expected gain greater than the return provided by the risk free interest rate.

The Black Scholes model, also known as the Black--Scholes--Merton model, is a model of price variation over time of financial instruments such as stocks that can, among other things, be used to determine the price of a European call option. The model assumes the price of heavily traded assets follows a geometric Brownian motion with constant drift and volatility. When applied to a stock option, the model incorporates the constant price variation of the stock, the time value of money, the option's strike price, and the time to the option's expiry.

The Black Scholes model is one of the most important concepts in modern financial theory. It was developed in 1973 by Fisher Black, Robert Merton and Myron Scholes and is still widely used now. It is regarded as one of the best ways of determining fair prices of options. The Black Scholes model requires five input variables: the strike price of an option, the current stock price, the time to expiration, the risk-free rate, and the volatility. Additionally, the model assumes stock prices follow a lognormal distribution because asset prices cannot be negative. Moreover, the model assumes there are no transaction costs or taxes; the risk-free interest rate is constant for all maturities; short selling of securities with use of proceeds is permitted; and there are no riskless arbitrage opportunities.

The Merton model is an analysis model – named after economist Robert C. Merton – used to assess the credit risk of a company's debt. Analysts and investors utilize the Merton model to understand how capable a company is at meeting financial obligations, servicing its debt, and weighing the general possibility that it will go into credit default. This model was later built out by Fischer Black and Myron Scholes to develop the Black--Scholes pricing model.

Robert C. Merton is a famed American economist and Nobel Memorial Prize laureate, who befittingly purchased his first stock at age 10. Later, he earned a Bachelor in Science at Columbia University, a Masters of Science at California Institute of Technology (Cal Tech), and a doctorate in economics at Massachusetts Institute of Technology (MIT), where he later become a professor until 1988. At MIT, he developed and published groundbreaking and precedent-setting ideas to be utilized in the financial world.

Prior to the development of the Black-Scholes options pricing model, numerical techniques were often used to estimate the fair price of an options contract. As we will see later, an options contract is considered fairly priced if there is no way to use some combination of buying/selling stocks and options to earn a risk free profit greater than the return on a risk free asset.

A common technique used to model the stock market is to assume that the price of a stock follows Geometric Brownian Motion with positive drift, meaning that the percent-wise movement of the stock price is assumed to be random and normally distributed with positive mean. It is also assumed that there exists a risk free interest rate, meaning that there is some percent-wise return that can be earned on invested money with no risk. Generally this is interpreted to mean the interest rate on U.S. Treasury bonds, or some other instrument which has a very low risk of default. Lastly, the stock is assumed to pay no dividend. The Black--Scholes model makes certain assumptions:

In the real world, extremely large movements in stock prices are much more common than these assumptions would suggest. Financial crises and large macroeconomic disruptions can often cause severe downward movements in stock markets that would be nearly impossible if the market followed a true Geometric Brownian Motion.

Despite this difference, and a few others, the above assumptions are still used as simplifying assumptions to develop the pricing model for options contracts.

We denote the stock price as a function of time S(t) and the price of the call contract as a function of both the stock price and time V(S,t). Given our assumption that the stock price follows a Geometric Brownian Motion, we can write that

\[ {\text d}S = \mu S\,{\text d}t + \sigma S\,{\text d}X , \]
where X follows Brownian Motion. Intuitively, this indicates that changes in the stock price are related to some positive drift term that is proportional to the stock price, and to some random movement, whose average magnitude is also proportional to the stock price.

Suppose we are to buy the call option and sell the stock short, meaning that we borrow shares of stock from someone else and sell them at the current market price, giving us the obligation to purchase the shares at a later date and return them. Our position thus has a net value of

\[ \Pi = V(S,t) - \Delta S(t) , \]
where delta is the amount of the stock that we short. (In the real world, it is not possible to short fractions of a stock, but these values can be scaled up to very close integer ratios, so our analysis is still viable.)

From this equation, we can get

\[ {\text d}\Pi = {\text d}V(S,t) - \Delta {\text d}S(t) . \]

The next step requires us to use Ito’s Lemma, a lemma that is used to calculate the differential of a function of time and a stochastic function, which is exactly what he stated the option price to be. See Evans [5] for more information on Ito’s Lemma. Using Ito’s Lemma, we can arrive at

\[ {\text d}V = \frac{\partial V}{\partial S}\,{\text d}S + \frac{\partial V}{\partial t}\,{\text d}t + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\,{\text d}t . \]
Plugging into our equation for the value of our portfolio, we get:
\[ {\text d}\Pi = \frac{\partial V}{\partial S}\,{\text d}S + \frac{\partial V}{\partial t}\,{\text d}t + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\,{\text d}t - \Delta {\text d}S. \]
It would be theoretically possible to constantly adjust delta so that
\[ \Delta = \frac{\partial V}{\partial S} . \]
We therefore cancel the first and last term to arrive at
\[ {\text d}\Pi = \frac{\partial V}{\partial S}\,{\text d}S + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\,{\text d}t . \]
Because it should not be possible to earn a risk free profit off of the option that is greater than the risk free rate of return (which we will denote here as r), we also have
\[ {\text d}\Pi = r\,\Pi \,{\text d}t = r \left( V - \Delta S \right) {\text d}t = r \left( V - \frac{\partial V}{\partial S}\,S \right) {\text d}t \]
because the risk free rate will earn the product of the value of our portfolio, the rate of return, and time elapsed. Therefore,
\[ \frac{\partial V}{\partial t}\, {\text d}t + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}\,{\text d}t = r \left( V - \frac{\partial V}{\partial S}\,S \right) {\text d}t . \]
Upon dropping dt, it follows
\[ \frac{\partial V}{\partial t} + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} = r \left( V - \frac{\partial V}{\partial S}\,S \right) . \]
This gives the Black--Scholes equation:
\[ \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} - r\, V =0 . \]

The price of an option V(S, t) is defined for 0 < S < ∞ and 0 &lel tT because a stock price is between 0 and infinity and there is a fixed time T until expiration. The boundary conditions are as follows:

\[ \begin{split} V(S,T) &= \max \left( S - K, 0 \right) , \\ V(0,t) &= 0 , \end{split} \]
where K is the strike price of the call. The first condition states that at expiration, the contract is worth the difference between the stock price and the strike price. This is because the contract gives the owner the right to buy stock at price K at time T, and so if S > K, then the contract is worth S - K. However, if S < K, then the contract is worthless because there is no point in exercising your right to buy at K when you can just buy at the market price S, which is less. The second condition states that the call is worthless if the stock price is 0. Because we assumed that the stock would follow Geometric Brownian Motion, a price of 0 would imply that the price would stay at 0 forever, since movements are proportional in magnitude to the stock price. Therefore there is no value to a call contract if the price of the stock is 0 (assuming the strike price is positive).

We apply the transformation \( u = V\,e^{ −rt} \) (or \( V = u\,e^{ rt} \) ) and accordingly \( e^{rt} \frac{\partial u}{\partial S} = \frac{\partial V}{\partial S} \) and \( e^{rt} \frac{\partial^2 u}{\partial S^2} = \frac{\partial^2 V}{\partial S^2} \) to the Black--Scholes equation

\[ \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} - r\, V =0 . \]
to obtain
\[ r\,e^{rt} u + e^{rt} \frac{\partial u}{\partial t} + \frac{1}{2}\,\sigma^2 S^2 e^{rt} \frac{\partial^2 V}{\partial S^2} + rS\,e^{rt} \frac{\partial u}{\partial S} - r\,e^{rt} u =0 . \]
This is simplified to
\[ e^{rt} \frac{\partial u}{\partial t} + \frac{1}{2}\,\sigma^2 S^2 e^{rt} \frac{\partial^2 V}{\partial S^2} + rS\,e^{rt} \frac{\partial u}{\partial S} =0 , \]
which upon eliminating the common exponential multiple yields
\[ \frac{\partial u}{\partial t} + \frac{1}{2}\,\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\, \frac{\partial u}{\partial S} =0 , \]
From here, we can attempt to make the substitutions \( S = e^y \) and \( t= T-\tau \) to transform u(S,t) to u(y,τ). This should likewise transform \( \frac{\partial}{\partial S} \) to \( \frac{1}{S}\,\frac{\partial}{\partial y} \) and \( \frac{\partial}{\partial t} \) to \( -\frac{\partial}{\partial \tau} . \) We can also derive
\[ \frac{\partial}{\partial S^2} = - \frac{1}{S^2}\,\frac{\partial}{\partial y} + \frac{1}{S} \,\frac{\partial}{\partial S}\,\frac{\partial}{\partial y} = - \frac{1}{S^2}\,\frac{\partial}{\partial y} + \frac{1}{S^2} \,\frac{\partial}{\partial y^2} . \]
Making these substitutions and grouping the \( \frac{\partial u}{\partial y} \) terms, we arrive at
\[ - \frac{\partial u}{\partial \tau} + \frac{1}{2}\,\sigma^2 \frac{\partial^2 u}{\partial y^2} + \left( r - \frac{1}{2}\,\sigma^2 \right) \frac{\partial u}{\partial y} =0 . \]
We use one last transformation,
\[ z = y + \left( r - \frac{1}{2}\,\sigma^2 \right) \overline{\tau} \qquad\mbox{and} \qquad \overline{\tau} = \tau , \]
and the corresponding relationships \( \frac{\partial}{\partial y} = \frac{\partial}{\partial z} \) and
\[ \frac{\partial}{\partial \overline{\tau}} = \frac{\partial}{\partial \tau} \, \frac{\partial \tau}{\partial \overline{\tau}} + \frac{\partial}{\partial y} \, \frac{\partial y}{\partial \overline{\tau}} = \frac{\partial}{\partial \tau} - \left( r - \frac{1}{2}\,\sigma^2 \right) \frac{\partial}{\partial y} . \]
We can substitute all these equations into Black--Scholes equation to obtain
\[ - \frac{\partial u}{\partial \overline{\tau}} - \left( r - \frac{1}{2}\, \sigma^2 \right) \frac{\partial u}{\partial z} + \frac{1}{2}\,\sigma^2 \frac{\partial^2 u}{\partial z^2} + \left( r - \frac{1}{2}\, \sigma^2 \right) \frac{\partial u}{\partial z} =0 , \]
which can be simplified further to the heat equation:
\[ \frac{\partial u}{\partial \overline{\tau}} = \frac{1}{2}\,\sigma^2 \frac{\partial^2 u}{\partial z^2} \]
The initial condition is transformed to
\[ u(z,\tau ) = e^{-r(T- \tau )} \, V \left( e^{z - (r - \sigma^2 /2 )} , T- \tau \right) \qquad \Longrightarrow \qquad u(z,0) = e^{-rT} \,V \left( e^z , T \right) = e^{-rT} \,\max \left( 0, e^z - K \right) . \]
Our final result is simply the heat diffusion equation \( \frac{\partial u}{\partial t} = k\,\frac{\partial^2 u}{\partial x^2} \) subject to the initial condition u(x,0) = φ(x). Its solution is known to be
\[ u(x,t) = \frac{1}{\sqrt{4\pi kt}} \, \int_{-\infty}^{\infty} e^{-(x-y)^2 /(4kt)} \, \varphi (y) \,{\text d}y . \]
If we plug our values into this solution, and do some algebra to replace some of the transformed variables with the originals (using different variable names to avoid confusion with the ones used before), we have
\[ V(S,t) = \frac{e^{-r\tau}}{\sqrt{2\pi \sigma^2 \tau}} \,\int_{-\infty}^{\infty} \exp \left\{ - \frac{\left( r- \sigma^2 /2 \right) \tau + \ln (S) - \eta}{2\sigma^2 \tau} \right\} \max \left( 0, e^{\eta - K} \right) {\text d} \eta , \]
which evaluates to the Black--Scholes Formula
\[ V(S,\tau ) = S\,\Phi \left( d_1 (S, \tau ) \right) - K\,e^{-r\tau} \Phi \left( d_2 (S, \tau ) \right) , \]
\[ \begin{split} d_1 (S, \tau ) &= \frac{1}{\sigma\sqrt{\tau}} \left[ \ln \frac{S}{K} + \left( r + \frac{1}{2}\, \sigma^2 \right) \tau \right] , \\ d_2 (S, \tau ) &= \frac{1}{\sigma\sqrt{\tau}} \left[ \ln \frac{S}{K} + \left( r - \frac{1}{2}\, \sigma^2 \right) \tau \right] , \end{split} \]
where Φ is the normal cumulative distribution function and τ = T - t.

Using the solution to the Black-Scholes equation, we can simulate the price of a call or put contract expiring on Sept 1, 2018 with various strike prices, starting a year before that on Sept 1, 2017. MATLAB code used for the approximation is shown below:

K = 2500:100:3100; %strike prices
n = length(K);

for i = 1:n
cblackscholes(K(i), false)
hold on
xlabel(’Time Elapsed (years)’)
ylabel(’Call Price ($)’)
title(’Black-Scholes Simulated Call Prices for S&P500, 8/17-8/18’)
hold off

% put
for i = 1:n
cblackscholes(K(i), true)
hold on
xlabel(’Time Elapsed (years)’)
ylabel(’Put Price ($)’)
title(’Black-Scholes Simulated Put Prices for S&P500, 8/17-8/18’)
hold off

function cblackscholes(K,isput)
% constants
r = 0.02; % risk free interest rate
sigma = xlsread(’VIX’,’sheet’); % volatility
T = 1; %day
S = xlsread(’SP500’,’sheet’);
n = length(S);
t = linspace(0,1,n);
call = linspace(0,0,n);
for i = 1:n
d1 = (1/(sigma(i)*sqrt(T-t(i))))*(log(S(i)/K)+(r+(sigma(i)^2/2))*(T-t(i)));
d2 = d1 - sigma(i)*sqrt(T-t(i));
call(i) = N(d1)*S(i) - N(d2)*K*exp(-r.*(T-t(i))); % black scholes forumla
if isput
call(i) = call(i) - S(i) + (K * exp(-r * (T - t(i))));
function wt = N(d) % cdf of normal distribution
pd = makedist(’Normal’);
wt = cdf(pd,d);
Using this code and data VIX and SP500, we generate the following graphs: ../../image
Black Scholes simulated Put Prices S&P500 8/17--8/18 Black Scholes simulated Call Prices S&P500 8/17--8/18
  1. Black, Fischer and Scholes, Myron (1973). "The Pricing of Options and Corporate Liabilities". Journal of Political Economy. 81 (3): 637--654. doi: 10.1086/260062
  2. Bodie, Zvi, Kane, Alex, and Marcus, Alan, Investments, 11th Edition, McGraw-Hill Education, 2018.
  3. Dunbar, Steven R., Solution of the Black-Scholes Equation, Stochastic Processes and Advanced Mathematical Finance.
  4. Esekon, Joseph Eyang’an, Analytic solution of a nonlinear Black--Scholes equation, International Journal of Pure and Applied Mathematics, Volume 82 No. 4 2013, 547--555.
  5. Evans, Lawrence, An Introduction to Stochastic Differential Equations, Version 1.2, Department of Mathematics of UC Berkeley. Berkeley, CA.
  6. Hull, John C. (1997). Options, Futures, and Other Derivatives. Prentice Hall. ISBN 0-13-601589-1.
  7. Kumar, S, Yildirim, A, Khan, Y, Jafari, H, Sayevand, K, and Wei, L, Analytical solution of fractional Black--Scholes European option pricing equation by using Laplace transform, ournal of Fractional Calculus and Applications, Vol. 2. Jan 2012, No. 8, pp. 1-9. ISSN: 2090-5858.
  8. Manafian, Jalil and Paknezhad, Mahnaz, Analytical solutions for the Black--Scholes equation, Applications and Applied Mathematics: An International Journal, Vol. 12, Issue 2, (December 2017), pp. 843--852.
  9. Merton, Robert C. (1973). "Theory of Rational Option Pricing". Bell Journal of Economics and Management Science. The RAND Corporation. 4 (1): 141–183. doi: 10.2307/3003143.
  10. Stecher, Michael, Converting the Black-Scholes PDE to the heat equation
  11. Stehlíková, Beáta, Black-Scholes model: Derivation and solution
  12. Yalincak, Orhun Hakan. Criticism of the Black-Scholes Model: But Why is It Still Used?: (The Answer is Simpler than the Formula). New York University. New York, NY, 2005.
  13. Zhu, S.P., A closed-form analytical solution for the valuation of convertible bonds withconstant dividend yield, ANZIAM J.47:477–494 (2006)
  14. Zhu, S.P., An exact and explicit solution for the valuation of American put options,Quant. Finan.6:229–242 (2006).


Return to Mathematica page

Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions