Mathernatica and Diffusions
J. Michael Steele and Robert A. Stine
9.1 Introduction
A central aim of this chapter is to illustrate how symbolic computing can simplify or eliminate many of the tedious aspects of the stochastic calculus. The package Di f f usion. m included with this book provides a suite of functions for manipulating diffusion models, and individuals with a basic knowledge of Mathematica should be able to use this package to expedite many of the routine calculations of stochastic calculus. After demonstrating the basic features of this package, we give an extensive example that applies the functions of the package to a problem of option pricing.
This application offers a derivation of the well known Black Scholes formula for options pricing. The derivation exploits the idea of a self financing portfolio (Duffie 1988). Our goal is show how Mathematica simplifies the derivation of the central expression in this problem. This task is to be distinguished from that engaged in Miller (1990) which shows a variety ways to describe and use the Black-Scholes formula.
In the next section we give a brief overview of the central idea of diffusions as
required by options pricing theory~ and we introduce the notation that is needed,
This section is indeed very brief, and readers unfamiliar with diffusions should
consider a text like that of Arnold (1974) or Duffy (1988) for the missing details
and intuition. Section 9.3 introduces the new Mathematica functions that "0
be needed in Section 9.4 where we present a Mathematica derivation of the
Black Scholes formula. Section 9.5 demands more knowledge of Mathematiffl
since it describes issues of implementation and includes further examples of tile
methods used in the extended example. Of
This package uses features found in Version 2 of Mathematica. Since sonie these features are not available in earlier versions of the software, problems Will occur if one attempts to use the package without having upgraded to at least Version 2.0.
9. Mathematica and Diffusions 193
9.2 Review of Diffusions and lto''s Formula
Diffusions form a class of stochastic processes that allow one to apply many of the modeling ideas of differential equations to the study of random phenomena. In simplest terms, a diffusion X_t is a Markov process in continuous time t > 0 that has continuous sample paths. A key feature of diffusions is that there exist two functions 1A and a which together with an initial value Xo completely characterize each diffusion. Moreover, in many applications ft and a have physical or financial interpretations.
For this article we will assume that the scalar valued functionsp and a are locally bounded, and let Wt denote the standard Weiner process. Under these conditions, a fundamental result of the It6 calculus is that there exists a welldefined process Xt having the suggestive representation
t t
xt = X0 + p (X, s) ds + a (X, s)dW,.
10 10
The intuition behind this representation for Xt is that the function p determines the instantaneous drift of Xt whereas a controls its disperson or variability. If, as in later examples, Xt denotes the price of a stock at time t, then p and (T may be interpreted as the rate of return and the instantaneous risk.
Since the process Wt is not of bounded variation, the second integral on the right hand side of (1) cannot be interpreted naively. For example, suppose that we were to attempt to evaluate fit W, dW, as the limit of an approximating sum over partitions of the interval [0, t],
w~,i, (Wti Wti '), 0 < 4 1 < si < 4 :5
Since Wt is not of bounded variation, the value for the integral suggested by such approximation turns out to depend upon how we choose to locate si in the interval [ti l,ti]. Throughout this chapter, we adopt the convention si = ti j. The resulting stochastic integral has a unique solution known as the It6 integral. A very important reason for choosing si = ti 1 is that with this choice a stochastic integral with respect to Wt is always a martingale.
One of the most common and convenient ways to express a diffusion is to use a notation that is analogous to that of differential equations. Rather than use the stochastic integral to represent Xt as in (1), one often uses a shorthand notation:
dXt 1A (Xt, t) dt + a (Xt, t) dWt. (2)
This compact description of the process Xt resembles a differential equation, but it must be interpreted appropriately as just shorthand for (1).
An important property of diffusions is their behavior under smooth transformations. Perhaps the central result in the theory of diffusions is that a smooth function of a diffusion produces another diffusion. Moreover, an explicit formula due to It6 identifies the new process by giving expressions for the drift and disperson functions. Suppose that the function f maps (R, C) R,
194 J. Michael Steele and Robert A. Stine
and consider applying this function to a diffusion, say Yt = f (Xt, t). Under modest conditions on f, the process Yt is also a diffusion. One set of sufficient conditions is to require continuous first and second derivatives of f in its first argument and a single continuous derivative in the second. Denote these derivatives by f., = Of (x, t) lax, fxx = O'f (x, t) lax', and ft Of (X, t) lat.
Under these restrictions on f, It6's formula guarantees that Yt is a diffusion and gives the functions that characterize this new diffusion. Specifically, the stochastic differential for Yt = f (Xt, t) is
d Yt = {ft (Xt, t) + p (Xt, t) f, (Xt, t) + f.x 0,2 (Xt, t) /2}dt + {f. (Xt, t)o. (Xt, t)} dWt (3) This expression is not in the canonical form of (2) since the drift and dispersion functions of Yt are functions of Xt rather than Yt, but this problem is easily remedied. If the transformation f has an inverse g so that
Y = f (X, t) t=> X = M' t), then we can use this inverse to obtain the desired form. Substitution in (3) gives
dIYt = {ft (g (Yt, t), t) + p (g (Yt, t), t) fx (g (Yt, t), t) +
2(g(y
f., X 0, t, t), t) / 2 } dt + {f. (g (Yt, t), t) a (g (Yt, t), t) } dWt
So indeed,
dYt (Yt, t) dt + & (Yt, t) dWt (4)
for the indicated values of Tt and a. The complexity of (4) suggests that we can avoid some tedium if we put Mathematica to work managing the transformation of the pair {IA (x, t), or (x, t) } into {TA (x, t), 5~ (x, t) }
It6's formula (3) can be made more revealing if one also introduces the notion of the infinitesimal generator of the diffusion Xt. The infinitesimal generator of a diffusion measures the instantaneous rate of change in the expected value Of a transformation of a diffusion as a function of its starting value. Under soille natural conditions on the function f , we can define the infinitesimal generator as
Axf (x) limt,o Ex f (Xt, t) f (x, 0) t
where E,, denotes the expectation conditional on the starting value Xo = x. This limit turns out to be closely related to the drift expression in It6,s formula, and one can show under mild conditions that Ax is the differential operator that has the explicit form
AX + 1A (X, t) a + 1 a (X, t) a2
Yt jx 2 aX2
By using the infinitesimal generator associated with the diffusion Xt, we obtaill a more compact form of It6's formula,
dYt = # (Xt, t) = Axf dt + fx (Xt, t)a(Xt, t)dWt.
9. Mathematica and Diffusions
9.3 Basic Mathernatica Operations
9.3.1 Introduction to the Package
This section introduces basic functions that permit us to build and manipulate diffusions. In order to follow the operations described here, one must first import the package that defines the needed functions. The package itself is located in the file named "Dif f usion. m" which should be placed in a directory which is easily, if not automatically, searched. The following command imports the package.
In[l]:= << Diffusion'
The names of most of the functions in this package begin "'Diffusion", though some also have shorter, more convenient names. This convention allows for an easy listing of the available functions via the built in Names function.
In121:= Names["Diffusion*"]
Out[21:= {DiffusionA, DiffusionChangeSymbol, DiffusionDispersion, DiffusionDrift, DiffusionExpand, DiffusionExpandRules,
DiffusionFeynmanKac, DiffusionInitialValue,
DiffusionIto, DiffusionLabel, DiffusionMake,
DiffusionPrint, DiffusionSimulate, DiffusionSymbol, DiffusionWeinerProcess, DiffusionWeinerProcessMake, DiffusionWeinerProcessSymbol)
A brief synopsis of each function is available via the standard Mathematica request for help; one precedes the name of the function of interest with a question mark. For example, the following command reveals the syntax and a brief description of the function which prints diffusions.
In[31:= ?DiffusionPrint
DiffusionPrint[d] prints a formatted version of the diffusion object d using subscripts and symbolic names.
The next section describes how to build the diffusion objects that this function requires as arguments.
9.3.2 Building a Diffusion
The first order of business in building a diffusion is to specify a Weiner process. In this and many other examples, we begin by calling the Clear function. This function removes any value that might already be bound to the named symbol, as might occur when experimenting with Mathematica outside of the range of commands shown here. Thus, we first clear the name W which we intend to use to represent a Wiener process.
196 J. Micha& Steele and Robert A. Stine
In[l]:= Clear[W]
The next command associates the name W with a Weiner process. Until we say otherwise, the name W will denote a Weiner process to the system.
In[2]:= WeinerProcessMake[W]
We can now use this Weiner process to define more complex diffusions. To build each diffusion, our convention requires a symbolic identifier (or name) for the diffusion and for the underlying Weiner process, two expressions (the drift ft and the dispersion a), and an initial value. As an example, we first consider the lognormal diffusion. This process is used extensively in financial modeling and it will play a central role in Section 9.4 when we derive the BlackScholes formula. In the classical notation of the stochastic calculus, a lognormal diffusion can be specified by
dXt = aXtdt + 0XtdWt, XO = v,
for scalars a and 0 > 0. To represent this diffusion in Mathematica, use the command Dif fusionMake. Again, we clear the name X before we make the new process.
In[31:= ClearM In[41:= DiffusionMake[X, W, alpha X, beta X, v]
To test the success of this command, we recall that typing the name of any object in Mathematica reveals the value of that object. In the case of a diffusion, the value is a list that holds the definition of the process.
In[S1:= X
Out:[51:= diffusion[x, W, alpha X, beta X, v]
While revealing to the Mathematica programmer, this form is probably not very appealing to one more familiar with the standard mathematical notation. The function Dif fusionPrint ameliorates this problem by producing a more familiar rendering of a diffusion, complete with subscripting. Here are examples for the two diffusions created so far. The printed form for a Weiner. process is pretty simple.
In[6]:= DiffusionPrint[W] W = Weiner Process with scale 1 t
The format in general for printing a diffusion matches the notation of equa
tion (2).
In[71:= DiffusionPrint[X]
dX = alpha X dt + beta X dW X = v
t t t t 0
9.3.3 Simulating Diffusions
In order to experiment with a diffusion model, it is often useful to simulate realizations of the process. The function Dif fusionSimulate generates sill"
9. Mathematica and Diffusions 197
ulated realizations by building a discrete time approximation to the continuous~ time process. The method implemented in this software is a little naive it resembles Euler's method for approximating solutions of a differential equationbut it functions quite usefully in most instances. To illustrate, we begin by asking Mathematica to reveal its summary of Dif fusionSimulate.
In[l]:= ?DiffusionSimulate
DiffusionSimulate[d diffusion, n, dt] simulates a diffusion. It returns a list of n pairs
{ {i*dt,x[i*dt]} } of the diffusion d. Each realized
series is an independent realization. Responds with an
error if the diffusion contains symbolic paramters.
As a simple example of the use of this function, we will simulate 100 values of a Weiner process at times {O, 0.01, . . . , 0.99}. Assuming that we are simulating a diffusion labelled Xt, the simulated realization begins with the pair {O, Xo}
and adds n 1 pairs of the form {i dt, Xidt}, i = 1, 2,. . ., n 1, separated in time by the chosen step size dt which is the last argument to the function. In this example, the semicolon at the end of the command suppresses the printing of the complete simulated realization. The function Short reveals the start and end of the list, and the bracketed number indicates that 97 of the items in the list are hidden.
In[21:= simW = DiffusionSimulate[W,100,0.011; In [3] : = Short [simW]
Out131:= ~{0, 0}, {0.01, 0.0984664}, <<97>>, {0.99, 1.55817}}
The list structure of the simulated diffusion makes it quite easy to have Mathematica plot this artificial realization versus time.
In[41:= ListPlot[simW, Plotioined True]
. AA
o.2 0.4 0.6 0.8 1
0.5
Out[41:= Graphics
198 J. Michael Steele and Robert A. Sti
Simulations of other diffusions can be obtained just as easily. For example, suppose that we wished to plot the simulation of the diffusion
dYt = Ytdt + MWt, Yo = 10.
First we build the diffusion using Dif fusionMake.
In[51:= DiffusioriMake[Y,W, Y,Y,10]
Next we generate the partial realization using Dif fusionsimulate, and again we note that Short reveals just the extremes of the list. Since the time spacing is not specified in the call to Dif fusionSimulate, this program sets the spacing to a default value of 0.01.
1n161:= simY = DiffusionSimulate[Y,100];
Short[simY]
Out:[61:= {{O, 10.}, {0.01, 9.99454}, <<97>>,{0.99, 37.4779}}
In[71:= ListPlot[simY,
Plotioined True]
30
25
20
15
10
0.4 0.6 0.8
Out:[71:= Graphics
The function RandomSeed makes it possible to generate correlated realizations. By default, each realization of a diffusion is independent of other realizations. This structure may not be desired for many problems. As an example, consider the relationship of the Weiner process realization and the diffusion for the process dYt = Ytdt + YtdWt considered above. Independent realizatiOlls fail to capture the relationship between the series and suggest little relationshiP between the two. After all, these simulated realizations are independent. I Jere is a plot of Yt on Wt with the plots joined in time sequence order.
In181:= ListPlot[Table[{simW[[i,2]], simY[[i,2111,(i,1,100)lf
AxesLabel {"W[t]", Ovy[tIVY),
Plotioined True],
9. Mathematica and Diffusions 199
22
5
.5
%, W 1 t 1
1' Po 0 .5
Out[S1:= Graphics
By explicitly setting the random seed used by Mathematica, we force the real~ ization of Yt to be based on the same random values used to simulate the Weiner process Wt. This connection gives the realizations we would have expected.
In[91:= SeedRandom[732712];
simW = DiffusionSimulate[W,100];
SeedRandom[732712);
simY = DiffusionSimulate[Y,100];
The plot of Yt on Wt now exhibits the strong relationship between the two simulated re~tlizations that is missing without matching the simulation seeds.
In[M1:= ListPlot[Table[(simW[[i,21], simY[[i,2]]},{i,1,100)], AxesLabel >("W[t]", Vvy[t111), PlotJoined > True] Y[t]
25
20
POW ,
W[t]
0.2 0.4 0.6
Out[10]:= Graphics
200 J. Michael Steele and Robert A. Stine
9.4 Deriving the Black Scholes Formula
To illustrate the use of these tools in a problem of considerable historical interest, we will use our tools to derive the Black Scholes option pricing formula. Our approach parallels the development in Duffle (1988) and uses the notion of a selffinancing trading strategy. Following the tradition in elementary option pricing theor)~ we will ignore the effects of transaction costs and dividends. Readers who are curious for further details and intuition should consider Duffle (1988) for a more refined discussion than space permits here.
The model begins with a simple market that consists of two investments, a stock and a bond. The bond is taken to be a risk free asset with rate of return p, so in our differential notation the bond price Bt at time t obeys dBt = PBtdt. The stock is assumed to have rate of return IA as well as some risk, so we can model the stock price St as a diffusion for which dSt = MStdt + oStdWt with the scalar a > 0. We can set up these processes and clear the needed symbols as follows:
In[l]:= Clear[W, S,so, B,BO, V, mu,sigma,rho, g,PI
WeinerProcessMake[W]
DiffusioriMake[S, W, mu S, sigma S, SO]
DiffusionMake[B, W, rho B, 0, BO]
Before continuing, we should perhaps view our processes in more conventional form by getting the printed version of each of the three diffusions.
In[21:= DiffusionPrint[W] DiffusionPrint[S] DiffusionPrint[B]
W Weiner Process with scale 1
t
dS mu S dt + sigma S dW S= SO
t t t t 0
dB rho B dt B = BO
t t 0
With these processes as building blocks, we can define the equations that permit us to give the explicit value of the European option. Suppose that there exists a function V (s, t) such that V (St, t) is the value at time t, 0 :5 t :5 T, of an option to purchase the stock modeled by St at a terminal time T at the strike price P. Our goal is to find an expression for this function in temls of p, a, P, T, and t. Since the value of the option is a function of the difft" sion St, It6's formula gives an expression for V(St, t). Even though we d, not explicitly know V(St, t), we can use the Ito function to identify this nely diffusion symbolically in terms of derivatives of V. In this case, we want tO avoid putting the diffusion into the canonical form (4) since doing so Would conceal how V(St, t) depends upon St. The use of the itoinvert opti011 avoids the inversion to the canonical diffusion form, thereby retaining the
1 1
9. Mathematica and Diffusions 201
stock symbol S in the expressions for the drift and dispersion of the new diffusion.
Inj31:= Ito[V[S,t], ItoInvert >False];
(011) (110)
dV[S, t] = (V [S, t] + MU S V [S, tl +
2 2 (2,0)
sigma S V [S, tl (110)
) dt + (sigma S V [S, tl) dW
2 t
Now equate the value of the option to that of a self financing trading portfolio. That is, assume that the value V (St, t) of this stock option can be reproduced by a portfolio consisting of at shares of stock and bt shares of bond, V (St, t) = at St + bt Bt. Here, of course, we are assuming that at and bt are both stochastic processes. An argument for the existence of such a portfolio in this problem appears in Duffy (1988). The matching of the value of the option to that of a portfolio gives a second expression for V in addition to that from It6's formula. If we equate these two expressions for V(St, t), we can solve for at and bt. The manipulations that support the elimination of at and bt are very easy in Mathematica since the diffusion package defines some simple algebraic operations for diffusions.
By default, diffusions retain their symbolic form in algebraic expressions. For example, if we enter a sum of two diffusions., the sum is retained and no attempt is made to combine the drift of one with that of the other.
In[41:= at S + bt B
Out[41:= bt B + at S
This behavior is consistent with the way Mathematica handles many other symbolic expressions, such as the way a product is not expanded unless the user makes an explicit request:
In[51:= (a + b) (c + d)
out[51:= (a + b) (c + d) in[61:= Expand[%]
Out[61:= a c + b c + a d + b d
The analogous behaviour is needed in the algebra of diffusions. In order to equate the two expressions for V (St, t), we require the drift and dispersion of a diffusion that is the sum of two diffusions. The function DiffusionExpand combines several diffusions, though the resulting "diffusion" is not of the canonical form. Since we are only interested in the drift and dispersion, the absence of the canonical form is not a problem.
In[71:= diff = DiffusionExpand[V (at S + bt BM
202 J. Michael Steele and Robert A. Stine
Out[71:= (011)
cliffusion[$3, W, (bt rho B) at mu S + (V) [S, tl
2 2 (2,0)
(110) sigma S (V) [S, tl
+ MU S (V) IS, tl +
2
(110)
(at sigma S) + sigma S (V) [S, t11
(bt BO) at SO + V[SO, tl 1
In order for the difference V (St, t) (at St + bt Bt) to be zero, both the drift and dispersion of this new diffusion must be identically zero. This gives two equations in two unknowns, and thus the number of stock and bond shares in the matching portfolio, at and bt. The built in function Solve gives a set of rules that define the solution of the system of two equations.
In[B1:= roots = Solve[(DiffusionDrift[diff]==0,
DiffusionDispersion[diff]==0}, (at,bt}]
Out[el:= (1,0)
Hat > (V) [S, tI, bt
(110)
MU S (V) [S, t]
rho B
(011) (110)
(2 (V) [S, t] + 2 mu S (V) [S, tl +
2 2 (2,0)
sigma S (V) [S, tl) / (2 rho B)}}
For convenience, we next extract the values of at and bt implied by these rules and assign them appropriately. The additional simplification eases later rnanipulations and seems unavoidable in some implementations of Mathematica. The function First in the next two expressions extracts the solution from the
single element list in which it is embedded.
In[91:= at = simplify[ First[at roots]
Out191:= (1,0)
(V) [S, ti
In1101:= bt = Simplify[ First[bt roots]
Out[lol:= (011) 2 2 (2,0),
2 (V) [s, ti + sigma s (V) [S, tl
2 rho B
In more conventional notation,
V
t (St, t) + 10,2 St2 V
at = V. (St, t), bt 2 (St, t)
2pBt
9. Mathematica and Diffusions 203
where V
., denotes the partial derivative of V(x,u) with respect to its first argument, and Vxx denotes the second partial derivative in the first argument. Similarly, Vt is the first partial in the second argument.
To find a partial differential equation for the value of the option, we substitute these expressions for at and bt back into the relation V (St, t) = at St + bt Bt. The use of the function Expand assures that the function Coefficient extracts the proper term.
In[111:= pde = Expand[ at S + bt B V[S,t]
(011)
(V) [S, tl (110)
V[S, t] + + S (V) IS, tl +
rho
2 2 (2,0)
sigma S (V) [S, tl
2 rho
To set things up for the next step, normalize this PDE so that the coefficient of Vt is 1.
In[121:= pde = Expand[ pde / Coefficient[pde, D[V[S,t],t]]
Out[121 (011) (110)
(rho V[S, t]) + (V) [S, t] + rho S (V) [S, tl
2 2 (2,0)
sigma S (V) [S, tl
2
These equations and the boundary conditions discussed shortly are sufficient to determine V, thus solving the option pricing problem. The problem now faced is the purely mathematical one of solving our PDE. A variety of tools exist for solving PDE's that arise in the application of diffusions, and one of the most powerful is based upon the Feymann Kac theorem. For our purposes, the Feynman Kac theorem expresses the solution of a certain second order PDE as an expectation with respect to a related diffusion. The first question one has to resolve before applying this method is whether the PDE of interest is of the appropriate type. The second issue is to make an explicit correspondence between one's PDE and the form of the Feynman Kac result. The function FeyrmanKac performs both of these tasks. As a side effect, it also prints out the terms used in its matching using the notation of Duffle (1988).
The use of this function also requires identifying boundary conditions. Let g (x) denote the payout function for the option. For example, the payout function for a European option with exercise price P is the piecewise linear function g(x) = (x P)+ where (x P)+ = x P if x > P and is zero otherwise. The payout function determines the needed boundary condition, V(xJ) = g(x). Here we apply the FeynmanKac function, using the symbol g to denote
204 J. Michael Steele and Robert A. Stine
the payout function. In the following output, the symbol Ave stands for the expected value operator since the symbol E denotes the base of the natural log e in Mathematica.
In [131 soln = FeyrimanKac [pde, g]
f = V; rho = rho; u = 0
dX = rho X dt + sigma X dW X= X
t t t t 0
Out[131:= g[X[ t + T1]
Ave[ 1
rho ( t + T)
E
The ".'solution" given by the Feynmart Kac theorem is rather abstract and the task of rendering it concrete is not always easy The result of this function indicates that the solution of our PDE is the expected discounted payout
e p(T t)E.q(XT t),
where Xt is the diffusion that satisfies dXt = pXtdt + aXtdWt with initial value XO = x and Weiner process Wt. The process Xt is just the familiar lognormal diffusion, as confirmed by use of 1t6's formula.
We will next confirm that the exponential of the normal diffusion dYt a dt + b dWt is indeed a lognormal diffusion, and we use Mathematica to make the required identification of coefficients in both drift and dispersion.
In[141:= Clear[Y];
DiffusionMake[Y,W, a,b,O];
Ito [Y] [EXP [Y1 1
2 Y
Y (2 a + b ) E Y
dE dt + (b E ) dW
2 t
Inversion rule ... {Y > Log[Z]}
Out[141:= 2
(2 a + b ) Z
diffusion[Z, W, b Z, 11
2
Clearly, b corresponds to a. Equating the drift coefficients (2a + P 2) p shows that a = p (s2/2).
In1151:= a = Simplify[a/.First[ Solve[(2 a+sigma^2)/2 == rho, an]
Out[I51:= 2
sigma
rho 2
This calculation implies that the diffusion Xt in the solution from the Feynrilar,Kac theorem is the exponential of a normal diffusion with constant drift P
Mathematica and Diffusions 205
(0,2 /2) and dispersion a. Recalling that Xo = x, at any time t > 0 Xt satisfies
j~2)t+0.Wt,
Xt = xe(p 2 _
Hence for any t, Xt has the same distribution as xe(P (,2 12))t+,'I'_Tz where Z is a standard normal random variable with mean zero and variance one. Finally we see that the value of the option at time t = 0 is
,T p pT 00 (P ~2 )T+v/Taz _Z2 /2
V(X,O) = e Eg(XT) = fl g(xe 2 )e dz.
For the European option with exercise price P, we have g (x) = (x P)+ and this integral becomes
6 pT oc )T+.VT az P)e_,2/2
= f (xe (P '2 dz,
,/27r ZO
where zo solves
xe(p jz 2 )T+vrT azo p.
2
Calculation of this integral is somewhat tedious, but its evaluation makes for a nice illustration of the integral solving capabilities of Mathematica. To make the coding a little more modular, we first define a function h(z) xe(P (c' 2 11))T+,1 T az P and use it to locate the lower bound zo.
In[161:= h[z_l X EXp[(rho sigma^2/2.)T+Sqrt[T] sigma z] P
rule = Solve[ h[z]==O, z];
zo = Simplify[ First[z /. rule] 1
Out[161:= 2 p
(rho T) + 0.5 Sigma T + Log[ ]
X
sigma Sqrt[T1
To simplify the input further, we define the Gaussian kernel gauss (z)=e7
In[171:= gauss[z_l := E^( (z^2)/2)/Sqrt[2 Pil
Given these definitions of the functions h and gauss, it is quite simple to describe the needed integral. We call the integration function with specific limits a and b rather than symbolic infinite limits. The integration routines seem to behave more robustly in this application with these limits rather than infinite limits.
In[M1:= Clear[a,b]
intab = Integrate[E^( rho T) h[z] gauss[z], {z,a,b)l
Out[M1:= a b
P Erf[ 1 P Erf[ 1
Sqrt[21 Sqrt[2]
rho T rho T
2 E 2 E
206 J. Michael Steele and Robert A. Stine
Out [181 (con t
2 2
(rho T) + (sigma T)/2 + (rho 0.5 sigma ) T
(E X
a sigma Sqrt[T]
Erf[ 2 +
Sqrt[2]
2 2
(rho T) + (sigma T)/2 + (rho 0.5 sigma T
(E X
b sigma Sqrt [T]
Erf[ 2
Sqrt[2]
Now we can apply two substitution rules that specify the range of integration, a = ZO and b = oc, and we have the desired integral.
In[191:= int = intab /. (a >zO, b >lnfinity};
Sitaplify[int]
Out[191:= 2 p
(rho T) 0.5 sigma T + Log[ ]
X
x Erf[
X Sqrt[2] sigma Sqrt[T1
rho T 2 2
2 E
2 p
(rho T) + 0.5 sigma T + Log[ ]
X
P Erf[ 1
Sqrt[21 sigma Sqrt[T1
rho T
2 E
This expression involves the error function defined in Mathematica as erf(z) =
fz e 1 2 dx. We get a more familiar result by using the equivalence erf(z) ~ 2 4) (v/'2 z) 1, where ~D (z) denotes the cumulative standard normal distribution, (D (Z) = f (e x 2 / 2) / v/"2 7r) dx.
In[201:= bs Simplify[int /. Erf[i~_1 >
(2 NormalCDF[Sqrt[2]x] 1)]
Out[201:= 2 p
(rho T) 0.5 s T + Log[ ]
p X
( ) + x x Normal= +
rho T s Sqrt[T)
E
9. Mathematica and Diffusions 207
Out[201(cont.)
2
(rho T) + 0.5 s T + Log[ ]
X
P NormalCDF[ s Sqrt[T]
rho T
E
We can express this result in yet more familiar form by using a rule that expresses the relationship (P (x) = 1 P( x).
In[2l]:= Simplify[ bs/.NormalCDF[x ] >1 NormalCDF[ x]]
Out:[211:= 2 p
(rho T) 0. 5 s T + Log [ ]
X
x NormalCDF[ ( s Sqrt[T]
2 p
(rho T) + 0.5 s T + Log[ ]
X
P NormalCDF[ ( s Sqrt[T]
rho T
E
This is the Black Scholes formula for pricing the European option (Duffle, 1988, p. 239). Miller (1990) discusses this function at length and shows various manipulations of it using Mathematica.
9.5 More Mathematica Details and Examples
The material of this section focusses on some issues that can be omitted at a first reading, but that might prove useful for the user who wishes to extend the methods or examples. The discussion begins with the underlying data structure used to represent a diffusion. Many of the ideas come from objectoriented programming. We next consider in some detail the use of our functions that produce the infinitesimal and It6's formula. Further programming details appear as comments embedded in the code of the package itself.
The examples in this section use both of the Mathemaffica diffusion objects built in Section 9.3. In case these are not available, the following command rebuilds the Weiner process Wt and the associated lognormal diffusion Xt.
in[l]:= Clear[X,W1;
WeinerProcessMake[WI;
DiffusionMake[x,w,alpha X, beta X, v];
The data structure we use to represent a diffusion parallels the notation of Section 9.2. To represent an arbitrary diffusion, one must describe each of the
208 J. Michael Steele and Mobert A. Sth
distinguishing components: drift, dispersion, Weiner process, and initial value In addition, we need to supply a symbol that names the diffusion itself. As ~ result, the list that we use to represent a diffusion has these five items: name oi the diffusion, name of Weiner process, drift expression, dispersion expression, and initial value. This data structure is particularly simple for a Weiner process, The first two symbols which are the names of the diffusion and its underlyinS Weiner process match for a Weiner process. The drift and dispersion are thE constants 0 and 1, respectively. The matching of the leading symbols identifies in the software the presence of a Weiner process.
Several accessor functions permit one to extract components of a diffusion without requiring detailed knowledge of its data structure. Although these functions simply index the list that holds the components of the diffusion, use of these accessors frees one from having to remember the arrangement of the list. Such abstraction offers the opportunity to change the data structure at a later point without having to rewrite code that uses accessor functions. Here are a few examples of the accessor functions.
In[21:= DiffusionWeinerProcessSymbol[Xl
Out[21:= W
Notice that this function results in the symbol associated with the Weiner process rather than the Weiner process itself. A different accessor extracts the process.
In[31:= DiffusionWeinerProcess[X]
Out[31:= diffusion[W, W, 0, 1, 0]
The next three examples extract the remaining components of the diffusion.
In[41:= DiffusionDrift[X]
Out[41:= alpha X
1n151:= DiffusionDispersion[Xl
Out[51:= beta X
In[61:= DiffusionInitialValue[Xl
Out[61:= v
It is important to notice that the results of some of these functions include the symbol X which represents the diffusion. A sleight of hand is needed to obtaill this behavior, and all is not quite as it appears. The built in function Ful1F01" reveals that the results of both Dif fusionDrift and Dif fusionDisper5100 contain so called "held expressions" that one must use in order to keep the system from trying to evaluate the diffusion symbols that appear in the drift and dispersion expressions.
in[71:= DiffusionDrift[xl // FullForm
Out[71:= Times[alpha, HoldForm[X11
1 1
9. Mathematica and Diffusions 209
Were it not for holding the evaluation of the symbol x in this expression, Mathematica would descend into an endless recursion, continually substituting the list representing the diffusion each time it encountered the symbol X in the list. Further discussion of this recursion appears in Steele and Stine (1991).
We have found that it is most easy to manipulate the drift and dispersion functions as expressions rather than Mathematica functions. Still, occasions arise when one wants a function. For example, one might want to differentiate or plot the drift. The accessors normally extract these components in a manner that preserves their symbolic content. That is, they return an expression which includes the symbol for the diffusion. When a function is desired, the accessors Dif fusionDrift and DiffusionDispersion have an optional argument that forces the output to be returned as a function. The returned function has two arguments, the first for the diffusion and the second for time, as in U (x, t) and o,(x, t) respectively.
As example, consider generating a plot of the drift from the lognormal diffusion. First, extract the drift as a Mathematica function and give it a suggestive name.
In[B1:= mu = DiffusionDrift[X,Function]
Out[B1:= Function[{x$, t$}, alpha x$1
We can treat this function just like any other, differentiating or plotting as we choose.
In [_91 : = D [mu [x, t] , xl
Out[91:= alpha
Of course, if we expect to plot the drift function, we have to make sure that all of the symbolic terms have a value. Here we set ce 2 so that fi(x, t) 2x.
In[101:= Plot[ mulx,tl/.alpha >2, {x,0,5)
10
8.
6
4
2
2 3 4 5
Out[I01:= Graphics
210 J. Michael Steele and Robert A. Stin(
In general, we would need to plot the drift ti(x, t) as a surface over the plane, bu, in this and many other common circumstances this elaborate plot is not needed
As we noted in our review of 1t6's formula, the infinitesimal generator oi a diffusion process has an intimate relationship to the stochastic differentia': representation of the diffusion. The argument given to the function which findE the infinitesimal is an expression involving a diffusion. For example, the nexl example determines the infinitesimal of the process defined by an arbitrary function g of Wt. We see that the infinitesimal is half of the second derivative of g.
InIll]:= Afg[W] I
Out [1l]:= g,' [W]
For the lognormal diffusion Xt, the result is somewhat more complex, but the syntax of the commands is the same. We saw this expression in Section 9.4 in the option pricing problem.
In[121:= A[g[X] I
Out [121:= 2 2
beta X g" [X]
alpha X g'[Xl +
2
The function associated with 1t6's formula provided in the accompanying package is more potent than the other functions of the package. The function it 0 (or more elaborately, Dif fusionito) builds the new diffusion associated with the input expression which again is an expression which includes a diffusion. The program also assigns a default name to the new diffusion unless a new name is chosen.
As an example we show how to use It6's formula to find the diffusion associated with half the square of a Weiner process, Yt = Wt/2. The optional argument sets the name of the new diffusion to be the symbol Y which has been cleared of any prior value.
In [131 Clear [Y]
Ito[l/2 W2, ItoSymbol >Y]
d dt + (W) dW
2 2 t
Inversion rule... {W > Sqrt[2] Sqrt[Y]}
Out:[131:= 1
diffusion[Y, W, , Sqrt[2] Sqrt[Y], 0]
2
Two pieces of intermediate output precede the final result in this example. The first portion of the output shows the diffusion associated with the transformation
9. Mathematica and Diffusions 211
in the form of (3) before conversion to canonical form. This portion of the output is occasionally useful in solving various stochastic integrals. In this example, the output expression dWt2/2 = (l/2)dt + WtdWt suggests the solution to a stochastic integral,
t
WdW, = Wt2/2 1/2
Of course, one would have to know to look at It6's formula applied to Wt' in
order to find fot' WdW, in this way. However, since f,'xdx = t2 /2, this is not
0
such a bad place to start looking for an answer.
The second piece of intermediate output in the example beginning "Inversion rule. . ." (above) gives the inverse transformation used to convert the result of 1t6's formula (which is a function of Wt) into a function of Yt. This is the inverse function g described in the introductory review.
Here is another example. This example shows that the exponential of a Weiner process is a lognormal diffusion. Notice as always that it is important to begin with an unbound symbol to use for the new diffusion.
in [141 Clear [Y1
Ito[Exp[W],ItoSymbol >Y]
W
W E W
dE ( ) dt + (E ) dW
2 t
Inversion rule ... {W > Log[Y1}
Solve::ifun:
Warning: Inverse functions are being used by Solve, so some solutions may not be found.
Out[141:= Y
diffusion[Y, W, , Y, 11
2
The warning message in this example will often suggest a problem in the inversion process. In this example the inversion via logarithms works since Yt = ewt, implies Wt = log, Yt. Were the process Wt complex, for example, such inversion might not be appropriate. Since Mathematica does not assume Wt is real valued, it displays a warning.
9.6 Summary and Concluding Remarks
This chapter began by reviewing some of the basic notions of the theory of diffusions such as the local drift g, the local dispersion a, and the specification of a diffusion process though the formalism of stochastic integraton. We then illustrated how the fundamental formula of 1t6 for functions of diffusions can lead to calculations that can be profitably performed by symbolic methods.
1
212 J. Michael Steeie and Robert A. StinE
Much of the chapter served to illustrate how our package for stochastic calculuss can be applied.
The problem chosen for this illustration was that of the derivation of the famous Black Scholes formula for pricing of European options. This problem offers a natural first hurdle for any stochastic calculus package, and the successful clear~ ing of that hurdle provides honest reassurance that the package implementation is rich enough and robust enough to engage problems of genuine interest. Since the practical and theoretical success of the Black Scholes theory has served as a central motivation for much of the spreading of the technology of stochastic calculus, any symbolic system for stochastic calculus must sooner or later meet the challenge offered by this theory.
The implementational details given in Section 9.5 provide some of the design criteria that were used in the development of our stochastic calculus package, but it should be clear that even this second formulation (cf. Steele and Stine 1991) leaves room for further work. One immediate concern is that any serious stochastic calculus package needs to be able to handle the multivariate versions of It6's formula. The most recent version of our package does handle multivariate diffusion, and, in fact, the current distribution of the package owes much of its design to the demands of multivariate diffusions. Because we chose to provide the easiest possible introduction to symbolic stochastic calculus, we have kept this chapter and the accompanying software focused on one dimensional diffusions. Even these diffusions offer many important applications, and for those with pressing need for mutivariate tools the best recourse is to go directly to the package internals.
One natural development that has not been addressed in any version of the package is that of integrating the tools of stochastic differential equations with the differential equation solving features of Mathematica. This issue is more complex than it might seem at first thought, and at the very minimum one needs to come to grips with the fact that even simple problems in stochastic differential equations can lead to partial differential equations that test the limits of modern theory. One twist that could evolve to be of considerable interest is that of reversing the natural course of investigation and using the stochastic calculus to assist in the development of symbolic tools for PDE's. This may seem an unlikely turn of events, but the use of the Feynman Kac formula here to solve the key PDE provides a hint of the practicality of the process.
Symbolic computation has been available in some special environments for a good number of years, but the widespread use of symbolic computation is still rather new. The number of users of symbolic computation is growing at more than 50% per year, and the community of sophisticated users is already many times larger than one could have imagined just a few years ago. The vital importance of diffusion processes in many areas of scientific and social modelling teamed with the explosive growth in symbolic computation virtually guarantee that symbolic stochastic calculus can expect to have many subsequent developments of far reaching influence.
9. Mathematica and Diffusions 213
9.7 References
Amold, L. (1974). Stochastic Differential Equations: Theory and Applications. Wiley, New York.
Duffie, D. (1988). Security Markets. Academic Press, New York.
Steele, J. M. and R. A. Stine (1991). "Applications of Mathematica to the stochastic calculus." In American Statistical Association, Proceedings of the Statistical Computing Section., 11 19, Amercian Statistical Association, Washington, D.C.
Miller, R. (1990). "Computer aided financial analysis: an implementation of the BlackScholes model." Mathematica journal, 1, 75~79.