6. 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:

where stands for the differential operator . The method that we developed in Sections 3 and 5 can be implemented in Mathematica with the following three definitions (in addition to the definition of from Section 5), in which the interpolation region is the interval , the number of (equally spaced) interpolation nodes is , and the time step is . The object represents the list of all values of at the entire set of interpolation nodes and, for a given object , which Mathematica treats as a function, returns the list of values of for all .

Note that we have replaced with . 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 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 gives the value of an European call option with a strike price of $80 that expires in years, assuming that the volatility of the stock is 20 percent/year and the interest rate is 10 percent/year. We will take years to maturity, time periods (5 iterations), and 61 interpolation nodes () in the interpolation interval .

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.

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

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

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

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.

The only reason the last two results are different is that applies the Feynman-Kac formula to the interpolated payoff, while 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 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 and the terminal data , the solution of the PDE in (6.1) 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 when the stock price is away from the strike price.

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.

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

Notice that we have chosen the interpolation interval 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 of the equation

with boundary condition .

The interpolation interval is , the number of interpolation nodes is , and the time step is .

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

Sometimes one must use the dynamic programming procedure in a nontrivial way and work with some small time step, even if the coefficients and 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 ; 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 with 41 interpolation nodes and use 20 iterations, which amounts to a time step of 0.15.

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

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 , , starts from a point that is close to the end of the interpolation interval, a sizeable part of the distribution of 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 goes, the error from the extrapolation can be significant. In this particular example, one possible remedy would be to set to for all sufficiently large values of , 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 , that is, as the time step converges to 0, the solution produced after iterations will converge to the actual solution in the entire interpolation interval. This is because when the process , , starts from one of the end-points of that interval, for a very small value of , the distribution of 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).

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.

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.

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.

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).

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.