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

Andrew Lyasoff

By using specific examples—mostly from computational finance—this paper shows that Mathematica is capable of transforming Feynman-type integrals on the pathspace into powerful numerical procedures for solving general partial differential equations (PDEs) of parabolic type, with boundary conditions prescribed at some fixed or free time boundary. Compared to the widely used finite difference methods, such procedures are more universal, are considerably easier to implement, and do not require any manipulation of the PDE. An economist or a financial engineer who works with very general models for asset values, including the so-called stochastic volatility models, will find that this method is just as easy to work with as the classical Black–Scholes model, in which asset values can follow only a geometric Brownian motion process.

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  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 []. 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 \(\(\(]]\) \(a\)\), \(\(b\)\([[\)\)\) and closed intervals are denoted as [[a, 
      b]].

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 \(TraditionalForm\`Intel\^®\ Pentium\^®\)4 processor with 768MB of RAM running Mathematica 5.0.

. Some Background: Solutions to Second-order Parabolic Partial Differential Equations as Path Integrals

This section outlines some well-known principles and relations on which the method developed in this paper is based. A complete, detailed, and rigorous description of these principles can be found in the ground-breaking work of Stroock and Varadhan [] and also in [, ]. The computational aspects of the methods developed in these works have not been explored fully, and to a large extent the present work is a modest attempt to fill this gap.

The Feynman integral of a real function f[] defined on the space (Subscript[𝒞, 0])[[t, 
        T]], which consists of continuous sample paths \(TraditionalForm\`ω : \(\([[\)\(t, 
          T\)\(]]\)\) ↦ ℝ\) with the property (Subscript[ω, t]) = 0, can be expressed as the following formal expression

\(TraditionalForm\`\(\(∫\_{ω  ∈ \ 𝒞\_0[[t, T]]\ }f[ω]\ e\^\(\(-\(1\/2\)\) \(∫\_t\%T\(\((dω\_s\/ds)\)\^2\) ds\)\)\ \((\ ∏\_\(s ∈ \(\([[\)\(t, T\)\(]]\)\)\)dω\_s\/\@\(2  π\)\ )\)\)\(\ \)\(,\)\)\)

which is understood as the limit of the following sequence of multiple integrals

\(TraditionalForm\`\(\(\(\((2  π)\)\^\(-\(k\/2\)\)\) \(∫\_ℝ … \(∫\_\(\(ℝ\)\(\ \)\)f\_k[ω\_0, ω\_d, …, ω\_\(k d\)]\ e\^\(-\(∑\_\(i = 1\)\%k\((ω\_\(i d\) - ω\_\(\((i - 1)\) d\))\)\^2\/\(2  d\)\)\)×d\((ω\_d - ω\_0)\) … d\((ω\_\(k d\) - ω\_\(\((k - 1)\) d\))\)\)\)\)\(,\)\(\ \)\)\)

where k⩾1, d := ((T - t)/k) and \(TraditionalForm\`f\_k[ω\_0, ω\_d, …, ω\_\(k d\)]\) is the result of the evaluation of f[] at the piecewise linear sample path from the space (Subscript[𝒞, 0])[[t, 
        T]], whose values at the points t + i d, 0 ⩽ i ⩽ k, are given by the vector \(TraditionalForm\`\((ω\_0, ω\_d, …, ω\_\(k d\))\)\). It would not be an exaggeration to say that the convergence of the sequence in (.)—and therefore the fact that the integral in (.) is meaningful—is a miracle, without which very little of what we currently know about the universe would have taken place. Indeed, “miracle” is the only term that can explain the fact that, as was pointed out in [], the infinities that arise under the exponent in (.) not only offset the infinities that arise under the product sign, but also offset them in such a way that the result is a perfectly well defined and highly nontrivial probability measure on the space (Subscript[𝒞, 0])[[t, 
        T]], known as the Wiener measure. Not only is the integral in (.) the perfectly meaningful “expected value of f[],” but when f[] is given by f[ω] := Λ[
          x + (Subscript[ω, T])], \(TraditionalForm\`ω ∈ 𝒞\_0[[t, T]]\), for some fixed x ∈ ℝ and some reasonably behaved function \(TraditionalForm\`Λ : ℝ ↦ ℝ\), then, as a function of x and t, this expected value coincides with the solution of the standard heat equation

\(TraditionalForm\`∂\_t u[t, x] + ℒ\ u[t, x] = 
      0\ ,  ℒ := \(1\/2\) ∂\_\(x, x\) , \ \ \ \ x ∈ ℝ\ ,      t ⩽ 
      T\ , \)

with boundary condition u[T, x] ≡ Λ[x]. In fact, this statement admits a generalization: if ℒ is replaced by a general elliptic differential operator of the form

ℒ = (1/2) (σ^2)[
          x] (∂_(x, x))+a[
              x] (∂_x),

assuming that the coefficients σ[] and a[] are reasonably behaved, the solution of the boundary value problem (.) can still be expressed as a path integral similar to the one in (.). The only difference is that in this case the path integral is understood as the limit of the sequence

\(TraditionalForm\`\(\(\ \ \ \ \ \)\(∫\_ℝd\((ω\_d - x)\) 
        p\_d[x, ω\_d]×\ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \(∫\_ℝd\((ω\_\(2  d\) - ω\_d)\) 
            p\_d[ω\_d, ω\_\(2  d\)]\ ×⋯×\(∫\_ℝd\((ω\_\(k d\) - ω\_\(\((k - 1)\) d\))\) 
                p\_d[ω\_\(\((k - 1)\) d\), ω\_\(k d\)] Λ[
                  x + ω\_\(k d\)]\)\), 
													k ≥ 
      1\ , \)\)\)

where (Subscript[p, t])[x, y] is the transition probability density of the Markov process (Subscript[X, t]), t ⩾ 0, determined by the stochastic equation

d(Subscript[X, t]) = σ[
            (Subscript[X, t])] d(Subscript[ß, t]) + 
        a[(Subscript[X, t])] dt,

in which (Subscript[ß, t]), t ⩾ 0, is a given Brownian motion process. In other words, (t, x, y)⟶ (Subscript[p, t])[x, y] is characterized by the fact that for any test function φ[] one has

Ε[φ[(Subscript[X, t + h])] |  (Subscript[X, t]) = 
          x ] = (∫_ℝ) φ[y] 
          (Subscript[p, h])[x, y] dy.

We will say that the process (Subscript[X, t]), t ⩾ 0, is a sample path realization of the Markovian semigroup \(TraditionalForm\`e\^\(t ℒ\)\), t ⩾ 0, ℒ being the differential operator from (.), and will refer to the transition density (t, x, y) ⟶(Subscript[p, t])[x, y] as the integral kernel of the semigroup \(TraditionalForm\`e\^\(t ℒ\)\), t ⩾ 0. This integral kernel is also known as the fundamental solution of the equation ((∂_t)-ℒ) u[t, x] = 
        0, that is, a solution which satisfies the initial condition u[0, x] = (Subscript[δ, y])[x], where  (Subscript[δ, y])[] stands for the usual Dirac delta function with mass concentrated at the point y ∈ 𝒟. This means that, formally, (Subscript[p, t])[x, y] can be expressed as \(TraditionalForm\`p\_t[x, 
          y] = \((\(e\^\(t ℒ\)\) δ\_y)\)[x]\).

A thorough exploration of this generalization of the Feynman integral can be found in the seminal work of Stroock and Varadhan []. The terminology is somewhat different from the one used here, that is, instead of “path integrals,” it comes down to a study of probability distributions on the pathspace associated with equations similar to the one in (.). Such a framework is much more rigorous and powerful. However, our goal is to develop a numerical procedure that computes the solution of the PDE in (.) as an expected value, and in this respect Feynman’s formalism proves to be quite useful.

. Path Integrals as Approximation Procedures

The multiple integral in (.) is actually an iterated integral, which can be transcribed as a recursive rule. More specifically, if for some fixed T ⩾ 0, the function (t, x) ⟶ u[t, x] satisfies the equation

\(TraditionalForm\`\((∂\_t\(+ℒ\))\) u[t, x] = 
      0\ , \ \ \ \ \  ℒ ≡ \(1\/2\) \(σ\^2\)[
          x] ∂\_\(x, x\)\(+a[
              x]\) ∂\_x\ , \)

for t ⩽ T and for x ∈ ℝ, assuming that the map x ⟶ u[T, x] is given, one can compute the maps x ⟶ 
        u[T - i d, x] consecutively for i = 1, 2, …, according to the following procedure

\(TraditionalForm\`u[T - i d, 
        x] = \(∫\_ℝ p\_d[x, y]×
            u[T - \((i - 1)\) d, 
              y] dy = Ε[ u[
                T - \((i - 1)\) d, X\_d]\  | \ X\_0 = 
            x\ ]\)\ , 
										i = 1, 
    2, \(\(…\)\(.\)\)\)

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

\(TraditionalForm\`u[t, 
        x] = ∫\_ℝ p\_\(T - t\)[x, y]×
            u[T, y] dy ≡ Ε[ u[T, X\_\(T - t\)]\  | \ X\_0 = 
            x\ ] ≡ \(\(Ε[ Λ[X\_\(T - t\)]\  | \ X\_0 = 
              x\ ]\)\(.\)\)\)

Unfortunately, the integral kernel \(TraditionalForm\`\((t, x, y)\)\ ⟶
          p\_\(T - t\)[x, 
            y] ≡ \((\(e\^\(\((T - 
                        t)\) ℒ\)\) δ\_y)\) \([x]\)\) 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 a[]. The idea, then, is to use the recursive rule (.) with some “sufficiently small” d > 0 and hope that for such a choice of d, (Subscript[p, d])[x, y] can be replaced with some approximate kernel, which, on the one hand, makes the integral (the expected value) on the right side of (.) computable and, on the other hand, does not allow the error in the calculation of u[T - i d, x] caused by this approximation to grow larger than i d o(d). We will do what appears to be completely obvious and intuitive: replace σ[] and a[] with their first-order Taylor expansions about the point x ∈ ℝ. Intuitively, it should be clear that when the process (Subscript[X, t]), t ⩾ 0 is governed by the stochastic equation in (.) over a short period of time, its behavior cannot be very different from the behavior of a process that starts from the same position (Subscript[X, 0]) = x and is governed by the first-order Taylor approximation of the stochastic equation in (.) around the point x ∈ ℝ. The reason why one has to settle for a first-order approximation is that, in general, only when the coefficients σ[] and a[] in (.) 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 (.) 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 a[] plus the variables t, x, and y) 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 x, the time t, and the position of the driving Brownian motion at time t. After this second level of approximation, the computation of the integral on the right side of (.) 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 a[] are linear functions, this second level of approximation of the solution of the stochastic equation in (.) 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 y ⟶u[T - i d, y], i ⩾ 0, 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 (.) (with the approximate integral kernel) will be computed only for

\(TraditionalForm\`\(\(x ∈ 𝒫 := {Σ + j  h\/m; \ \(-m\) ⩽ j ⩽ 
          m}\)\(,\)\)\)

where Σ, h, and m are given parameters. Furthermore, in the actual computation of the integral, the map y ⟶
        u[T - (i - 1) d, y] will be replaced by the map constructed by way of polynomial interpolation (and extrapolation) from the values

u[
      T - (i - 1) d, Σ + 
        j  (h/m)],         -m ⩽ 
      j ⩽ m.

In other words, our strategy is to store each of the functions u[t, ·], for t = T - i d, i = 1, 2, …, as a list of 2  m + 1 real values and interpret u[t, ·] as the result of polynomial interpolation from that list. This means that we will be constructing the solution u[t, x] only for \(TraditionalForm\`x ∈ \ \(\([[\)\(Σ - h, Σ + h\)\(]]\)\)\) and t = T - i d, i = 1, 2, …, for some appropriate choice of d > 0. In most practical situations, one needs to construct u[t, ·] only on some finite interval and only for a fixed t < T. To do this with our method, one must set d = ((T - t)/k) for some sufficiently large \(TraditionalForm\`k ∈ \(ℤ\_+\)\), and then iterate (.) k times after choosing Σ and h so that the interval of interest is covered by the interval [[Σ - 
        h, Σ + h]]; the number of interpolation nodes (2  m + 1) 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 (.) 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 (Subscript[X, t]), t ⩾ 0, governed by the stochastic equation in (.) and the parabolic PDE (.) runs very deep (see [, ] for a detailed study of this relationship). The point is that the distribution on the space of sample paths associated with the SDE (.) encodes all the information about the elliptic differential operator ℒ from (.) and, therefore, all the information about the function

\(TraditionalForm\`u[t, 
        x] = \((e\^\(\((T - t)\) ℒ\))\) Λ[
          x]\ , \ \ \ \ \ x ∈ ℝ\ ,    t ⩽ \(\(T\)\(.\)\)\)

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 \(TraditionalForm\`𝒟 ⊂ ℝ\), if t ⟶ (Subscript[X, t]) were an integral curve for ℒ starting from (Subscript[X, 0]) = x, and if (t, x) ⟶ u[t, x] has the property ((∂_t)+ℒ) u[t, x] = 
        0, then one would have

u[T, (Subscript[X, T - t])] - 
        u[t, x] = (∫_t^T)(∂_s) u[s, 
                (Subscript[X, s - t])] ds = (∫_t^T) 
              0 ds = 0.

In other words, the integral curves of the vector field ℒ can be characterized as “curves along which any solution of ((∂_t)+ℒ) u[t, x] = 
        0 remains constant.” This is interesting as well as quite practical. Indeed, if ((∂_t)+ℒ) u[t, x] = 
        0 must be solved for t ⩽ T and x ∈ 𝒟 with boundary condition u[T, x] = Λ[x], x ∈ 𝒟, one can compute the solution u[t, x] by first solving the ODE (d/d)t 
          (Subscript[X, t]) = ℒ[(Subscript[X, t])] with initial data (Subscript[X, 0]) = x. Here we treat ℒ as a “true” vector field (i.e., as a map of the form \(TraditionalForm\`ℒ : 𝒟 ↦ ℝ\)) and then by evaluating the function Λ[] (assuming that it is given) at the point (Subscript[X, T - t]). Thus, the computation of u[t, x], for a given x ∈ 𝒟 and t < T, essentially comes down to constructing the integral curve [[0, 
        T - t]] ∋ 
        s ⟶ (Subscript[X, s]) for the vector field ℒ with initial position (Subscript[X, 0]) = x.

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 (.) is essentially an analog of the relation in (.) 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 ((∂_t)+ℒ) u[t, x] = 
        0 remains constant on average along the sample path of the process (Subscript[X, t]), t ⩾ 0, which is exactly the sample path representation of the semigroup \(TraditionalForm\`e\^\(t ℒ\)\), t ⩾ 0. Clearly, for any given starting point (Subscript[X, 0]) = x, the probability distribution of the sample path of the process (Subscript[X, t]), t ⩾ 0, 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 [] and has been known for quite some time. The only “novelty” here is our attempt to utilize its computing power. The recursive rule (.) 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 ((∂_t)+ℒ) u[t, x] = 
        0. 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 \(TraditionalForm\`𝒟 ⊂ ℝ\) and a number \(TraditionalForm\`T ∈ \(ℝ\_+\)\), find a function \(TraditionalForm\`\(\(\(\(]]\)\(-\)\(∞\)\), T\)\( \)\(]]\)\)×𝒟 ∋ \((t, x)\)\) ⟶ u[t, 
            x] ∈ ℝ and also a free time-boundary \(TraditionalForm\`\(\(\(\(]]\)\(-\)\(∞\)\), T\)\( \)\(]]\)\) ∋ t\) \(TraditionalForm\`\(\(⟶\)\(\ \)\(\(x\^*\)[
            t]\)\) ∈ 𝒟\) so that

({{1), u[t, x]⩾Λ[x],for x ∈ 𝒟 ;}, {2), u[t, x] = Λ[x],for\(\(\(\(x ∈ \(\(𝒟\)\(⋂\)\)\)\(\ \)\(]]\)\) - ∞, \(x\^*\)[t]\)\(]]\)\) ;}, {3), ((∂_t)+ℒ) u[t, x] = 
                0,for\(\(\(x ∈ \(\(𝒟\)\(⋂\)\)\)\(\ \)\(]]\)\) \(x\^*\)[
                        t], \(\(∞\)\([[\)\)\);}, {4), \(lim\_\(x\  ↘ \(x\^*\)[t]\)u[t, x]\)=\(Λ[\(x\^*\)[t]]\) and \(\((∂\_x u)\)[
                        t, \(x\^*\)[
                          t]] = \(Λ\^′\)[\(x\^*\)[t]]\);}, {5), u[T, x]=Λ[x]forx ∈ 𝒟.}})

This is an optimal stopping problem: at every moment t ⩽ T, one observes the value of the system process (Subscript[X, t]), t ⩾ 0, which is governed by the stochastic equation in (.), and decides whether to collect the termination payoff Λ[(Subscript[X, t])] or take no action and wait for a better opportunity. The function (t, x) ⟶ u[t, x], which satisfies the conditions in (.), gives the price at time t (t < T) of the right to collect the payoff at some moment during the time period [[t, 
      T]], assuming that at time t the observed value of the system process is (Subscript[X, t]) = x. Also, at time t (t < T), the set \(\(\(\(]]\)\(-\)\(∞\)\), \(x\^*\)[t]\)\(]]\)\)⋂𝒟 gives the range of values of the system process for which it is optimal to collect the payoff immediately, while the set \(\(]]\) \(\(x\^*\)[t]\)\),∞[[ ⋂𝒟 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 Λ[(Subscript[X, t])], (Subscript[X, t]) being the price of the underlying asset at time t, is a classical example of an optimal stopping problem. Of course, here we assume that the exercise region has the form \(\(\(\(]]\)\(-\)\(∞\)\), \(x\^*\)[t]\)\(]]\)\). 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 (.) can produce, in the limit, the solution ]]-∞, 
                T]]×𝒟∋(t, x)  ⟶ u[t, x] of the above optimal stopping problem with the following trivial modification:

u[T - i d, x] = 
      Max[Λ[
          x], (∫_ℝ) (Subscript[p, d])[x, y]×
            u[T - (i - 1) d, 
              y] dy],       i = 1, 
    2, ….

This is simply a rephrasing of Bellman’s principle for optimality in discrete time, and the conditions in (.) 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,

ℒ = (1/2) (σ^2)[
          x] (∂_(x, x))+a[
              x] (∂_x)-r

where r ∈ ℝ is a given parameter. The reason why such a generalization is actually trivial is that when (t, x) ⟶ u[t, x] satisfies the equation ((∂_t)+ℒ) 
        u[t, x]=0 with ℒ = (1/2) (σ^2)[
            x] (∂_(x, x))+a[x] (∂_x), then the function

(t, 
          x) ⟶ (e^(-(T - t) r)) 
      u[t, x]

satisfies the same equation, but with ℒ = (1/2) (σ^2)[
            x] (∂_(x, x))+a[x] (∂_x)-r. In any case, if ℒ is chosen as in (.) and not as in (.), 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 (.) and (.), respectively, with

\(TraditionalForm\`u[T - i d, 
        x] = ∫\_ℝ\( e\^\(\(-r\) d\)\) p\_d[x, y]×
          u[T - \((i - 1)\) d, y] dy, \ \ \  i = 
      1, 2, …, \)

and

\(TraditionalForm\`u[T - i d, x] = 
      Max[Λ[
          x], ∫\_ℝ\( e\^\(\(-r\) d\)\) p\_d[x, y]×
            u[T - \((i - 1)\) d, y] dy]\ , \ \ i = 1, 
    2, \(\(…\)\(.\)\)\)

. The Black–Scholes Formula versus the Feynman–Kac Representation of the Solution of the Black–Scholes PDE

In the option pricing model developed by F. Black and M. Scholes [] the asset price S is assumed to be governed by the stochastic equation

\(TraditionalForm\`dS\_t = ς\ \(S\_t\) dß\_t + μ\ \(S\_t\) dt\ , \ \  t ⩾ 0\ , \)

in which ς > 0 and μ > 0 are given parameters. At time t > 0 the price of an European call option with payoff Λ[S] = Max[S - K, 0], which expires at time T > t, can be identified with the expression u[t, (Subscript[S, t])], where (t, x) ⟶ u[t, x] is the solution of the PDE:

(∂_t) u[t, 
                x] + (1/2) (σ^2) (x^2) ((∂_x)^2) 
                  u [t, x] + 
            r x (∂_x) u [t, x] - 
            r u[t, x] = 0,   t ⩽ T,  x > 0,

with boundary condition (lim_(t ↗ T))u[t, 
            x] = Λ[x], where r > 0 is a parameter representing the (constant) short-term interest rate and Λ[] gives the termination payoff as a function of the price of the underlying asset. The fact that the option price must satisfy the PDE (.) may be seen as a consequence of Bellman’s principle for optimality, but one has to keep in mind that the use of this principle in finance is somewhat subtle. Notice, for example, that the Feynman–Kac formula (.) links the PDE (.) not with the stochastic equation in (.), but, rather, with the stochastic equation

\(TraditionalForm\`dX\_t = ς\ \(X\_t\) dß\_t + 
        r\ \(X\_t\) dt\ , \ \ \ t ⩾ 
      0. \)

Therefore, the solution of the Black-Scholes equation in (.) can be expressed as an expected value of the form

\(TraditionalForm\`u[t, 
        x] = Ε[\(e\^\(-\(r (
                      T - t)\)\)\) Max[X\_\(T - t\) - K, 0]\  | \ X\_0 = 
          x], \ \ t < 
      T,   x ∈ \(ℝ\_+\), \)

and not as

\(TraditionalForm\`u[t, 
        x] = Ε[\(e\^\(-\(r (
                      T - t)\)\)\) Max[S\_\(T - t\) - K, 0]\  | \ S\_0 = 
          x], \ \ t < 
      T,     x ∈ \(\(ℝ\_+\)\(.\)\)\)

In other words, in the context of option pricing, one must use Bellman’s principle as if the predictable part in the system process has the form r (Subscript[X, t]) dt, regardless of the form that the predictable part in the system process actually has. This phenomenon is well known and its justification, of which there are several, can be found in virtually any primer on financial mathematics. Roughly speaking, this is caused by the fact that one must take into account the reward on average which investors collect from the underlying asset (see []).

Since the stochastic equation in (.) admits an explicit solution given by

\(TraditionalForm\`X\_t = \(X\_0\) e\^\(ς ß\_t\  + \ \((r - \(1\/2\) ς\^2)\) t\), \ \ \ t ⩾ 0, \)

and since this solution is simple enough, it is possible to implement the recursive procedure from (.) with a single iteration. This means that, in this special case, the Feynman–Kac formula in (.) is directly computable without the need to approximate the integral kernel of the associated Markovian semigroup. This also eliminates the need for interpolation, if all that one wants to do is compute the solution u[t, x] for some fixed t and some fixed x. Of course, this assumes that the payoff can be expressed in a way that Mathematica can understand. In our first example, we will compute u[0, 50], assuming that ς = 0.2, r = 0.1 and that the boundary condition is u[3, x]= Λ[x], where Λ[x] = Max[x - 40, 0].

Here is the result.

OptionPrice[50] // Timing

What we have just computed was the price of an European call option with a strike price of $40 and T = 3 years to maturity, assuming that the current price of the underlying asset is $50, that the volatility in the stock price is 20 percent/year and that the fixed interest rate is 10 percent/year. Of course, the price of such an option can be computed directly by using the Black–Scholes formula. First, we will express the Black–Scholes formula in symbolic form

and then compute the price of the same option.

BS[r, ς, 50, 40, T] // Timing

One can see from this example that, from a practical point of view, the Feynman–Kac representation of the solution of the Black–Scholes PDE in Mathematica is just as good as the “explicit” Black–Scholes formula. (Ironically, the Black–Scholes formula takes somewhat longer to type, not to mention the time it might take to retrieve it from a book.) However, the Feynman–Kac representation is much more universal. Indeed, with no extra work—and no additional complications whatsoever—it allows one to compute the value of an option with a very general payoff. Consider for example an option that pays at maturity the smallest of the payoffs from a call with a strike price of $40 and a put with a strike price of $80.

With all other parameters kept the same as before (3 years left to maturity, 20 percent volatility and 10 percent interest rate), assuming that the stock price is $50, the price of this option is

OptionPrice[50] // Timing

Obviously, the computing time is negligible, but what is much more important is that the only thing that we had to do in this case was to simply change the definition of the payoff. Here is another example.

Λ[x_] := 
      Max[50 - 0.05 (((x - 40))^2), 0]; 
    Plot[Λ[x], {x, 0, 100}]

OptionPrice[50] // Timing

We will conclude this section with an example of what is known as a cash-or-nothing option.

Off[NIntegrate::"slwcon"]; 
    OptionPrice[50] // Timing

The overwhelming majority of the options that are currently traded are some combination of calls, puts, and cash-or-nothing options. The origins of this tradition are not exactly certain, but it might have something to do with the fact that these are the options for which popular and easy-to-use computing tools have been widely available. It should be clear from the examples presented in this section that, in the realm of the Black–Scholes model, essentially any reasonably defined option is, in fact, trivial to price, and that having a Black–Scholes formula is no longer necessary. Of course, the well-known binomial model does allow one to compute the value of an option with a very general payoff; however, this popular approach does require some work, and its accuracy and overall efficiency is far from satisfactory. Furthermore, this method is tied exclusively to the assumption that stock prices follow the law of a geometric Brownian motion; the technology that we are developing is not.

One must realize that, in keeping with the method outlined in Section , solving the PDE in (.) is an entirely trivial matter. The point is that the integral kernel of the associated Markovian semigroup happens to be directly computable, in the sense that it can be written in symbolic form that can be understood by Mathematica. However, our main goal is to demonstrate that this same method is just as easy to implement in much more general situations. In terms of applications to computational finance, that would mean that considerably more sophisticated models for asset values—for example, models with stochastic volatility and stochastic return of a rather general form—would be just as easy to use as the basic Black–Scholes model. The main step in this program is to develop an adequate approximation of the sample path representation of the semigroup \(TraditionalForm\`e\^\(t ℒ\)\), t ⩾ 0, for a general elliptic differential operator ℒ (see the next section).

. Approximation of the Fundamental Solution by Way of a Discrete Dynamic Programming Procedure

In general, there are many ways in which the solution of a given SDE can be approximated by a discrete time series (see [] for a thorough exploration and a comprehensive list of references on this topic). Unfortunately, none of the available approximation schemes was developed for a purpose similar to ours, so we will develop one which actually is. Another extremely interesting approximation method, whose computational aspects are yet to be explored, can be found in [].

To begin with, notice that the stochastic equation in (.) can be written in integral form as

\(TraditionalForm\`\(X\_t\% \) = 
      x + ∫\_0\%t σ[\(X\_s\% \)] dß\_s + ∫\_0\%t 
              a[\(X\_s\% \)] ds\ , \ \ \(X\_0\% \) = 
      x ∈ ℝ,   t ⩾ 0. \)

One possible approximation of this equation can be obtained by replacing the coefficients σ[] and a[] by constants, which gives the following SDE

(Subscript[Y, t]^ ) = 
      x + (∫_0^t) σ[
              x] d(Subscript[ß, s]) + (∫_0^t) a[
              x] ds,   (Subscript[Y, 0]^ ) = 
      x ∈ (ℝ^n).

This approximation, which we will call the (0^th)-order approximation, is so crude that the SDE it leads to admits a trivial solution:

(Subscript[Y, t]^ ) = 
          x + σ[x] (Subscript[ß, t]) + a[x] t,  t⩾0.

Intuitively, it should be clear that, when t is sufficiently small, the distribution of the random variable (r.v.) (Subscript[Y, t]) is not very far from that of the r.v. (Subscript[X, t]). If we were to implement the plan outlined in Section  with this level of approximation of the Markovian semigroup associated with the equation ((∂_t)+ℒ) u[t, x] = 
        0, the recursive procedure in (.) would come down to

\(TraditionalForm\`\(\(\ \)\(u[T - i \( T - t\)\/k, 
        x] = \(\((2  π)\)\^\(-\(1\/2\)\)\) \(∫\_ℝ u[T - \((i - 1)\) \(T - t\)\/k, 
              x + σ[x] \(T - t\)\/k\ γ + 
                a[x] \(T - t\)\/k] \(e\^\(-\(γ\^2\/2\)\)\) dγ\), 
											i = 
      1, …, \ \(\(k\)\(.\)\)\)\)\)

Assuming that u[T - (i - 1) ((T - t)/k), ·] is the result of polynomial interpolation, the above integral is certainly computable for every choice of x: in fact, it only has to be computed for x ∈ 𝒫, where 𝒫 stands for the set of interpolation nodes. One can show that as k -> ∞ and as the number of interpolation nodes increases to ∞, u[T - t, ·] converges to the actual solution uniformly in the interpolation region. However, this convergence is known to be rather slow. The (0^th)-order approximation of SDEs is as old as the theory of SDEs itself and can be traced back to the original works of Itô—in fact, this approximation is crucial for his construction of the stochastic integral. One may say that the (0^th)-order approximation is “optimal” if all that one needs to do is show the existence of solutions to SDEs and derive various properties of the associated probability distributions, in which case the speed of the convergence is not an issue. In our case, however, it is and one can do considerably better by using the (1^st)-order Taylor approximation of the coefficients σ[·] and a[·]. Our plan, then, is to replace the stochastic equation in (.) with the following SDE

(Subscript[Z, t]^ )=x+(∫_0^t)(({{}, {}})σ[x]+(σ^′)[
                          x] ((Subscript[Z, t]^( x)) - x) )d(Subscript[ß, s])+(∫_0^t)(({{}, {}}) a[x]+(a^′)[
                          x] ((Subscript[Z, t]^( x)) - x))ds,  t ⩾ 0,

which admits this explicit solution

\(TraditionalForm\`\(\(\ \)\(\(Z\_t\% \) = e\^\(\(σ\^′\)[x] ß\_t + \((\(a\^′\)[x] - \(1\/2\) \(σ\^′\)[x]\^2)\) t\)\ ×\((\((a[x] - \(a\^′\)[x] 
                    x)\) \(∫\_0\%t\( e\^\(\(-\(σ\^′\)[x]\) ß\_s - \((\(a\^′\)[
                                x] - \(1\/2\) \(σ\^′\)[x]\^2)\) 
                          s\)\) ds\)\ 
	 + \ \((σ[x] - \(σ\^′\)[x] 
                    x)\) \(∫\_0\%t\( e\^\(\(-\(σ\^′\)[x]\) ß\_s - \((\(a\^′\)[
                                x] - \(1\/2\) \(σ\^′\)[x]\^2)\) 
                          s\)\) dß\_s\) + 
            x)\)\ 

    			\(-\(σ\^′\)[
          x]\)×\((σ[x] - \(σ\^′\)[x] 
            x)\)\ \(\(t\)\(.\)\)\)\)\)

In spite of the “explicit” form of this solution, the density of the r.v. (Subscript[Z, t]) is still difficult to compute. Therefore, we need to develop yet another level of approximation. In order to approximate the random variable on the right side of (.) with something computable, notice that if the integrand in the first integral is replaced by the constant 1 this will produce an error of order o(t) (a small o of t). The approximation of the second integral is somewhat more involved. A straightforward application of the Itô formula yields

\(TraditionalForm\`\(\(	\)\(∫\_0\%t\( e\^\(\(-\(σ\^′\)[
                    x]\) ß\_s - \((\(a\^′\)[
                      x] - \(1\/2\) \(σ\^′\)[x]\^2)\) 
                s\)\) dß\_s\ 

    		\(\(\(=\)\(\(1\/\(σ\^′\)[
                x]\) \((1 - e\^\(\(-\(σ\^′\)[x]\) ß\_t - \((\(a\^′\)[x] - \(1\/2\) \(σ\^′\)[x]\^2)\) t\))\)\  + \((\(σ\^′\)[
                x] - \(a\^′\)[x]\/\(σ\^′\)[x])\) \(∫\_0\%t\( e\^\(\(-\(σ\^′\)[
                          x]\) ß\_s - \((\(a\^′\)[
                            x] - \(1\/2\) \(σ\^′\)[x]\^2)\) 
                      s\)\) ds\)\)\)\(,\)\(\ \)\)\)\)\)

which shows that

\(TraditionalForm\`∫\_0\%t\( e\^\(\(-\(σ\^′\)[x]\) ß\_s - \((\(a\^′\)[
                        x] - \(1\/2\) \(σ\^′\)[x]\^2)\) 
                  s\)\) dß\_\(\(s\)\(\ \)\) = \(1\/\(σ\^′\)[
                x]\) \((1 - e\^\(\(-\(σ\^′\)[x]\) ß\_t - \((\(a\^′\)[x] - \(1\/2\) \(σ\^′\)[x]\^2)\) t\))\) + \(σ\^′\)[x] 
          t - \(\(a\^′\)[x]\/\(σ\^′\)[x]\) 
          t + \(\(o(t)\)\(.\)\)\)

Thus, (.) can be rewritten as

\(TraditionalForm\`\(Z\_t\% \) = e\^\(\(σ\^′\)[x] ß\_t + \((\(a\^′\)[x] - \(1\/2\) \(σ\^′\)[x]\^2)\) t\)\)
\(TraditionalForm\`	\(\(×\)\(\ \)\((\((a[x] - \(a\^′\)[x] 
                    x)\) t\  + \((σ[x] - \(σ\^′\)[x] 
                    x)\)
	×\((\(1\/\(σ\^′\)[
                          x]\) \((1 - e\^\(\(-\(σ\^′\)[x]\) ß\_t - \((\(a\^′\)[x] - \(1\/2\) \(σ\^′\)[x]\^2)\) t\))\) - \(\(a\^′\)[x]\/\(σ\^′\)[
                          x]\) t)\) + x)\)\)\  + \(\(o(
          t)\)\(.\)\)\)

Of course, this identity holds only if (σ^′)[x] ≠ 0, but this is the only truly interesting case, since (σ^′)[x] = 0 implies that

\(TraditionalForm\`\(\(\(Z\_t\% \) = e\^\(\(a\^′\)[x] t\)×\((\((a[x] - \(a\^′\)[x] x)\) 
              t + σ[
                x] \(∫\_0\%t\( e\^\(\(-\(a\^′\)[x]\) s\)\) dß\_s\) + 
            x)\)\)\(,\)\)\)

where the second integral is a Gaussian r.v. of variance \(TraditionalForm\`\(\(\(1\/\(2\ \(a\^′\)[
                  x]\)\) \((1 - e\^\(\(-2\)\ \(a\^′\)[x]\ t\))\)\)\(=\)\)\) t -  (a^′)[x] (t^2) + O((t^3)). Thus, we will approximate expressions of the form Ε[f[(Subscript[Z, t]^ )]] with

\(TraditionalForm\`\(1\/\@\(2  π\)\) \(∫\_\(-∞\)\%∞ 
            f[ψ[t, 
              x, γ ]] \(e\^\(-\(γ\^2\/2\)\)\) dγ\)\)

where

In particular, the recursive procedure in (.) will be replaced with

\(TraditionalForm\`u[T - i \( T - t\)\/k, 
        x] = \(1\/\@\(2  π\)\) \(∫\_\(-∞\)\%∞ u[
              T - \((i - 1)\) \(T - t\)\/k, ψ[\(T - t\)\/k, 
                x, γ ]] \(e\^\(-\(γ\^2\/2\)\)\) dγ\), 
										i = 1, …, k, \)

in the case of a terminal value problem, and with

\(TraditionalForm\`u[T - i \( T - t\)\/k, x] = 
      Max[Λ[
          x], \(1\/\@\(2  π\)\) \(∫\_\(-∞\)\%∞ u[T - \((i - 1)\) \(T - t\)\/k, ψ[\(T - t\)\/k, 
                  x, γ ]] \(e\^\(-\(γ\^2\/2\)\)\) dγ\)\ ], 
										i = 1, …, 
    k, \)

in the case of a free boundary value problem. In Section 6 we will show that the computation in Mathematica of the last two integrals for a fixed x ∈ 𝒫 is completely straightforward when the object u[T - (i - 1) ((T - t)/k), ·] is the result of polynomial interpolation.

It is not hard to check that in the special case where σ[x] ≡ ς×x and a[x] ≡ r×x one has

\(TraditionalForm\`\(\(ψ[Δ, x, y ] = 
      x \ e\^\(ς y\ \@Δ + \ \((r - \(1\/2\) ς\^2)\) Δ\)\)\(,\)\)\)

which means that, with this choice for σ[] and a[], the solution to the stochastic equation in (.) is simply (Subscript[X, t]) = ψ[t, x, (Subscript[ß, t]) ].

. Dynamic Programming with Mathematica: The Case of American and European Stock Options

First, we will consider several examples of terminal-value problems of the form:

((∂_t)+ℒ) u[t, x] = 
      0 ,  x ∈ ℝ ,  
    t ⩽ T ,  
    u[T, x] = Λ[x]

where ℒ stands for the differential operator (1/2) (σ^2)[
          x] (∂_(x, x))+a[x] (∂_x)-r. The method that we developed in Sections  and  can be implemented in Mathematica with the following three definitions (in addition to the definition of ψ[t, x, γ] from Section ), in which the interpolation region is the interval [[Σ - 
        h, Σ + h]], the number of (equally spaced) interpolation nodes is 2  m + 1, and the time step is ((T - t)/k). The object \(TraditionalForm\`U[
        i, \(\(·\)\(\(,\)\(k\)\(,\)\(\(·\)\(\(,\)\(\(·\)\(\(,\)\(·\)\)\)\)\)\)\)]\) represents the list of all 2  m + 1 values of u[T - i (( T - t)/k), ·] at the entire set of interpolation nodes 𝒫 and, for a given object f, which Mathematica treats as a function, \(TraditionalForm\`A[
        f, \(\(·\)\(\(,\)\(\(·\)\(\(,\)\(\(·\)\(\(,\)\(·\)\)\)\)\)\)\)]\) returns the list of values of Ε[f[ψ[d, x, γ]]] for all x ∈ 𝒫.

Note that we have replaced (∫_(-∞)^∞) with (∫_(-7)^7). This makes the integration procedure more efficient, while the loss in precision is well within the accuracy of the method. One must be aware, however, that the interval [[-7, 
      7]] may be too narrow for some problems.

Since our goal here is not to improve Mathematica’s performance, we will ignore the following error messages.

In our first example, we take

The solution u[0, x] gives the value of an European call option with a strike price of $80 that expires in T years, assuming that the volatility of the stock is 20 percent/year and the interest rate is 10 percent/year. We will take T = 3 years to maturity, k = 5 time periods (5 iterations), and 61 interpolation nodes (m = 30) in the interpolation interval [[80 - 75, 
      80 + 75]].

μ = 
      U[5, 3, 5, 80, 75, 30] // 
        Timing; μ[[1]]

We can compare the result (the red curve) with the “exact” solution (the blue curve) produced by the Black–Scholes formula; the green curve represents the termination payoff.

Plot[{Λ[
          x], ListInterpolation[μ[[2]], {{80 - 75, 80 + 75}}][x], 
        BS[0.1, 0.2, x, 80, 3]}, {x, 0, 350}, PlotPoints -> 200, 
      PlotRange -> All, 
      PlotStyle -> {Hue[0.3], Hue[0. ], Hue[ .7]}]

One can see that the numerical solution (the red curve) is reasonably close to the exact solution in most of the interpolation region. To develop some sense for the accuracy of the method, we will compare the value of the option when the stock price is $50, produced by the dynamic programming procedure

ListInterpolation[μ[[2]], {{80 - 75, 80 + 75}}][50]

with the following value produced by the Black–Scholes formula.

BS[0.1, 0.2, 50, 80, 3]

Of course, one can also compute this value by using the Feynman–Kac formula, just as we did in Section .

r = 0.1; ς = 0.2; 
    OptionPrice[50]

As was pointed out earlier, when the coefficients are linear, ψ[] gives the exact solution and one can get away with a single iteration, that is, with a time step equal to the entire time interval. Essentially this comes down to using the Feynman–Kac formula directly.

μ = 
      U[1, 3, 1, 80, 75, 
        30]; ListInterpolation[μ, {{80 - 75, 80 + 75}}][
      50]

The only reason the last two results are different is that U[] applies the Feynman–Kac formula to the interpolated payoff, while OptionPrice[] applies the Feynman–Kac formula to the actual payoff. Of course, there is no reason to use the dynamic programming procedure when the coefficients are linear and one deals with a fixed boundary value problem; here we did so only for the purpose of illustration.

Now we will consider an example in which the coefficients σ[] and a[] are not linear, and therefore the Feynman–Kac formula cannot be used directly, that is, with a single iteration.

With this choice of the coefficients σ[] and a[] and the terminal data Λ[], the solution of the PDE in (.) will give the price of a call option with a strike price of $40, but with stochastic volatility for the underlying asset. For the purpose of illustration, we choose a volatility coefficient that depends on the price of the stock: it decreases to 0.1 when the price of the stock approaches the strike price of the option and stays constant at 0.2 when the stock price is away from the strike price.

Plot[σ[x]/x, {x, 0.05, 80}, 
      PlotRange ->  All, PlotPoints -> 200]

mk = U[20, 3, 20, 50, 48, 50] // Timing; 
    mk[[1]]

Now we can compare the price of this option with the one produced by the Black–Scholes formula (the blue line), which ignores the fact that the volatility decreases when the stock price approaches the strike price. Just as one would expect, when the stock price is away from the region where the volatility decreases, the price of the option looks the same as the price produced under the assumption of constant volatility.

Plot[{BS[0.1, 0.3, x, 40, 
          3], ListInterpolation[
            mk[[2]], {{50 - 48, 
                50 + 48}}][x]}, {x, 0, 80}, PlotPoints -> 400, 
      PlotStyle -> {Hue[0.7], Hue[0. ]}]

When the stock price is $50, the price of this option is

ListInterpolation[
        mk[[2]], {{50 - 48, 50 + 48}}][
      50]

Notice that we have chosen the interpolation interval [[50 - 48, 
      50 + 48]] so that it is centered around the stock price of interest, namely $50; remember that the accuracy of the method is the highest at the center of the interpolation interval. Notice also that, while it took much longer to compute, the calculation of the option price under stochastic volatility did not require any additional work.

Our next example, which is unrelated to finance, will compute the solution x ⟶ u[0, x] of the equation

(∂_t) u[t, 
            x] + (e^(-(x^2))) (∂_(x, x))u[t, 
              x] + Sin[x] (∂_x) u[t, x] = 0

with boundary condition u[3, x] = Sign[x].

The interpolation interval is [[0 - 5, 
      0 + 5]], the number of interpolation nodes is 2×25 + 1 = 51, and the time step is (3/30) = 0.1.

g0 = U[30, 3. , 30, 0, 5, 25] // Timing; 
    g0[[1]]

Now we can compare the solution that we just found, that is, the function u[0, · ] (pictured as the blue line below) with the terminal data u[3, ·] (pictured as the red line).

Plot[{Λ[
          x], ListInterpolation[
            g0[[2]], {{-5, 5}}][
          x]}, {x, -4.5, 4.5}, PlotStyle -> {Hue[0], Hue[ .7]}, 
      PlotPoints -> 200]

Sometimes one must use the dynamic programming procedure in a nontrivial way and work with some small time step, even if the coefficients σ[] and a[] happen to be linear; such a case is an optimal stopping problem, which leads to boundary conditions prescribed along a free boundary. Now we will consider several examples that involve stock options in which early exercise is possible (the so-called American options). We only need to change the definition of the function U[]; everything else in the procedure will remain the same.

Under the Black–Scholes assumption for the price process, American call options are the same as European call options because early exercise of such options is never optimal. Thus, in our first example we consider an American put option with a strike price of $40.

First, we will compute the price of the option under the assumption that—just as in the Black–Scholes model—the stock price follows a geometric Brownian motion process (of course, the Black–Scholes formula cannot be used in this case as is, for it is only valid for European options).

We will interpolate the solution in the interval \(TraditionalForm\`\(\([[\)\(40 - 38, 
        40 + 38\)\(]]\)\) ≡ \(\([[\)\(2, 78\)\(]]\)\)\) with 41 interpolation nodes and use 20 iterations, which amounts to a time step of 0.15.

m = U[20, 3, 20, 40, 38, 20] // Timing; 
    m[[1]]

Now we can plot the solution (the red curve below) together with the termination payoff function (the blue curve).

Plot[{Λ[
          x], ListInterpolation[
            m[[2]], {{40 - 38, 
                40 + 38}}][x]}, {x, 0, 80}, PlotPoints -> 200, 
      PlotStyle -> {Hue[0.7], Hue[0. ]}]

To those familiar with the basics of American options, it should be clear from this graph that the approximation that we just produced is not very accurate in the last 1/4 of the interpolation interval (it is well known that the actual stock price decreases to 0 monotonically). This is caused by the fact that, when the process (Subscript[X, t]), t ⩾ 0, starts from a point that is close to the end of the interpolation interval, a sizeable part of the distribution of (Subscript[X, t]) might be spread in the region where the procedure is forced to use extrapolation. As a result, depending on how far inside the extrapolation region the essential part of the distribution of (Subscript[X, t]) goes, the error from the extrapolation can be significant. In this particular example, one possible remedy would be to set u[T - (i - 1) ((T - t)/k), x] to 0 for all sufficiently large values of x, regardless of the value produced by the interpolation (or the extrapolation) procedure. This is a common problem in essentially all numerical methods for PDEs, and the common way around it is to work with a grid in the state-space which covers a considerably larger region than the one in which an accurate representation of the solution is being sought. It is to be noted, however, that, as k  -> ∞, that is, as the time step converges to 0, the solution produced after k iterations will converge to the actual solution in the entire interpolation interval. This is because when the process (Subscript[X, t]), t ⩾ 0, starts from one of the end-points of that interval, for a very small value of t, the distribution of (Subscript[X, t]) will be concentrated mostly around the end-point where the extrapolation procedure is still reasonably accurate. Thus, by using the method developed here, one can construct an arbitrarily accurate representation of the solution in the entire interpolation interval by simply choosing a sufficiently large number of interpolation nodes and a sufficiently small time step. In general, this is not true for the finite difference method. In order to illustrate this phenomenon, we will now repeat the above calculation with a time step, which is three times as small (accordingly, the calculation will run about three times longer).

mm = U[60, 3, 60, 40, 38, 20] // Timing; 
    mm[[1]]

Plot[{Λ[
          x], ListInterpolation[
            mm[[2]], {{40 - 38, 
                40 + 38}}][x]}, {x, 0, 80}, PlotPoints -> 200, 
      PlotStyle -> {Hue[0.7], Hue[0. ]}]

Let us remember that here we are also solving an optimal stopping problem, whose solution is given by the point, which separates the “exercise” region, that is, the region where the price of the option coincides with the termination payoff, from the “hold” region, where the price of the option is larger than the termination payoff.

FindRoot[ListInterpolation[
            mm[[2]], {{40 - 38, 
                40 + 38}}][x] == Max[40 - x, 0], {x, 38, 
        40}]

It is interesting that the solution, which was produced earlier and about three times faster by using a three times larger time step, gives the same result.

FindRoot[ListInterpolation[
            m[[2]], {{40 - 38, 
                40 + 38}}][x] == Max[40 - x, 0], {x, 38, 
        40}]

This illustrates perfectly the point that we made earlier: the result produced with the larger time step was already reasonably accurate in the middle of the interpolation interval.

Now we will compute the price of the same put option, but under the assumption that the stock price exhibits stochastic volatility of the type we introduced earlier in conjunction with European-type options.

mx = U[60, 3, 60, 40, 38, 30] // Timing; 
    mx[[1]]

Now we can compare the last result (the blue curve below) with the one that was produced in the case of constant volatility (the red curve).

Plot[{Λ[
          x], ListInterpolation[
            mm[[2]], {{40 - 38, 
                40 + 38}}][
          x], ListInterpolation[
            mx[[2]], {{40 - 38, 
                40 + 38}}][x]}, {x, 30, 70}, PlotPoints -> 200, 
      PlotStyle -> {Hue[0.3], Hue[0. ], Hue[0.7]}]

We can see that, just as one would expect, the price of this option differs from the one produced under the constant volatility assumption only in the region where the volatility decreases. Also, the decrease of the volatility near the exercise price actually decreases the range of stock prices for which immediate exercise is optimal, as the following computation shows.

FindRoot[ListInterpolation[
            mx[[2]], {{40 - 38, 
                40 + 38}}][x] == Max[40 - x, 0], {x, 23, 
        30}]

. Concluding Remarks

Although the examples presented in Section  were borrowed mostly from the realm of computational finance, it is important to realize that neither the general method developed in Sections  and , nor the Mathematica code used in Section  takes advantage of the special nature of the problems associated with financial applications in any way. The same method—and, indeed, the same computer code—can produce an approximate solution to a rather general one-dimensional parabolic PDE with boundary conditions imposed on some free or fixed boundary. Furthermore, this method does not require any manipulation whatsoever of the PDE.

As is well known, the “no free lunch” principle reigns in the world of computing, too, and the simplicity, the generality, and the precision with which we were able to solve some classical option valuation problems in the previous section did not come easily by any means. In terms of this metaphor, not only was the lunch not free, but it was, in fact, a rather expensive one. Indeed, an efficient integration procedure was absolutely crucial for the implementation of the method in Section . As explained in The Mathematica Book (see [] p. 1071) the function Integrate[ ] “uses about 500 pages of Mathematica code and 600 pages of C code.” The role played by the function ListInterpolation[ ] was just as crucial. Indeed, it allowed us to work on an exceptionally coarse grid in the state-space and then “restore” the values at all remaining points by a fairly sophisticated interpolation procedure. More importantly, we were able to feed NIntegrate[ ] with the object produced by ListInterpolation[ ]. There are intrinsic obstacles to such “shortcuts” in essentially all modifications of the finite difference method.

Acknowledgment

This work is dedicated to D. W. Stroock and S. R. S. Varadhan. All numerical procedures used here are a mere transcript of their insights.

References

[]  F. Black and M. Scholes, “The Pricing of Options and Corporate Liabilities,” Journal of Political Economy, 81, 1973 pp. 637–659.

[]  D. W. Stroock and S. R. S. Varadhan, “Multidimensional Diffusion Processes,” Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], 233, New York, Berlin: Springer-Verlag, 1979.

[]  D. W. Stroock, “Lectures on Stochastic Analysis: Diffusion Theory,” London Mathematical Society Student Texts, 6, 1987, Oxford: Cambridge University Press.

[]  D. W. Stroock, Probability Theory, an Analytic View, Oxford: Cambridge University Press, 1993.

[]  A. Lyasoff, “Three Page Primer in Financial Calculus.” (2000) Boston: Boston University. www.bu.edu/mathfn/people/primer1.pdf.

[]  P. E. Kloeden and E. Platen, Numerical Solutions of Stochastic Differential Equations, New York: Springer-Verlag, 1997.

[]  D. Stroock and S. Taniguchi, “Diffusions as Integral Curves, or Stratonovich without Itô,” in Progress in Probability #34, ed. by M. Freidlin, Birkhäuser, 1994, pp. 331–369.

[]  S. Wolfram, The Mathematica Book, 5th ed. Champaign, Oxford: Wolfram Media/Cambridge University Press, 2003.