APPLICATIONS OF MATHEMAITCA TO TBE STOCHASTIC CALCULUS

J.1dichaelS~ 1 & Robert A. Stine'

Robert A. Stine, Dept. of Stadsdes~lhe Wharton School, University of Pa, PA 19104~

Key Words: BlackScholes, diffusion, Assume that the functions g a: (R xR+) >R

s~lic computing. are members of L2 and let Wt denote the stan

Introduction dard Weiner process. A fundamental result

of Itt's theory of stochastic integrals is that

Symbolic computing has established its there exists a welldefined process Xt that can

value in dealing with many of the more te be written suggestively as

dious aspects of differentiation and integra t t

tion. Our goal here is to show that symbolic ;Q =X0 + s) dt + s) dWs . (1)

computing is also a powerful tool for the study illXSI JCIIXRI

of the stochastic calculus. The vital intuition behind Xt and this repre

As an illustration of the power of this ap sentation is that the functionp determines the

proach, we consider the classic option pricing drift of Xt whereas a controls its variability.

problem of mathematical finance and use For example, if Xt represents the price of a

symbolic computations to obtain the Black stock, then p and a could be interpreted as the

Scholes formula. Our objective is to derive associated rate'of return and risk.

this formula for pricing an option rather than Generally (1) is expressed in differential

manipulate the wellknown final result. form as

After a brief review of diffusions, we de dXt = XXt, 0 d t + a (Xt, 0 d Wt. (2)

scribe a symbolic representation for a diffu Note that the process Wt is not of bounded

sion and illustrate how this representation is variation so that the integral in (1) can not be

sufficient for simulating and plotting real interpreted naively. The interested reader

izations of diffusions. We then consider can consult, for example, the book of Arnold

transformations of diffusions via the Itt (1974) or Duffle (1988) for further motivation

formula. The Itt formula is the basis of our and details.

symbolic algorithm that determines the dif

fusion defined by a transformation of an

other. g's Formula

While this paper uses Mathematica, these

same ideas could also be implemented in

other symbolic environments, such as d's formula is the key ingredient that under

Macsyma. Our choice of Mathematica was lies many of our symbolic manipulations of dif

based on its graphics, programmability, and fusions. Assume for the function f mapping

large user base. If one has access to this sys (RR+)* R that the derivatives

tem, the commands shown here can be repli fx=df(x,t)ldx, fxx =d2f(x,t)ldx2,

cated after installing the definitions of the and

functions. ft =df(X, t) / d t

A catalog of the details about the use of exist. If Yt = f(Xtt), then It~'s Formula shows

Mathematica appear in Wolfram (1991); the that Yt is also a diffusion,

more gentle tutorial of Gray and Glynn

(1991) provides a good introduction. Other ap dYt

plications of symbolic computing in statistics 2

appear in Steele (1985) and Stine (1990). (ft (Xt, t) + gXt, t) fx(Xt, t) + fxx a (Xt, t) 12) dt

+ t fx(Xt,, 0 o(Xt) ) dWt (3)

This expression is not, however, in the

A diffusion Xt is identified by two func canonical form (2) since the drift and disper

tions g and a and its initial value X0. sion are functions of Xt rather than 1 Yt. If we

assume that f has an inverse g so that tion at time t .5 T where St denotes the price of

Y = fix) <> X = g(y), the stock at time t.

we obtain the desired form by substituting The key feature of the BlackScholes formu

g(Yt) for Xt in (3), obtaining lation of this problem is to introduce two diffu

sions. Assume that the market consists of two

dYt = (ft(g(Yt),t)+p(g(Yt),t)fx(Xtt) + investments, stocks and bonds. The prices of

f= J&Yt),012) dt shares of stock and bond St and Bt are treated

+ ( fx(k(Yt), 0 o(Xt, 0 ) d Wt as diffusions of known form,

dSt = p St dt + aSt dWt

= 9 Xt,0 dt + ~& Xt,0 dWt . (4) dBt= r Bt dt .

The bond process represents a riskftee asset

A key task for Mathematica is to manage the with guaranteed rate of return r, whereas the

mapping ( p(xt), o(xt) (i~yt), 5(Ylt) stock diffusion is a more risky asset offering

higher potential returns with some chance of

less. Here we ignore dividends and transac

The Infinitesimal Generator tion costs and assume that the market. does not

permit arbitrage. This brief introduction is

The infinitesimal generator simplifies ftom Duffle (1988, §22) which contains farther

the evaluation of the It~ formula. Given a motivation and details.

diffusion Xt, the infinitesimal generator A The BlackScholes approach duplicates the

measures the infinitesimal rate of change in value of the option V(St,0 with a portfolio con

the expected value of f(Xt, t) given that Xt sisting of at shares of stock and bt shares of

starts from x, bond. We need to determine at and bt. The

EX f (Xt 1 t) fix 1 0) portfolio representation for the value of the op

Af(x) = lim t tion implies that V(St,0 is itself a diffusion,

as t~0, where EX denotes the expectation having the form

conditional on XO = x. Viewed as an operator, V(St,0 = at St + bt Bt .

It follows that

the infinitesimal associated with Xt is the dV(St, 0 = at A + bt dBt (6)

sum = (at g St + bt r Bt) dt + at aSt dWt

A=.~+P(Xlt)i+1120(~,t) so that the drift of V(Stt) isatMSt+btrBt

a 49x ar2 and the dispersion is at aSt. On the other hand,

This expression is useful since the drift the It~ formula also gives expressions for the

function in (3) is the result of applying the in mean and dispersion of V(Stt),

finitesimal associated with Xt to the function dV(St , 0 = A1V1 dt + Vx(St, 0 a St d Wt , (7)

f. Substituting (5) into (3), we obtain a much where Vx denotes the derivative of V with re

more manageable expression for Ith"'s for spect to its first argument. Equating the drift

mula (3), and dispersion functions of (6) and (7) yields a

dYt = df(Xt) system of equations which determine the coef

= Affi dt + fx(Xt, 0 otXt, 0 dWt ficients at and bt of the matching portfolio,

at = Vx(St,0

. bx VXX(St, t) CV2St2 + Vt(St, t)

BlackScholes Q2tion aicing Mpdel rBt

where V= denotes the second derivative of V

All of these concepts are captured in the fol with respect to its first argument, and Vt is the

lowing simple variation of the options pricing first derivative with respect to its second argu

problem. The problem is to determine the cur ment. Substituting these solutions into the port

rent market value of an option to purchase folio representation (6) gives a partial differ

shares of a stock at some future time T for a set ential equation for V,

price P. Such options are known as European W(XW + Vt(xt) + rXMXW

options. Let V(St, 0 denote the value of the op + cp2x*Vxx(xt)12 = 0.

12

This PDE can be solved by various methods, printing diffusions, are imported ftom an

including the FeyninanKac theorem. external file. If one has access to

If A is the infinitesimal of diffusion Xt, the Mathematica, these expressions can be en

FeynnianKac theorem states that the solu tered directly into the program. Note that

tion of the partial differential equation terminating a statement with a semicolon

suppresses printing the results of that state

AIV(x,t)lrV(x,t)=0 (0i5t.5T) (8) ment.

We begin the session by importing the ex

with boundary conditionV(xT) g(xT) is ternal definitions from the file Ito.m, and

P(Tt) then create two diffusion objects.

WX,0 = e E[g(XTT) 1 Xt xl. (9)

rn[l]:

For the BlackScholes problem, the payoff <

ary condition. For example, the payoff func In (21:

tion for a European option is the segmented bin MakeDiffusion[woorl,Ol;

function exp MakeDiffusion[X,r X,s X,v]

g(x) = (X P)' Out[21:

where P is the exercise price of the option. 1 diffusion[X,r X,s X,v]

Although (9) looks rather impressive, We

later show how it simplifies in the case to a The items bm and exp are Mathematica ob

lognormal integration which is solved by jects representing Brownian motion

standard methods. dWt = Odt + 1 dWt

and the diffusion

A Data Structure fQr Diffusions Xt = rXt dt + sXt dWt; XO = v. (10)

Notice that the name exp identifies the diffu

In order to use Mathernatica when work sion object, whereas the argument X in the

ing on problems like deriving the Black second input is a symbol identifying the dif

rm'ula, we first need the ability to fusion in the expressions for the drift and

Scholes fo ' dispersion functions p and a.

create and manipulate diffusions. Our While it is tempting to use the symbol X as

choice here is simple, yet contains all of the the name of the diffusion object at line 2, we

needed information. Basically, the data cannot do so here. Such statements lead to in

structure is a list with four elements: the finite recursion.

symbol which identifies the diffusion, the

functions p and a, and finally the initial InM:

value. The functions p and a are entered as X=MakeDiffusion[X,r X,s X,v]

expressions of the identifying symbol and t.

Rather than just put these items together in a General::recursion:

list, however, we make the head of the expres Recursion depth of 256 exceeded.

sion the identifier diffusion. This subtle

change permits some abstraction and pro OutM:

vides the opportunity for typechecking in $Interrupted

functions that operate on diffusions. In order

to insure proper initialization, each diffu The system repeatedly substitutes the entire

sion is created by a call to the creator function expression for each X appearing in the right

MakeDiffusion. hand side of W31.

To demonstrate the use of Mathematica, Automatic evaluation of arguments of

the following dialog shows boxed snapshots functions can cause other problems. For ex

from a Mathematica session. Mathematica ample, if x had previously been assigned the

expressions are shown in Courier type and value 2, as by the expression x = 2, then

are numbered in italics. The dialog contains Mathematica would interpret the second line

definitions for many of the functions that we

use. Othersl/including those for making and

of M21 as

exp MakeDiffusion[2,r 2,s 2,v),

and confusion would rein. To insure against this problem, the definition of the function MakeDif fusion begins

MakeDiffusion[s_Symbol, mu,._, disp, init_].

The leading argument s_Syrnbol provides type checking: the function is only called when the argument a is a symbol, and not, for example, a number.

To enhance this objectoriented approach we use a set of "methods" for manipulating diffusions. Although Mathematica does not provide a full set of objectoriented programming facilities, we can extend the abstraction by providing accessor methods. For example, the functions

symbol[d diffusion] :=d[[111;

drift[ddiffusion] :dH211;

dispersion[d,_diffusion) :d[[311;

initialValue[ddiffusion]:=d[[4]1;

extract the components of a diffusion by in

dexing into the data structure. As long as all

other routines access the drift function p via

the function drif t (rather than indexing di

rectly), we can * change the internal represen

tation of a diffusion by changing these acces

sor functions rather than finding every place

that we extract the drift component.

To illustrate the use of the diffusion data structure, consider how to display a diffusion. While the expression

diffusion[X,r X,s X,v] in Out[21 is indicative of the mathematical notation in (10), it leaves a lot to be desired. Mathematica offers, however, some limited formatting capabilities that permit us to define the function P rintDi f f us ion which yields a more familiar rendering,

In[C:

PrintDiffusion[exp]

Out [41 :

dX (r X ) dt + (a X ) dW ; X = v

t t t t0

The data structure also contains sufficient information for producing simulated realization&. A discrete approximation to Xt suggests that we can "walY the diffusion out from its starting point by incrementing time by small amounts and applying the functions p and a recursively. To do this easily in Mathematica, break the problem into two parts. First, build a function that generates the next value in the discrete approximation given the current value:

In[SI:

nextValue[(t,current_)] : Block xl Mu[current, tl dt; x2 = sig[current, tInRand[dt]; {t+dt,current+xl+x2) // N 1;

In[C:

1 dataNestList[nextvalue,(0,xO),100];

Here nRand [ x 1 produces a zeromean Gaussian random variable with variance x. Like MakeDi f f us i on and P rintDi f f us i on, nRand is defined in the external file Ito.m. The simulated realization is generated by applying nextvalue recursively via the builtin function NestList. NestList returns the recursively defined list

(nextvalue[(0,xO)],

nextValue[nextValue[(0,xO}]],

The function Realize manages the associated variable and function initializations:

In[71:

Realize[d diffusion, n, dt_l

Block[(x,mu,sigma,symsymbol[d],x01,

mu[x t_l (drift[d] /.SYM>x);

,sig[2~_,t_](dispersion[dll.sym>x);

xO N[initialValue[d]3;

NestList[nextValue,10,xO},nl 1;

Realize creates the functions mu and sig used in nextvalue, with the symbol of the diffusion replaced by the generic identifier x. Had we instead attempted todefine mu with the line mu [sym t 1 drift[d]; Mathematica wo7uld look for the symbol sym in drif t, rather than look for the symbol which is the value of sym

7

As an example, the following steps generate three realizations of the diffusion

1 dNt = Nt dt + Nt dWt CV0 10). The mapping function (identified by avoids looping, and Show plots the three realizations together on the same axes.

Xn[BI:

dif MakeDiffusion[n,n,n,10];

SimDif ' Realize[#,200,.005]& M Table[dif,(3)1;

plta=ListPlot[#,PlotJoined>Truel& 1@ simDif;

Show[pltsl

Out [81 : =

35.

30

25

20

is

110

5

1k

0:2 0.,4 0,6 0,8 '1

Finally, we need to describe how to add, subtract, and multiply diffusions. These instructions are given as rules associated with the diffusion type. For example, the following rule describes how to add two diffusions:

diffusion /:

diffusion[sl_,fl_,gl_,il_l +

diffusion[s2_,f2_,g2_,i2_1 =

diffusion[SUM, fl+f2, gi+g2,

il+i21;

The leading 'W' associates this rule for addition with diffusion objects rather than the protected internal function "+". It is implicit in this rule that both diffusions are defined in terms of the same Weiner process. When the drift and dispersion functions do not depend on the diffusion symbols, as with Brownian motion, this simple rule works fine. Notice, however, that the symbolic name for the new diffusion, sum, does not appear in the resulting drift and dispersion functions. Thus, this representation does not

conform to the representation (2) sincep and a are not functions of the leading symbol SuK In many cases it is not possible to obtain this form. For example, consider

InM:

xDiffMakeD!ffusion[X,Sqrt[X],2,0]; yDiffMakeDiffusion[Y,Sqrt[Y1,3,01; zDiff xDiff + yDiff;

Out [91: m l diffusi'on[SUM,Sqrt[XI+Sqrt[Y],,5,,01

In this case, it is not possible to express the

drift X1 12+yl 12 as a function ofZ = X+Y

alone. Fortunately, in many cases such as

the BlackScholes illustration below, we can

leave M and a in the form produced by this

simple rule. Our expanded system notes that

it cannot make this substitution and marks

the diffusio 1 n object in such a way that other

routines, such as that for finding the in

finitesimal, are aware of the problem.

A

Implementating the Ito Form la

We begin with the infinitesimal generator. In keeping with the form of (5), we can implement the infinitesimal as a function which consists of derivative operators. The second two derivatives of (5) are captured by the function derivop, which is the key element of evaluating the infinitesimal function. The #l& pairs in derivop denote abbreviated function definitions, similar to the lambda functions of Lisp.

InUOl:

derivOp[var Symbol]:=

(D[#,,varl&, DE#,var,varl&);

In[I11: '

A[d diffusion][f]

Block[

(dim,dft,disp,fd,sd,,ggp,,dt,,Ops,,

driftop, dispop),

dftdrift[d];dispdispersion[d];

ops derivOp[ symbol[d] ];

fd=ops[[113[f]; sd.ops[[211[f];

dft dft fd;

disp 1/2 disp disp ad;

dt D[f,t];

Simplify[drift + disp + dt]

is

The function A uses an operator notation, As the number of dimensions increase, it

with A [x 3 denoting the generator A for diffu becomes more and more easy to overwhelm

sion Xt. As an example, consider applying the system with expressions that it is unable

this operator to an arbitrary function g(Wt to simplify quickly.

of Brownian motion Wt. Since the symbol With the infinitesimal in hand, the im

associated with the Brownian motion bm is W, plementation of the It~ formula is straight

we denote this function g (w). The following forward, where we again use an operator no

output shows thatAW =g"12. Withinthis tation.

Mathematica session, g has not been defined

and so is treated abstractly. To obtain a more In[l6I:

explicit result, we use "%" to identify the re ito[(~__diffusionl[f_,opts Rule]:

sult at Out[ 121 and substitute Log for g, giv Block[(drift,disptoldSyrnbol,

ing us the infinitesimal for Log[Wj. invRule, ns,invert),

invertitoInvertl.foptell.

In[I2J: options[itol;

oldSymbol symbol[d];

A[bml[g[WII drift A[dIM;

dispD[f,oldSymboll.dispersion[d];

out[121: Print[I1Prior to inversion "];

g111WI Print ["d (", f, 11) (t)

(,drift,")dt C,disp,

2 11) dW (t) n 1 ;

If[Not[invert],

In[I31:= ne f,

% /. g>Log ns=itoSyrnbol /. (opts}

options[itol;

Out[M1:= invRulesFiret

1 [Solve[nsf,oldSymboll);

Print["Inversion rule

2 invRules);

1 2 W __j drift Simplify[drift/.invRules);

disp Simplify[disp/. invRuleal;

We have also experimented with vectors 1 MakeDiffusion[nsdriftdisp) j;

whose elements are diffusions. For example,

vBm is a diffusion in Rd whose elements

W[1]i W[21i ... w[d] areindependent

scalar Brownian motions. With d=2, we find Using the ItOA Formula

that the infinitesimal is the Laplacian

The intermediate printed output from this

1 ( c12g + a2g program is useful in solving stochastic inte

2 axl ax2 grals. For example, what is the value of

t

in[141: sdWs ?

dim 2; p

vBM makeDiffusion[Array[w,dim], Since

Table[0,(i,dim}l, t

identitymatrix[dimll; 1x dx = t212,

In[ISI: 0

A[vBMI [g[w[l],W[211) we might suspect that the stochastic integral

out[ISI: is related to the function x212. If we apply the

(0,2) (2,0) function ito, we obtain

g [w[l],W[2]]+g [w[l],W[211

16

ply the Itt formula directly to the undefined

In [171 functionvts,ti. Setting the option itoIn

ito[bmi[i/2 wA2,itoSymbol>yl vert to False prevents the system from per

Prior to inversion forming the substitution Xt = g(Yt) leading

W 2 1 from (3) to (4). As a result~ the drift and dis

d ()(t) () dt + (W) dw(t) persion of vofs remain functions of the stock

2 2 symbol S rather than introducing some new

Out[M: 1 symbol. Next, at In[211, define the difference

diffu:sion[y,,Sqrt[2]Sqrt[yl,O] between the two representations of V given in

.2 (6) and (7).

In[MI: In[l9J:

ito[bmj[Exp[Wll stock MakeDiffusion[stmu S,sig S];

Out[I8J: bond MakeDiffusion[B, r B, 0];

Y

diffusion[Y, , Y, 1) In[M:

2 vOfs ito[stock][v[S,,t),,

itoInvert>False);

The intermediate result of In[171 implies

that In[211:

t diffvOfs (at stock + bt bond);

wt2 t sdWS In[221:

2 fw

0 roots Solve[(drift[diffl0,

so that the desired stochastic integral is 1 dispersion[diffl=0), fat, bt)l;

t wt2 t The resolution of In[211 requires the rules de

fWsdW, 2 fined earlier for addition (and subtraction)

0 of diffusions. Again, it is not necessary to

The option itoSymbol>y permits the user to resolve the drift and dispersion into func

specify what symbol to use for the new difru tions of, say, the symbol di f f . We want to

sion. The default symbol, as seen in In[181, retain the forms appearing in (6) as func

is Y. An additional option itoInvert can be tions of St and Bt. Finally, to solve for the so

used to suppress the substitution ofg(Yt) for lution in terms of at and bt, apply the builtin

Xt, leaving the new diffusion in the interme function solve to the equations based on the

diate form (3). The second example shows drift and dispersion functions of the differ

that exp(Wt) is the diffusion ence.

dYt = Xt 12) dt + Yt dWt , At this point we need to do some

which is similar to the proces's simulated in Mathematica housekeeping. The result of

the plot shown earlier. solve is a list of rules identifying all solu

tions known to the system. For example, the

roots of the equation x2 1 = 0 shown with

Deriving the BlackScholes &Uression In[231 are returned as a list of lists of substitu

tion rules. Thus, in order to use roots, we

Our implementation of diffusions permits first use root a to define a value for at and

us to derive the solution to the BlackScholes then extract the result from the extraneous

problem symbolically by replicating the ar list. This result is achieved by the somewhat

gument outlined in the introduction. The cryptic command at line 24. For bt, we also

logic is essentially the same, only all of the simplify the resulting expression. The ex

calculations are performed symbolically by plicit simplification of the resulting expres

Mathematica. sion deletes redundant terms and keeps the

Begin by defining the stock and bond dif system from being overwhelmed later in the

fusions, stock and bond. Then solve for the evaluation.

at and bt portfolio sequences. As before, ap

in [231: where dZt = r Zt dt + s & dWt, ZO = x.

Solve[x^2 1 01 In order to compute the expectation in (11),

Out[231: notice that the calculation at line 28 implies

((X > 1), 1X that the diffusion Zt in (11) is

In[241: & = Exp[ (ro2 12) dt + a dWtl

at First[at/.rootal That is, Zt is the exponential of a Gaussian

Out[241: random variable with mean (ro212) and

(1,0) variance A That is, Zt is lognormal. If we

v IS, tl assume the payoff function g(x) = (xP)+, then

In[251: a standard lognormal integration reduces

bt=Simplify[First[bt/.rootsll (11) to

Out[251:

(0,1) 2 2 (2,0) V( (r+0212)(Tt)

2 v [S, tl + S Sig v [S, t 1 Xt) = X 100X /P) + (12)

Offt)112

log(x 1p) + (ra212)(Tt)

P er(Tt) W O(Tty12 j

With the expressions for at and bt in hand, where OW is the standard cumulative nor

we can solve the resulting partial differen mal distribution. The expression (12) for

tial equation via the FeynnianKac theorem. V(xt) is wellknown and has been explored

In order to use the FeynnianKae theorem, we by Miller (1990) using Mathematica.

need to expand the infinitesimal A[Wx,01 in

(8) and identify a diffusion Xt such that (8)

holds. We first set up the PDE assuming that Directions

the payoff function is g, and then use pattern

matching rules to apply the FeynnianKac With this machinery in place, our plans

theorem. include exploring the breadth of input diffu

In[261: sions that can be handled in this manner.

pde Expand[at S + bt B v[S,t11; For example, what other stock and bond dif

fusions are amenable to these operations?

In[271: We also plan to add other methods which use

soln FeynmanKac[pdet g] diffusions to solve partial differential equa

Out[271:= tions, such as the Girsanov transformation.

g[Z,Tt] At some point, we would like to integrate

Ave[diffusion[Z,r Z,z sig,xll[ 1 these techniques into the PDE solving capa

r(Tt) bilities of version 2 of Mathematica. We also

E need better methods for solving expectations

In[281: with respect to diffusions as seen in (11). The

xt=MakeD!ffusion[x,r(s^2)/2,.9,start] solution from the FeynnianKac theorem is

ito [xt 1 [Exp [X] 1 convenient, but still requires considerable

Out[281: work and insight to reach the useful expres

diffusion[Y, Y r, Y sExp[start) sion in (12).

We gave a single example of the use of

In the result Out[271, Ave[(~_diffusioni vector diffusions. Vector models increase

denotes the expected value with respect to the the complexity of the symbol management,

diffusion d. Translated into more standard and we intend to expand our functions in this

notation, the solution at Out[271 is direction.

r (T4)

Wx, 0 = e E [ g(ZTd 1 (11)

18

Refgrence.g

Arnold, L.U974). Stochastic Differential Equations: Theory and Applications. Wiley, New York.

Duffle, D. (1988). Security Markets. Academic Press, New York.

Gray, T.W. and J. Glynn (1991). Exploring Mathematics with Mathematica.

AddisonWesley, Redwood City, CAL

Miller, R. (1990). Computeraided financial analysis: an implementation of the BlackScholes model. Mathematica Journal, 1, 7579.

Steele, M. (1985). MACSYMA as a tool for statisticians. American Statistical

Association Proceedings of the Statistical Computing Section, 14. American

Statistical Association, Washington. Stine, R. A. (1990). Mathematica in time series analysis. American Statistical Association Proceedings of the Statistical Computing Section. American

Statistical Association, Washington.

Wolfram, S. (1991). Mathematica: A System for Doing Mathematics by Computer (2nd Edition). AddisonWesley, Redwood

City, CA.

1 Supported by NSF DMS8812868, ARO

DAAL03.89G0092, AFOSR890301 and

NSAMDA904892034.

2 Equipment and software provided by a

grant ftom Merck, Sharp and Dohme.