Volume 9, Issue 2

Articles
In and Out
Trott's Corner
New Products
New Publications
Calendar
News Bulletins
New Resources
Classifieds

Editorial Policy
Staff
Submissions
Subscriptions
Back Issues
Contact Information

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

1. Introduction

Parabolic partial differential equations (PDEs), of which the standard heat equation is a special case, are instrumental for many applications in physics, engineering, economics, finance, weather forecasting, and many other fields. In practice, the solution to most parabolic PDEs can be computed only approximately by way of some special numerical procedure. Currently, in most if not all practical implementations, this procedure happens to be some variation of the well-known finite difference method.

This method involves tools and terminology that are familiar to most specialists in numerical methods. Its underlying principle is relatively easy to understand: it comes down to solving sequentially a very large number of linear systems. Unfortunately, in most applications the size of each system is overwhelming and the practical implementation of the method is not straightforward, especially in the case of non-constant coefficients.

It has been well known for nearly as long as the finite difference method has been in existence that the solution of the heat equation can be expressed as a path integral (also known as a Feynman integral). For a mathematician, this is nothing but an expected value relative to Wiener's probability distribution on the space of sample paths. This construction also admits a generalization: the solution of a fairly general parabolic PDE can be expressed as an expected value (alias: a Feynman integral) relative to a special probability distribution on the pathspace (alias: the Wiener measure). Although possible in principle, the numerical calculation of path integrals has always been seen as a formidable task and a rather inefficient way to produce a solution to a parabolic PDE. As the examples included in this article demonstrate, not only is this no longer the case, but, with the advent of Mathematica 4, the solution to a rather general parabolic PDE can be produced by merely expressing the associated path integral in Mathematica.

The method developed in this paper is a good illustration of the way in which significant new developments in the common man's computing tools transform the type of mathematics that is commonly perceived as interesting and important. We will show with concrete examples (see Section 4 below) that computing the price of an European call option--a classical problem in computational finance--in terms of path integrals is just as easy, just as accurate, and virtually just as fast as calculating the price directly by using the Black-Scholes formula. In other words, for all practical purposes, the path integral representation of the solution of the Black-Scholes PDE in Mathematica is just as good as a closed form solution. Unlike the true closed form solution, however, the path integral representation can be used with any reasonably behaved payoff. Of course, this approach would have been unthinkable when Black and Scholes developed their method [1]. Their calculation was interesting and profoundly important, although the really important discovery that they made was the principle that led to the Black-Scholes PDE--not the Black-Scholes formula, as is commonly believed.

It is remarkable how Mathematica's ability to combine advanced symbolic computing with advanced numerical procedures can turn a fairly abstract mathematical concept (in this case, integration on the pathspace) into a rather practical computing tool. This brings a method to consumers of computing technology that is both universal and "smart," in the sense that it is capable of taking advantage of the structure of the PDE. For example, the method will need fewer iterations if the coefficients are close to linear, and can work on a very coarse grid in the state space if the solution is reasonably close to a polynomial. Above all, its implementation is extremely straightforward, involves very little programming, and does not require any knowledge of sophisticated numerical methods, which Mathematica uses implicitly. For a financial engineer, this means that pricing an option with any reasonably-behaved payoff function on an asset that follows a rather general time-homogeneous Markov process, including options that allow early exercise, would be just as easy as pricing a vanilla European call or put option on an asset that follows a standard geometric Brownian motion process. Such technology can have a considerable impact on the financial industry. Indeed, it has long been recognized that actual stock prices exhibit self-regulatory patterns that are inconsistent with the properties of the geometric Brownian motion. Yet the Black-Scholes model remains immensely popular mainly because, with the current level of widely available computing technology, any other model would be virtually intractable for an average practitioner.

Finding the most efficient way to present the computational aspects of the path integral method and its implementation in Mathematica is a bit of a challenge. Indeed, on the one hand, specialists in (and consumers of) numerical methods for PDEs do not usually speak the language of path integrals. On the other hand, specialists in mathematical physics and probability, who do speak the language of path integrals, are not usually interested in numerical methods and computer programming. This article grew out of the desire to meet this challenge. In spite of its use of high-end mathematical jargon, it is mainly about Mathematica and the impact that Mathematica can make not only on the way in which important computational problems are being solved, but also on our understanding of which mathematical methods are "applied" and which are "theoretical."

Note about the notation: Following the common convention in Mathematica, in all equations and algebraic expressions, arguments of functions are wrapped with square brackets and parenthesis are used only for grouping. Double-struck square brackets are used to denote intervals: open intervals are denoted as and closed intervals are denoted as .

For the purpose of illustration, all time-consuming operations in this notebook are followed by //Timing. The CPU times reflect the speed of a 2.8Ghz 4 processor with 768MB of RAM running Mathematica 5.0.