Please see PDF version
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: Black Scholes, 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 well defined 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 JCIIX RI
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 well known 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 Black Scholes 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 risk ftee 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 Black Scholes 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)
Black Scholes 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 Feyninan Kac theorem. external file. If one has access to
If A is the infinitesimal of diffusion Xt, the Mathematica, these expressions can be en
Feynnian Kac theorem states that the solu tered directly into the program. Note that
tion of the partial differential equation terminating a statement with a semi colon
suppresses printing the results of that state
AIV(x,t)l rV(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(T t) then create two diffusion objects.
WX,0 = e E[g(XTT) 1 Xt xl. (9)
rn[l]:
For the Black Scholes problem, the payoff <function for the option determines the bound
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
log normal 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 type checking 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 object oriented approach we use a set of "methods" for manipulating diffusions. Although Mathematica does not provide a full set of object oriented programming facilities, we can extend the abstraction by providing accessor methods. For example, the functions
symbol[d diffusion] :=d[[111;
drift[d diffusion] : dH211;
dispersion[d,_diffusion) : d[[311;
initialValue[d diffusion]:=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 data NestList[nextvalue,(0,xO),100];
Here nRand [ x 1 produces a zero mean 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 built in 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,sym symbol[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:
xDiff MakeD!ffusion[X,Sqrt[X],2,0]; yDiff MakeDiffusion[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 Black Scholes 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),
dft drift[d];disp dispersion[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),
invert itoInvertl.foptell.
In[I2J: options[itol;
oldSymbol symbol[d];
A[bml[g[WII drift A[dIM;
disp D[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:= invRules Firet
1 [Solve[ns f,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 diff vOfs (at stock + bt bond);
wt2 t sdWS In[221:
2 fw
0 roots Solve[(drift[diffl 0,
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 built in
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 house keeping. 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 Black Scholes &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 Black Scholes 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[ (r o2 12) dt + a dWtl
at First[at/.rootal That is, Zt is the exponential of a Gaussian
Out[241: random variable with mean (r o212) and
(1,0) variance A That is, Zt is lognormal. If we
v IS, tl assume the payoff function g(x) = (x P)+, then
In[251: a standard log normal integration reduces
bt=Simplify[First[bt/.rootsll (11) to
Out[251:
(0,1) 2 2 (2,0) V( (r+0212)(T t)
2 v [S, tl + S Sig v [S, t 1 Xt) = X 100X /P) + (12)
Off t)112
log(x 1p) + (r a212)(T t)
P e r(T t) W O(T ty12 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 Feynnian Kac theorem. V(xt) is well known and has been explored
In order to use the Feynnian Kae 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 Feynnian Kac 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,T t] 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(T t) 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 Feynnian Kac 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(ZT d 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.
Addison Wesley, Redwood City, CAL
Miller, R. (1990). Computer aided financial analysis: an implementation of the Black Scholes model. Mathematica Journal, 1, 75 79.
Steele, M. (1985). MACSYMA as a tool for statisticians. American Statistical
Association Proceedings of the Statistical Computing Section, 1 4. 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). Addison Wesley, Redwood
City, CA.
1 Supported by NSF DMS 8812868, ARO
DAAL03.89 G 0092, AFOSR 89 0301 and
NSA MDA 904 89 2034.
2 Equipment and software provided by a
grant ftom Merck, Sharp and Dohme.