The Mathematica Journal
Volume 9, Issue 2


In This Issue
Tricks of the Trade
In and Out
Trott's Corner
New Products
New Publications
News Bulletins
New Resources

Download This Issue 

About the Journal
Editorial Policy
Back Issues
Contact Information

Path Integral Methods for Parabolic Partial Differential Equations with Examples from Computational Finance
Andrew Lyasoff

3. Path Integrals as Approximation Procedures

The multiple integral in (2.5) is actually an iterated integral, which can be transcribed as a recursive rule. More specifically, if for some fixed , the function satisfies the equation

for and for , assuming that the map is given, one can compute the maps consecutively for , according to the following procedure

As we will soon see, the implementation of this procedure in Mathematica is completely straightforward. Since the above identity holds for every choice of , this raises the question: why does one need the recursive procedure at all? After all, one can always take and then use the recursive rule in (3.1) only for , which comes down to what is known as the Feynman-Kac formula:

Unfortunately, the integral kernel can only be expressed in a form that allows one to actually compute this last integral--for example, by way of some numerical procedure--for some special choice of the coefficients and . The idea, then, is to use the recursive rule (3.1) with some "sufficiently small" and hope that for such a choice of , can be replaced with some approximate kernel, which, on the one hand, makes the integral (the expected value) on the right side of (3.1) computable and, on the other hand, does not allow the error in the calculation of caused by this approximation to grow larger than . We will do what appears to be completely obvious and intuitive: replace and with their first-order Taylor expansions about the point . Intuitively, it should be clear that when the process , is governed by the stochastic equation in (2.6) over a short period of time, its behavior cannot be very different from the behavior of a process that starts from the same position and is governed by the first-order Taylor approximation of the stochastic equation in (2.6) around the point . The reason why one has to settle for a first-order approximation is that, in general, only when the coefficients and in (2.6) are affine functions can one solve this equation explicitly, that is, write the solution as an explicit function of the driving Brownian motion. In general, the probability density of this explicit function is not something that is known to most computing systems and therefore its use in the computation of the integral on the right side of (3.1) is not straightforward. One way around this problem is to compute, once and for all, the transition probability density of a Markov process governed by a stochastic differential equation (SDE) with affine coefficients, and then add that function of seven variables (four parameters determined by the two affine coefficients and plus the variables , , and ) to the list of standard functions in the system. In this exposition we will not go that far, however. Instead, we will develop a second level of approximation, in which the solution of a general SDE with affine coefficients will be replaced by a computable function of the starting point , the time , and the position of the driving Brownian motion at time . After this second level of approximation, the computation of the integral on the right side of (3.1) comes down to computing a Gaussian integral, which, as we will soon see, is not only feasible, but is also extremely easy to implement. Note that when the coefficients and are linear functions, this second level of approximation of the solution of the stochastic equation in (2.6) actually coincides with the exact solution and is therefore redundant.

Before we can implement the plan that we just outlined, we must address an entirely mundane question: how is each function , , going to be represented in the computing system? The answer is that we will use polynomial interpolation; that is, in every iteration, the integral on the right side of (3.1) (with the approximate integral kernel) will be computed only for

where , , and are given parameters. Furthermore, in the actual computation of the integral, the map will be replaced by the map constructed by way of polynomial interpolation (and extrapolation) from the values

In other words, our strategy is to store each of the functions , for , , as a list of real values and interpret as the result of polynomial interpolation from that list. This means that we will be constructing the solution only for and , , for some appropriate choice of . In most practical situations, one needs to construct only on some finite interval and only for a fixed . To do this with our method, one must set for some sufficiently large , and then iterate (3.1) times after choosing and so that the interval of interest is covered by the interval ; the number of interpolation nodes () must be chosen so that the polynomial interpolation is reasonably accurate. This might appear as a daunting programming task, but, as we will soon see, in Mathematica it involves hardly any programming at all. The reason is that the integration routine can accept as an input the output from the interpolation routine--all in symbolic form. This illustrates perfectly the point made earlier: new developments in the common man's computing tools can change the type of mathematics that one can use in practice.

A good point can and should be made that approaching (from a purely computational point of view) a general parabolic PDE like the one in (2.3) in the way we did and with the type of mathematics that we have chosen, namely, integration on the path space, is more than just natural. Indeed, the connection between the distribution of the sample paths of the process , , governed by the stochastic equation in (2.6) and the parabolic PDE (2.3) runs very deep (see [2, 3] for a detailed study of this relationship). The point is that the distribution on the space of sample paths associated with the SDE (2.6) encodes all the information about the elliptic differential operator from (2.4) and, therefore, all the information about the function

This is analogous to the claim that the integral curves of a vector field encode all the information about the vector field itself. Actually, this is much more than a mere analogy. To see this, notice that if were a first-order differential operator (i.e., a vector field) in the open interval , if were an integral curve for starting from , and if has the property , then one would have

In other words, the integral curves of the vector field can be characterized as "curves along which any solution of remains constant." This is interesting as well as quite practical. Indeed, if must be solved for and with boundary condition , , one can compute the solution by first solving the ODE with initial data . Here we treat as a "true" vector field (i.e., as a map of the form ) and then by evaluating the function (assuming that it is given) at the point . Thus, the computation of , for a given and , essentially comes down to constructing the integral curve for the vector field with initial position .

Of course, the remark just made can be found in any primer on PDEs and demonstrates how natural our approach to parabolic PDEs actually is. Indeed, the Feynman-Kac formula in (3.2) is essentially an analog of the relation in (3.3) in the case where is a second-order differential operator of elliptic type. This amounts to the claim that the solution of the parabolic PDE remains constant on average along the sample path of the process , , which is exactly the sample path representation of the semigroup , . Clearly, for any given starting point , the probability distribution of the sample path of the process , , has the meaning of a "distributed integral curve" for the second-order differential operator . This interpretation of the Feynman-Kac formalism is due to Stroock [3] and has been known for quite some time. The only "novelty" here is our attempt to utilize its computing power. The recursive rule (3.1) is analogous to a numerical procedure for constructing an integral curve of a standard (first-order) vector field by solving numerically a first-order ordinary differential equation.

One should point out that the main concept on which all finite difference methods rest is quite natural, too. Indeed, the intrinsic meaning of a partial derivative is nothing but "a difference quotient of infinitesimally small increments." However, such a concept is intrinsic in general and not necessarily intrinsic for the equation . This is very similar to encoding an image of a circle as a list of a few thousand black or white pixels. There are various methods that would allow one to compress such an image, but they all miss the main point: a circle is determined only by its radius and center point, and this is the only data that one actually needs to encode. It is not surprising that, in order to achieve a reasonable accuracy, essentially all reincarnations of the finite difference method require a very dense grid on the time scale as well as the state-space.

As was pointed out earlier, path integrals can be used just as easily in the case of free boundary value problems of the following type: given an interval and a number , find a function and also a free time-boundary so that

This is an optimal stopping problem: at every moment , one observes the value of the system process , , which is governed by the stochastic equation in (2.6), and decides whether to collect the termination payoff or take no action and wait for a better opportunity. The function , which satisfies the conditions in (3.4), gives the price at time () of the right to collect the payoff at some moment during the time period , assuming that at time the observed value of the system process is . Also, at time (), the set gives the range of values of the system process for which it is optimal to collect the payoff immediately, while the set gives the range of values of the system process for which it is optimal to wait for a better opportunity. The valuation and the optimal exercise of an American stock option with payoff , being the price of the underlying asset at time , is a classical example of an optimal stopping problem. Of course, here we assume that the exercise region has the form . This assumption is always satisfied if is a decreasing function of the system process, provided that larger payoffs are preferable. As it turns out, the recursive procedure (3.1) can produce, in the limit, the solution of the above optimal stopping problem with the following trivial modification:

This is simply a rephrasing of Bellman's principle for optimality in discrete time, and the conditions in (3.4) are simply a special case of Bellman's equation.

We conclude this section with the remark that the above considerations can easily be adjusted to the case where the elliptic differential operator contains a 0-order term, that is,

where is a given parameter. The reason why such a generalization is actually trivial is that when satisfies the equation with , then the function

satisfies the same equation, but with . In any case, if is chosen as in (3.6) and not as in (2.4), the only change that one will have to make in the procedure outlined earlier in this section would be to replace the recursive procedure in (3.1) and (3.5), respectively, with


About Mathematica | Download Mathematica Player
Copyright © Wolfram Media, Inc. All rights reserved.