The Problem

The Solution

1900 	Not quite correct  solution given by Bachelier.

1973	Black and Scholes and Merton :   The Black-Scholes formula.

The Black-Scholes Model

First considered by P. Samuelson (Nobel prize winner in Economics)  in 1970 ,  
 Describes the movement of prices on a stock exchange by the stochastic differential equation:

((Subscript[dS, t])/(Subscript[S, t])) = μ  dt  +  σ (Subscript[dW, t])

(Subscript[S, t]) is  stock price at time t
(Subscript[dS, t])  change in price of stock in time dt
 (Subscript[dW, t])  a random increment in "standard Brownian motion"
μ and σ are the "drift" and the "volatility" of the stock.

Simulation of the Black-Scholes Model

<< "Statistics`NormalDistribution`";
SimulateBlackScholesModel(S_, μ_, σ_, T_, 
            n_ : 100):=Module[{dt = (T/n)},Show[ListPlot[{Range[0, T, dt], NestList[#1 Random[NormalDistribution(dt μ + 1, σ (Sqrt[dt]))] &, S, n]}^T, 
                    PlotJoined -> True, 
                    PlotStyle -> RGBColor[1, 0, 0], 
                    PlotRange -> All, AxesOrigin -> {0, 0}, 
                    Background -> RGBColor[0, 1, 1], 
                    DisplayFunction -> Identity],Graphics[Line[(({{0, 0}, {0, S}}))]]]];

Show[
        Table[SimulateBlackScholesModel(100, 0.2, 0.5, 1, 300), {20}], 
        DisplayFunction -> $DisplayFunction];

The Black-Scholes-Merton Formula

Current  value of a European Call  option is:

\(TraditionalForm\`V(S, 0) = 
      1\/2\ e\^\(\(-r\)\ T\)\ S\ \((erf(\(\(\@T\ σ\)\/2 + \(r\ T + log(S\/K)\)\/\(\@T\ σ\)\)\/\@2) + 1)\) - 
        1\/2\ e\^\(\(-r\)\ T\)\ K\ \((erf(\(\(r\ T + log(S\/K)\)\/\(\@T\ σ\) - \(\@T\ σ\)\/2\)\/\@2) + 
              1)\)\)

K -  strike price,
S  - current  price of  stock 
T -  strike time
r - risk free rate of interest
erf  the error function:   erf(z)=(2/(Sqrt[π]))(∫_0^z)(e^(-(t^2)))dt

The Black-Scholes PDE

Black and Scholes  reduced  problem to a boundary value problem.  
The Black-Scholes Partial Differential Equation has the form:

(1/2) (S^2) (((∂^2) V(S, t))/(∂(S^2))) (σ^2)-r V(S, t)+((∂V(S, t))/(∂t))+r S ((∂V(S, t))/(∂S))=0

V(S, t)  value of the option at time t
   S price of stock at time t. 
   Boundary conditions:

V(0, t) =  
      0  for all   0 ≤  t  ≤  
        T
V(S, T) = Max[S - K, 0]

The Proof.

The original proof of Black-Scholes was intuitively convincing but logically incorrect.
 Since then many other proofs have been given, using methods of probability theory. 
 Our main contribution here consists of two things:

1.  We give a logically correct version of the original Black and Scholes argument.
2. We perform virtually the entire proof by means of the computer algebra program Mathematica. 

Our work also demonstrates the rapidly increasing abilities of  mathematical computer software.

Implementation of Ito Calculus

The key tool:Stochastic Calculus (Ito Calculus) based on the notion of Ito Integral.

SetAttributes[d, Listable];
ItoRules = {((d(W(t)))^2) -> 
            d(t), d(W(t)) d(t) -> 0, ((d(t))^2) -> 
            0};
ItoProcesses = {};
ItoProcessQ(W(t)) = True;
ItoMake(
          Y_[t], μ_, σ_) := (ItoProcessQ(Y(t)) = True; 
          Last[Last[
              ItoProcesses = 
                ItoProcesses ⋃ {d(
                        Y(t)) -> μ d(t) + σ d(
                            W(t))}]]);
ItoGather(expr_) := 
        Coefficient[expr, d(t)] d(t) + 
          Coefficient[expr, d(W(t))] d(W(t));
ItoD(expr_) := 
        Module[{s, 
            itos = Cases[expr, _?(ItoProcessQ(#1) &), ∞]}, 
          ItoGather((ExpandAll[Normal[((O(s))^3) + (expr /. Flatten[{Thread[itos -> itos + s d(itos)], 
                                    t -> 
                                      t + s d(
                                        t)}])] /. s -> 1 /. ItoProcesses] /. ItoRules) - expr)];

Here we create an Ito process:

ItoMake(Y(t), μ(t), σ(t))

The the famous Ito Lemma:

ItoD(f(U(t), Y(t)))

% /. σ(t) -> 0

Usual chain rule.

Portfolios and their Value

ItoD(V(t)) = X(t). ItoD(p (t))

The Absence of Arbitrage Principle

The most important principle in modern finance is the principle of absence of arbitrage.  
Comparable to the principle of conservation of energy in physics.

The Idea of the derivation of the Black-Scholes PDE.

We want to find the value of a derivative security whose payoff depends only on the value of stock at some future time T. 
The idea  is to construct a portfolio which is instantaneously risk free and self financing.  Then:

Absence of Arbitrage  
 ⇓
  Instantaneous  return on such a portfolio  = Instantaneous return on the same amount invested in a bank account.

This will give the Black-Scholes PDE.

The Derivation

We start by creating two Ito processes, one representing our stock S[t] and the other a cash bond B[t] giving a risk-free return r.

ItoProcesses = {};
ItoMake(
        S(t), μ S(t), σ S(t));
ItoMake(B(t), r B(t), 0);

We can check the Ito processes we have created:

ItoProcesses

The following quantity called the Delta plays a major role in option pricing:

Δ(t_):=((∂V(S(t), t))/(∂S(t)));

Derivation continued

We define the value of a portfolio for a given price vector:

val(port_List, price_List) := port . price

We define the following portfolio process consisting of 1 option, -Δ(t) units of stock, and u(S(t), t) units of a cash bond, where u is some unknown function that we wan to determine.

portfolio(t_) = {1, -Δ(t), 
          u(S(t), t)};

We define the corresponding price process:

price(t_) := {V(S(t), t), S(t), B(t)}

Combine the definitions to get the portfolio value process:

p(t_) = val(price(t), portfolio(t))

Derivation continued

First we obtain the condition for our portfolio to be instantaneously risk free:

RF=First[Solve[Coefficient[ItoD(p(t)), d(W(t))] == 0,Derivative[1, 0][u](S(t), t)]]

For the portfolio to be self- financing the increase in value must be equal to that obtained only from the change in the  prices. In other word, SF below has to be zero (c.f. equation (5))

SF = 
      Simplify[ItoD(p(t)) - portfolio(t) . ItoD(price(t))]

We can now choose u to be a function of x and t that will make the portfolio risk free:

ss = RF /. S(t) -> x

dd = 
      First[ss] /. Rule -> Subtract

ss = 
      First[DSolve[dd == 0, u(x, t), x, t]]

u(x_, t_) = u(x, t) /. ss

This is the general form of u. Note that it consists of an amount of bond that will exactly cancel out the option and stock held, and an amount of cash (Subscript[c, 1])[t], which is a smooth function of t. In fact, as we shall see, (Subscript[c, 1])[t] depends on the path the stock has followed up to time t:, but has no "brownian" component. Observe that the total value of our portfolio depends only on (Subscript[c, 1]):

Simplify[portfolio(t)]

Simplify[p(t)]

Derivation completed

Now let's check what becomes of the self-financing condition with this value of u:

SF = Simplify[SF]

To make this zero we have to choose our function (Subscript[c, 1])[t] so that it's derivative with respect to time satisfies the above equation:

ss=First[Solve[SF == 0,(Subscript[c, 1])'(t)]]

Let us now assume that the bond grow at constant interest rate r with the initial value 1:

B(t_) := (E^(r t))

Since the portfolio is instantaneously risk-free and self-financing we have:

Simplify[
        ItoD(p(t)) == r p(t) d(t) /. ss, 
        d(t) > 0] /. a_ == b_ :> 
        Simplify[((b - a)/2)] == 0

We have derived the Black-Scholes PDE.

Solving the PDE.

Remove["Global`*"]

There are several ways to solve the Black-Scholes equation. At least one approach has been used to derive the famous formula of Black and Scholes for the value of vanilla Euopean options  with the help of Mathematica [Steele and Stine]. This is the method based on the use of the famous Feynman-Kac theorem and thus requires  a certain  knowledge of concpets of probablity and stochastic analysis. Our aim here  is to accomplish the same purpose but by another well known route which does not require any such knowledge; we use  Mathematica to reduce the Black-Scholes equation to the Heat (Diffusion)  equation using the standard methods that can be found in most texts on options pricing, e.g Willmott and Dawynne. We then use the standard "non-dimensionality" method to reduce the heat equation to an ordinary differential equation, which we solve using the mathematica function DSolve.  (We remark here that at this time (version 5) Mathematica's DSolve can  solve neither the general Black-Scholes equation nor the Heat equation, though NDSolve (as of version 5) can solve them numerically). The main point of what folows is to demonstrate hwo we can use Mathematica to preform the rather tedious algebraic transfromations needed to obtain the solutions. The methods are quite general and can of course be used in solving  other equations of this kind.

Before we start we shall need to add a few definitions to Mathematica that will simplify the manipulations that follow. We shall find it more convenient to use the total derivative Dt  we need to make certain modifications to it.

Mathematica Preliminaries

Unprotect[Dt];
Dt[True, _, ___] := True
Protect[Dt];
MakeConstants[l_List] := 
        Scan[SetAttributes[#, Constant] &, l];
Independent[x_, y_] := (x /: Dt[x, ___, y, ___] = 0; 
          x /: Dt[x, ___, {y, _}, ___] = 0; y /: Dt[y, ___, x, ___] = 0; 
          y /: Dt[y, ___, {x, _}, ___] = 0;);
Off[Solve::"ifun"];
Off[Eliminate::"ifun"];

Transforming Black-Scholes into Heat Equation.

We first make sure that the following pairs of variables will be treated as independent

\(TraditionalForm\`\(\(
\)\(\(Independent[S, t];\)

    \(Independent[x, τ];\)

    \(Independent[S, τ];\)

    \(Independent[x, t];\)\)\)\)

Next comes a list of constants in the following computation:

MakeConstants[{K, σ, T, α, β, 
        r, σ}]

Now we define the Black-Scholes equation:

BS=(1/2) (S^2) (((\[DifferentialD]^2) V(S, 
                      t))/(\[DifferentialD](S^2))) (σ^2)+r S ((\[DifferentialD]V(S, 
                      t))/(\[DifferentialD]S))+((\[DifferentialD]V(S, 
                    t))/(\[DifferentialD]t))-r V(S, t)==0;

We also have the following intital condition:

initials = {V[S, T] == 
            Max[S - K, 0]};

The idea is to use the following substitutions to reduce the equation to one with constant coefficients.

subs = {S == K (E^x), 
          t == T - (τ/(((σ^2)/2))), 
          V(S, t) == K v(x, τ)};

However, we need more that the above since we shall also need to substitute for the derivatives. Since Mathematica will not automatically deduce any needed identities involving derivatives from the above, we differentiate the above "manually". This gives us the complete set of identities which we will then apply to obtian equations in terms of the new variables.

allsubs = 
        Flatten[Through[{# &, Dt[#, {S, 1}] &, Dt[#, {S, 2}] &, 
                    Dt[#, {t, 1}] &}[#]] & /@ subs];

eq=Refine[Eliminate[Flatten[{BS, allsubs}],{V(S, t),((\[DifferentialD]V(S, 
                          t))/(\[DifferentialD]t)),((\[DifferentialD]V(S, 
                          t))/(\[DifferentialD]S)),(((\[DifferentialD]^2) V(S, 
                          t))/(\[DifferentialD](S^2)))},InverseFunctions -> False],allsubs]

eq = 
      FullSimplify[
        eq /. r -> k *(1/2) (σ^2), σ ≠ 
          0]

The next step is to substitute for the function v[x, t] the expression \(TraditionalForm\`e\^\(x\ α + β\ τ\)\ \(u(x, τ)\)\) where α and β are to be determined later.

v[x_, τ_] := (E^(x α + β τ)) u[
          x, τ]

eq = 
      FullSimplify[
        eq, σ ≠ 0 && Element[_, Reals]]

ss=First[Solve[{Coefficient[First[eq], u[x, τ]] == 0,Coefficient[First[eq],Derivative[1, 0][u][x, τ]]==0},{α, β}]]

eq = FullSimplify[eq /. ss]

We also would like to convert the intial conditon to one expressed in the new variables. While it is easily done by hand we can also do it using the above method.

FullSimplify[u[x, 0]/.First[Solve[{V(S, T)==(S - K) θ(S - K),V(S, T) ==  K v[x, 0]},u[x, 0],{V[S, T]}]/.S -> (E^x) K/.ss],{K > 0}]

A slight glitch in Matheamtica 5.0 caused to not to remove K.

(Subscript[u, 0])[
        x_] := u[x, 
          0] = ((E^((1/2) (k + 1) x)) - (E^((1/2) (k - 1) x)))*
          UnitStep[(E^((1/2) (k + 1) x)) - (E^((1/2) (k - 1) x))]

Solving the Heat Equation

From the general theory of th e Green function the solution is now given by:

(∫_(-∞)^∞) 
                (Subscript[u, 0])(s) (E^(-((((x - s))^2)/(4 τ)))) \[DifferentialD]s

FullSimplify[Integrate[(Subscript[u, 0])(
                  s) (E^(-((((x - s))^2)/(4 τ)))),{s, -∞, ∞},Assumptions->{τ > 0,x > 0,k∈ℝ}],{τ > 0,x > 0,k∈ℝ}]

FullSimplify[(Sqrt[π]) (Sqrt[τ]) ((E^((1/4) (k + 1) (2 x + k τ + τ))) (erf(((x + k τ + τ)/(2 (Sqrt[τ])))) + 
                    1) + (E^((1/4) (k - 1) (2 x + (k - 1) τ))) (erfc(((x + (k - 1) τ)/(2 (Sqrt[τ])))) - 2)) /. {k -> 2 (r/(σ^2)), τ -> 
            (1/2) (T (σ^2) - t (σ^2)), 
          x -> log((S/K))}, {σ > 0, T >= t, 
        t ≥ 0}]