<< "DiscreteMath`NewCombinatorica`";

<< "Finance`ItosLemma`";

Off[General::"spell"]; 
    Off[General::"spell1"]; $Line = 0;

\(TraditionalForm\`\(\(\ \)\(Derivative\ Pricing\ with\ Mathematica\)\)\)

Andrzej  Kozlowski

Abstract.

We present a Mathematica-based introduction to the theory and practice of financial derivative pricing and risk management and a quick survey of some Mathematica resources available for this purpose.

Key words: Stocks, optiosn, arbitrage, Ito's lemma, Black-Scholes equation

1. Brownian Motion and Securities.

We shall be dealing with securities. Some of these, e. g. stocks  (shares) will be risky, others (bonds) not (this is an ideal world).   So we start by defining a Mathematica view of a security.  We think of a security  as a stochastic process characterized by  several parameters rather than a tradable piece of paper. These parameters will be its volatility, drift, current value, and possibly the dividend it pays. Here is a function that defines a new security.

newSecurity[vol_, dr_, val_, div_:  0] := 
      Module[{volatility = vol, drift = dr, sec, value = val, 
          dividend = div}, sec[setVolatility, w_] := volatility = w; 
        sec[getVolatility] := volatility; sec[setDrift, d_] := drift = d; 
        sec[getDrift] := drift; sec[setDividend, d_] := dividend = d; 
        sec[getDividend] := dividend; sec[setPrice, v_] := value = v; 
        sec[getPrice] := value; sec]

We create "stocks" like this:

NEC = newSecurity[v, d, val];

NEC is now actually a function, but you can thing of it as an "object" with built in "methods", getVolatility, setVolatility, getDrift etc.

NEC[setVolatility, 0.29]; NEC[setDrift, 0.2]; NEC[setPrice, 550]; 
    NEC /@ {getVolatility, getDrift, getPrice}

We would like to be able to simulate and plot the time series representing the change in price of our security. Essentially all models  in finance are based on Brownian motion (Wiener process). There are almost countless ways to construct approximations to Brownian motion in Mathematica. The one I will use here is based on approximation by binomial random walks (see below).  I shall call these random walks Brownian paths, although they are of course only good approximations when n is large.

BrownianPath[
          W_ : W] := (Clear[W]; W[n_][0] = 
            0; W[n_][
              i_] := W[n][
                i] = W[n][
                  i - 1] + (((-1))^Random)[Integer]/N[(Sqrt[n])]);

Using the central limit theorem we can show that as n -> ∞ the distribution of W[n][t] does indeed tend to the normal distribution with mean 0 and standard deviation (Sqrt[t]) and that W[n] has all the other properties of the standard Brownian motion (e.g. see [1] or [7]). Note the usage of  "dynamic programming". Its purpose is to create a concrete Brownian path   (random walk) with BrownianPath[], which we can keep using until we reset the definition by evaluating BrownianPath[] again. Of course one can achieve the same thing by using SeedRandom[] as in [10], but our approach is more convenient, particularly when working with several independent random walks.
Now we  define a function  that generates a random walk  of the values of our security, based on the (approximate) geometric Brownian motion (the exponential of a Brownian motion) with volatility and drift taken from our security.

SetAttributes[diffuse, HoldAll]
diffuse[s_, n_][i_] := 
      s[getPrice] (E^(((i s[getDrift])/n) + s[getVolatility] W[n][i]))

Finally, we define a function which plots random walks representing the movement of prices of our security:

SetAttributes[diffusionPlot, HoldAll];
diffusionPlot[s_, n_, color_ : RGBColor[0, 0, 0]] := 
        ListPlot[Table[{(i/n), diffuse[s, n][i]}, {i, n}], 
          PlotJoined -> True, PlotStyle -> color];

BrownianPath[]; 
    plot1 = Block[{$DisplayFunction = Identity}, 
        diffusionPlot[NEC, 1024, RGBColor[1, 0, 0]]];

We can change the volatility of our security and compare the corresponding time series, based on the same BrownianPath.

NEC[setVolatility, 1]; 
    plot2 = Block[{$DisplayFunction = Identity}, 
        diffusionPlot[NEC, 1024, RGBColor[0, 1, 0]]]; 
    Show[plot1, plot2];

Much more about generating, manipulating and plotting diffusions (using a quite different approach)  can be found in [10].

In what follows we shall make the usual assumptions of about the market in which we are operating: completeness and the absence of arbitrage. The latter basically asserts that the best way to make money without risk is by putting it in a bank account and collecting interest. This is sometimes described as the principle that "there is no such thing as a  free lunch". (I have always thought that this ought to be "there is no such thing as a risk-free lunch", but that may sound too unnervingly true in the age of BSE etc.) 
 More precisely, it means that any two portfolios (see below) which always bring the same outcome should have the same value. If this was not true then by selling one and buying the other we could earn our lunch, risk free. In particular, the absence of arbitrage implies the following result. Suppose we have some risk free scheme that guarantees us income V after time T. Then the price of such a scheme ought to be: (e^(-r T)) V, where r is the bank rate of interest. Indeed, this is the amount of money that deposited in a bank at interest rate r will return V after time T.

2.Options and Arbitrage.

At this point we need to explain what we mean by a "portfolio". Our portfolios will be rather limited compared with those in the real world. They will consist of stocks and only one kind of derivatives, options on stocks.  We shall assume we can hold any real amount of these, both positive and negative. This latter fact, the ability to hold negative amount of stock and options may surprise those who have not had much to do with financial markets, but this is indeed possible. (This is related the other condition mentioned above, the completeness of our market). In effect having a negative amount of a share is the same as borrowing it. In reality borrowing shares (known as short selling) does incur a certain small cost but we shall assume that it can be done for free. In the same way as stock we can hold (in positive and negative quantities) options on stock. The options we shall first consider will be of two kinds: the so called European calls and puts (more precisely, what is known as "European vanilla calls and puts"). A European call option with a strike price E and maturity date T gives the holder the right (but not the obligation) to buy the underlying asset at time T for the price E. (In our case the underlying asset will be always one of the securities we defined above, but in the real world one can trade options on a variety of assets, including options on options etc. Our approach can be extended to include these cases).  European put option with a strike price E and maturity date T gives the holder the right to sell the underlying asset at time T for the price E. To hold a negative amount of an option is the same as "writing" the respective option, in other words, to undertake the obligation to purchase (if required)  Each financial instrument stock, put or call has a pay-off function,

SetAttributes[europeanCall, HoldAll];
SetAttributes[europeanPut, HoldAll];
SetAttributes[stock, HoldAll];
payoff[x_ + y_, t_] := payoff[x, t] + payoff[y, t];
payoff[a_?NumericQ x_, t_] := a payoff[x, t];
payoff[europeanCall[T_, e_, St_], t_] := 
        With[{s = St[getPrice]}, 
          If[t == T, (s - e) UnitStep[s - e], 0]];
payoff[europeanPut[T_, e_, St_], t_] := 
        With[{s = St[getPrice]}, 
          If[t == T, (e - s) UnitStep[e - s], 0]];
payoff[stock[St_], t_] := With[{s = St[getPrice]}, s]

Note that we distinguish the "abstract" security (essentially an index) like NEC and  a unit tradable amount of NEC stock, denoted by stock[NEC]. A portfolio is now is an element of the real vector space generated by the basic instruments (calls, puts and stock), that is a "formal sum" of these. 
We shall now use the functions we have developed to prove one of the oldest  theorems about options, the so called "Put-Call" parity.  The proof is a triviality but it illustrates a characteristic use of the "no free lunch" principle, and in some sense is a prototype of everything that follows. Let us first define a perfectly general stock

generalStock = newSecurity[vol, dr, val];

We now construct a "risk free" portfolio with generalStock as the underlying asset. By risk free I mean one that produced exactly the same payoff no matter what happens.

R = -europeanCall[T, e, generalStock] + 
          europeanPut[T, e, generalStock] + stock[generalStock];

Let's compute the payoff of this portfolio at time T. This is trivial to do by hand but let Mathematica do the work for us. There are three cases to consider, when on the expiration date s > e, s < e, and s == e. We check that the payoff of the portfolio is e in each case.

FullSimplify[payoff[R, T] == e, #[val, e]] & /@ {Less, Greater, 
        Equal}

We see that our portfolio gives us a guaranteed payoff of e (==s) at time T, no matter what happens to the stock! Hence by the principle of no arbitrage the values of the portfolio must be (e^(-r T)) e Thus we get the putt-call parity : value of a Eropean put =(e^(-r T)) (strike price)+value of a European call-price of stock.

3.Trees and Derivative Pricing

One of the most remarkable achievements in mathematical finance  was the development by Black, Scholes and Merton of the modern theory of pricing of derivative instruments. The origin of the subject  is the following  simple but profound observation. Let's imagine the following situation. Suppose a stock has current value of 2 dollars and we know that only two things can happen to it: it can either go up by 0.5 dollars or go down by the same amount. Suppose the probability of it going up is 2/3. Suppose we are offered the option of purchasing the stock for 2 dollars at the next instant. How much should such an option be worth?  Let us assume that the rate of interest is 0. The classic probabilistic answer would have been: if the price of the stalk goes up to 2.5 we shall buy it for 2.0 dollars gaining 0.5 dollar but if it goes down we do not exercise the option.Thus our expected gain is ((2 )/3) 0.5 + (0/3)=0.333333. However, this intuitive answer is wrong. The correct answer is obtained as follows. First of all we replace the "actual"  probabilities by probabilities which ensure that the expected value of our stock is exactly the same as the starting value (such a process is called a martingale). In other words we want to find p from:

Solve[1.5 (1 - p) + 2.5 p == 2, p]

The answer is the expected value based on this new "probability", i.e.  (0.5/2) + (0/2) =  0.25`. The reason why this is the "correct price" of the option is that any other price would allow possibilities of arbitrage. For example, suppose the price of the option was 0.2 dollars. We could buy the option for 0.2 dollars and we short sell 0.5 unit of stock for 1 dollar. Now suppose the stock goes up to 2.5. We use the option to buy the stock for 2.5, return it, gaining 0.5 and we loose 0.5(2 -2.5)= 0.25 on the short sale, giving us the net gain of 0.25. 
Or, if the stock falls to 1.5, we do not exercise the option but we gain  0.5(2-1.5) = 0.25 on the short sale. In both cases we have a net, risk free gain of 0.05 dollar. 
Similar argument can be used to show that risk free gain is possible if the option price is higher than 0.25 dollar.
Note that the actual probabilities play no role in determining the value of the option. This idea can ba extended to binary trees of arbitrary length. Here we illustrate it with a three stage process, with actual probabilities and martingale probabilities

prices = {80, 60, 120, 36, 72, 80, 180}; probs = Table[(1/2), {6}]; 
    gr1 = Block[{$DisplayFunction = Identity}, 
        ShowGraph[CompleteBinaryTree[8 - 1], VertexStyle -> Disk[0.04], 
          VertexLabel -> prices, EdgeLabel -> probs, 
          VertexColor -> Yellow, PlotRange -> All]]; 
    martingaleProbs = {(2/3), (1/3), (2/3), (1/3), (2/5), (3/5)}; 
    gr2 = Block[{$DisplayFunction = Identity}, 
        ShowGraph[CompleteBinaryTree[8 - 1], VertexStyle -> Disk[0.04], 
          VertexLabel -> prices, EdgeLabel -> martingaleProbs, 
          VertexColor -> Yellow, PlotRange -> All]]; 
    Show[GraphicsArray[{gr1, gr2}]];

This elementary example contains the essence of one approach to both the theory and practical computation of option values. We can use binary trees and discrete martingales to compute approximate value of any option. The virtue of this approach is that it is very easy to implement in Mathematica. For the American put option (see below) this was implemented in Mathematica by Ross.  Miller [6] and for several other well known options by S. Benninga, R. Steinmetz and J. Stroughair [2]. However, this approach is rather inefficient in general the differential equations alternatives offer more  accuracy and performance.

4. The Black-Scholes Differential Equation

There is another approach to option pricing which does not make the use of martingale theory. Let us consider a general stock (without dividends). In the Black-Scholes analysis it is assumed that such a stock performs an infinitesimal random walk given by the stochastic differential equation

\(TraditionalForm\`d\((X)\) == 
      X\ \((μ\ d\((t)\) + σ\ d\((W)\))\)\)

where σ and μ are constants called the volatility and the drift of the stock. Here W stands for the standard Brownian motion.  Of course S and W are random variables defined on some sample space which one never needs to bother about. (Essentially the space of all possible "continuous random walks").  It can be easily shown that this equation describes a limiting case of a discrete random walk of the type we defined above. Precise meaning can be given to it by considering it in terms of stochastic integrals. To get an intuitive idea of what is going on consider the fundamental theorem of classical calculus:

f(t)=f(0)+(∫_0^t)f'(s)ds

One wy to think of this is as a statement that every differentiable function can be approximated by straight line segments. Of course a Brownian motion is exactly something that can't be approximated in this way. It is therefore natural to consider a wider class of objects, stochastic processes that can be written in the form:

\(TraditionalForm\`X\_t = ∫\_0\%t\( μ\_s\) ds + ∫\_0\%t\( σ\_s\) dW\_s + 
        X\_0\)

Here (Subscript[W, s]) stands for Brownian motion,  μ and σ are random processes which are pre-visible with respect to the Brownian motion (meaning that (Subscript[σ, s])and (Subscript[μ, s])depend only on (Subscript[W, t]) for t ≤ sand the ordinary integral \(TraditionalForm\`∫\_0\%t\((\(\(σ\_s\^2\)\(+\)\) | μ\_s | )\) dt\)is finite for all t (with probability 1). Such a process (which can be defined precisely) is called an Ito process. It is often written in differential form:

d(Subscript[X, t])=(Subscript[μ, s]) ds+(Subscript[σ, s]) d(Subscript[W, s])

For example, the Ito  process used in the Black-Scholes model is given by the stochastic differential equation

d(Subscript[X, t])=(Subscript[X, t ])μ ds+(Subscript[X, t ])σ d(Subscript[W, s])

where μ and σ are constants. Let us now suppose that we have a derivative product which depends only on the value of the stock at time t. In other words, the value of the derivative is a function Val((Subscript[X, t]), t).The famous formula due to K. Ito ("Ito's Lemma") tells us that a smooth transformation of an Ito process is again an Ito process and tells us how to obtain a representation of the form (3). This leads to a whole subject called "Ito calculus", which is a generalization of traditional calculus  to Ito processes. A particularly powerful Mathematica implementation of the Ito calculus by Wilfrid Kendall can be found in [4] and another  one in [10] in the same volume.  Here we shall use of a simpler  package inspired by the latter, due to Mark Fisher and available from MathSource.

{TimeSymbol, TimeIncrement, BrownianIncrement, CorrelationSymbol} = {t, 
        dt, dW, ρ}; 
    dX = ItoMake[X[t], μ X[t], σ X[t]];

Now suppose at time t0 we hold a portfolio consisting of our derivative V and -Δ X of stock, where

Δ = (∂_X) Val[X, 
              t] /. t -> t0;

In other words, our portfolio P is

P = Val[X[t], t] - Δ X[t];

Let's compute the change of this portfolio at time t0,

ItoD[P] /. t -> t0

Note that the instantaneous change in the value of the portfolio is deterministic (no Brownian component). By the "no free-lunch" principle it must therefore be equal to what could be earned by an amount equal to P invested at the risk free interest rate for a period dt.

(r P dt /. {X[t] -> X, 
            t -> t0}) == %

ExpandAll[Thread[(#(1/dt) &)[%], 
            Equal] /. {t0 -> 
              t} /. a_ == b_ -> 
          b - a == 0]

This is the famous Black-Scholes differential equation.It is a second order linear partial differential equation of the parabolic type.Note that it is satisfied by every derivative security whose price value depends only on the current value of X and on t.A slight modification of the argument shows that if the stock pays a constant dividend q the differential equation takes the form

(1/2) (X^2) Derivative[2, 0][Val](X, t) (σ^2)-r Val(X, t)+Derivative[0, 1][Val](X, t)+(r - q) X Derivative[1, 0][Val](X, t)==0

This equation must be satisfied by the value of every derivative whose price depends only on the current value of X and on t. The differences between different kinds of options will express itself in the boundary conditions. 
At this point there are a number there are two basic methods to proceed if w try to solve the differential equation. One approach using Mathematica is described in the article by Steel and Stine. It is based on the so called Feynman-Kac formula which expresses the solutions of certain second order PDE's in terms of expectation of the corresponding diffusion (solution of a SDE). The other approach is preferred by those who do not know or do not wish to know much about stochastic techniques. By means of  certain standard transformations (for details see the book of Shaw [9] or if you prefer more explanation but no Mathematica code the book of Wilmott at al) one can reduce the Black-Scholes equation to one of the most famous and best studied PDE's in physics: the heat equation,  which can be solved using one of several well known techniques.  For the European call the explicit, Nobel prize winning solution is:

Norm[z_?NumberQ] := 0.5 Erf[(z/(Sqrt[2]))] + 0.5; 
    Norm[x_] := (1/2) (Erf[(x/(Sqrt[2]))] + 1); 
    (Subscript[d, 2])[s_, v_, k_, t_, r_, 
        q_] := (((r - q) t + Log[(s/k)])/(v (Sqrt[t]))) - ((v (Sqrt[t]))/2); 
    (Subscript[d, 1])[s_, v_, k_, t_, r_, 
        q_] := (((Sqrt[t]) v)/2) + (((r - q) t + Log[(s/k)])/(v (Sqrt[t])));
BlackScholesCall[s_, k_, v_, r_, q_, t_] := 
        s (E^(-q t)) Norm[(Subscript[d, 1])[s, v, k, t, r, q]] - 
          k (E^(-r t)) Norm[
              (Subscript[d, 2])[s, v, k, t, r, q]];

5.Risk Management

The purpose of financial derivatives is management of risk. This is done by constructing portfolio's stocks and options with a desired payoff. We know that we cannot guarantee gain above the free interest rate, but  if we have any hypothesis or opinion about the future of a stock (or a stock index) we can build a portfolio that will do well if this hypothesis turns out to be true. We can built risky portfolio's with a potential for high gain, or a safe portfolios, minimizing possible loss.  In fact it is possible to show [8] that any "contingent claim" can be matched by a portfolio of European options. 
Now that we can compute the value of a European call or a put we can easily create a Mathematica function that will compute the value of a portfolio, consisting of different calls, puts and stocks. This makes it easy to create a suitable portfolio for any specific purpose. For this purpose we define the function value, which computes the value of a portfolio, which we think of as a formal sum of different derivatives.

value[x_ + y_, r_] := value[x, r] + value[y, r];
value[a_?NumericQ x_, r_] := a value[x, r];
value[europeanCall[T_, e_, s_], r_] = 
        BlackScholesCall[s[getPrice], e, s[getVolatility], r, s[getDividend], 
          T];
value[europeanPut[T_, e_, s_], 
        r_] := (E^(-r T)) e + 
        value[europeanCall[T, e, s], r] - value[stock[s], r]
value[stock[s_], r_] := s[getPrice]

Let us construct a portfolio based on the NEC stock we have defined, which will make profits provided the stock price stays within a certain range.

NEC[setPrice, 550]; NEC[setVolatility, 0.1]; 
    myPortfolio = 
      europeanCall[1, 535, NEC] - 2 europeanCall[1, 550, NEC] + 
        europeanCall[1, 565, NEC]; price = value[myPortfolio, 0.05]; 
    NEC[setPrice, val]; 
    Plot[Evaluate[{price, payoff[myPortfolio, 1]}], {val, 500, 600}, 
      PlotStyle -> {RGBColor[0, 1, 0], {Thickness[0.01], 
            RGBColor[1, 0, 0]}}];

The part of the picture where the thick red graph (payoff) is above the thin green (cost) represents profitable range of final values of the stock, the rest represents a loss. This kind of portfolio is known as a straddle and is very popular with traders.

6. Other Types of Derivatives and Options

So far we have considered  only European options, but the great majority of options actually traded on stock exchanges are so called American options. They differ form European ones in that they can be exercised any time before the expiry date. In the case of call options without dividends this difference turns out to be illusory, for one can show that it is never optimal to exercise an such an option before expiry. However, in the case of American put options when the price of the stock falls low it may be optimal to exercise the option early. The result is that the Black-Scholes equation becomes only an inequality. The problem of finding the value of the American put option also reduces to a well known problem in physics  (the "obstacle problem") and it is known that no analytic solution to it exists. One has therefore to result to other methods. In [2] and [6] this is done by means of binary trees and the Monte carlo method, but much more attractive are the methods considered in [9], which include  approximate analytic solutions and the use of finite difference methods of solving differential equations. The same applies to other options, known as "exotic" Asian options,  look-back options etc.

7. Second Generation Option Pricing

We have only discussed here the use of Mathematica in the Black-Scholes theory , which is actually the "first generation"  of option pricing models. Among it's obvious limitations is the assumption of constant volatility of stocks, which no one who has paid the slightest attention to the behavior of the stock markets in recent years could believe to be realistic. New generation of theories that do not rely on this assumption are the objects of current research. Of course Mathematica is just as useful here as can be seen from the very interesting book [5].

Two freely available Mathematica packages were used in this notebook. They are: the new version of the Combinatorica package by Sriram Pemmaraju and Steven Skiena <http://www.cs.uiowa.edu/~sriram/Combinatorica/> and the "Ito's Lemma" package by Mark Fisher available from MathSource.