Comments on Benninga's article  and Shaw's Chapters 18-20.

General Remarks

The classic article by Simon Benninga, Raz Steinmetz and John Stroughair  "Implementing numerical option pricing models" first appeared in "The Mathematica Journal" in 1993. It gave the first implementation of the Cox and Rubinstein binomial and binary pricing model for valuing financial options.  The programs given in the paper are very elegant and rather easy to understand . There is also a short but clear discussion of the theoretical foundations of the binomial model. The article also discusses the Monte-Carlo approach to option pricing but i shall not  discuss it here but postpone it till later.
Chapters 18-20 of Shaw's book are in fact an  extension of the paper of Benninga et al. Here I am going to try to briefly explain the relation between the article and the chapters form Shaw's book and in what respect the book goes beyond the paper. 

The idea of the binomial (and also trinomial  etc. approach to option pricing) is to model the movements of the price of a share or stock  exchange index in a risk-neutral world by a discrete random walk rather than by a continuous random variable with lognormal distribution, as is the case with the Black-Scholes model. In the binary model the stock, whose value at time t is (Subscript[S, t]) is assumed to have a probability p of going up to u (Subscript[S, t]) at  time t +  δt and q of going down to d (Subscript[S, t]). For certain choices of the values of p, q,u,d valuation based on the discrete random walk will converge as the number of time steps increases to the Black-Scholes model. There is are a number of possible choices of the values of p, q,u, and d for which the valuation given by the discrete walk model converges to the Black-Scholes one. The Benninga-paper implements just one such model, the  due to Cox and Rubinstein. This model has the property that p * q  =  1. Another well known model is due to Jarrow and Rudd. In this model both p and q are taken to be (1/2). There are also other models, described in the books of Hull and that of Willmott-Howison-Dewynne. All the texts, including Benninga's article and the books give a "probabilistic" justification for the particular choices of p,q,u and d. The idea is to compare the mean and the variance of the discrete random walk with the continuous model. However, Shaw introduces a  different, non-probabilistic methods. In his approach these values are chosen so as to preserve certain basic solutions of the Black-Scholes differential equation. (We shall study this equation and various forms of solutions to it later in the course.) This method makes it possible to recover four known approaches (which Shaw calls the CRR model, the approximate CRR model, the JR model, the approximate JR-model), but also allows him to develop new models, which he calls the "supersymmetric tree model" and the "approximate supersymmetric tree model". These models have some advantages over the standard CRR and JR models. In particular, in the CRR model it is possible for the "probabilities" to become negative and the model to fail. Similarly, in the case of the JRR model (where the "probabilities" are (1/2)) it is possible for the abscissas to become negative.  Shaw's super-symmetric model avoids this problem.  In Chapter 18 of his book Shaw uses Mathematica to derive all the models using his approach. In Chapter 19 he gives Mathematica implementations of  pricing for European and American options using each model and the Benninga's programming approach. 
Show's programs in Chapter 19 use exactly the same method as Benninga, except that in addition to the value of the option they include the "Greeks", \(TraditionalForm\`Δ, \ Γ, \ Θ\) (delta, gamma, theta), using (essentially) the method described on page 394 of Hull's book. Also, Shaw includes the case of options on shares with a continuously paid dividend. (The case of discontinuous dividends is a bit more complicated, and discussed for example in Hull's book).
There is one other  interesting thing in Chapter 19 of Shaw's book. In Benninga's article we find a beginning of the analysis of the convergence of the binomial model as the number of steps increases. The authors noticed that increasing  the number of steps does not automatically insure a more accurate value, instead the values one obtains converge to the correct value in an oscillating manner. Shaw has studies this phenomenon further and found that one can obtain faster convergence if one chooses steps in a particular way. The best way to choose the steps depends on the model used and the type of option.

Benninga's Algorithm

Benninga's algorithm for a European option is rather obvious and very fast one. As is well known, the value of the option is given by the sum of the "state prices" at the terminal nodes of the tree times the option values at these nodes.

Sum[
      exercise[S u^j d^(n - j)] Binomial[n, j] P^
          j Q^(n - j), {j, 0, n}]

where exercise is the exercise function. As a "pure" function the exercise function for a call option can be written as

Max[# - K, 0] &

and for a put option as

Max[K - #, 0] &

Much more interesting is Benninga's algorithm or an American option. It is a recursive approach which is, in a sense, the opposite of the "standard" backward inductive approach. In the standard approach one first computes the

TreeOp[node_, level_] := TreeOp[node, level] = 
          If[level == n, exercise[S d^node u^(level - node)], 
            Max[{P, Q} . {TreeOp[node, level + 1], 
                  TreeOp[node + 1, level + 1]}, 
              exercise[S d^node u^(level - node)]]];
value = TreeOp[0, 0];

The computation starts from the root TreeOp[0,0]. and then proceeds up the tree. At each node of each level below the last one a new rule is established giving the value of the option in terms of successive levels and nodes. Of course no actual computations can be done until the final level is reached. Instead a rule is created for computing the values at each level and node  once the values at the next level and node are known. Only when the last level is reached these values are computed. Note the way both := (SetDelayed) and = (Set) are used in the program:

TreeOp[node_, level_] := TreeOp[node, level] =.  ..

This technique is known as "dynamic programming". It makes Mathematica remember the values once computed so they need not be recomputed again if the same value is again needed. You can try the same program with only SetDelayed

TreeOp[node_, level_] :=  ...

and check that it will be considerably slower.

Problem 1

The weakness of this approach is that it uses a lot of memory (so called "Stack memory") and is relatively slow. Moreover, Mathematica by default has a limit on the length of recursion that can be performed (which helps to avoid getting into certain infinite loops). The default is:

$RecursionLimit

If  you want to use Benninga's method with very large trees you had better sent $RecursionLimit to  larger value. You can see that at the beginning of Chapter 19 of Shaw's book (in which he uses Benninga's approach) and sets $RecursionLimit to 5000. However, for larger trees it is much better to use a different method. This is the method used in Chapter 21 of Shaw's book. The package OptionPricingWithTrees uses Benninga's method when the option Tree is set to true, and a much faster method (similar to Shaw's Chapter 21) otherwise.  This faster method is discussed below.

Symmetry and Super-Symmetry

The price of an option obtained by using a binomial tree agrees with the Black-Scholes model provided the parameters u,d,p and q are chosen appropriately. Usually this choice is justified by means of a probabilistic argument. One considers a discrete random walk represented by a tree and compares it with the "continuous random walk" of Black-Scholes (lognormal).  One then chooses the parameters for the discrete random walk so that it has the same expectation and variance in an infinitesimal time step as the continuous lognormal random walk. There are well known logical problems involved in thinking about binomial trees in probabilistic terms, including the fact that for certain size of time step the probabilities can become negative. Shaw adopts another approach, which is to consider the Black-Scholes equation, find some special, simple solutions and then require that the random walk preserves them. The Black-Scholes equation is:

(1/2) (S^2) (((\[DifferentialD]^2) Val)/(\[DifferentialD](S^2))) (σ^2)-r Val+(r - q) S ((\[DifferentialD]Val)/(\[DifferentialD]S))+((\[DifferentialD]Val)/(\[DifferentialD]t))==0

where Val is the value of a derivative security at time t, S the price of the underlying at time, t, r is the risk free rate of interest, and q a continuous dividend. Shaw considers so called "power solutions" of the Black-Scholes equations:

Val = ((Z ((((S/K)))^n) (E^(((r - q) n + (1/2) n (n - 1) (σ^2)) (T - t))))/(E^(r (T - t))))

Here Z and n are constant real numbers (n is not necessarily an integer).

Problem 2

Problem 3

Shaw defines a "tree rule" as the following:

Val(S, 
        t) = (E^(-r Δ t)) (∑_(i = 1)^n) (Subscript[w, i]) Val(S (Subscript[u, i]), t + Δ t)

If n=2 we have a binary (or binomial) rule, for n=3 ternary (or trinomial) rule, etc. The requirement that the rule preserves a power solution of the B-S equation gives and equation relating the w's and the u's. By choosing different sets of solutions to be preserved Shaw obtains all the known binomial schemes as well as a new one (or rather two, an "exact" and an "approximate" one) which he calls "super-symmetric".

Problem 4

Explain, with an example,  what is meant by an "exact" and by an "approximate" tree scheme.

Problem 5

Using solutions of the Black-Scholes equations different form those used by Shaw find your own binomial or trinomial scheme and investigate some of its properties

Problem 6

Explain why  Shaw's approach leads to schemes for valuing options which converge to the Black-Scholes value (you do not need to give a formal proof, although one would of course be welcome!).

Accuracy of the Tree Method

Benninga and his co-authors after implementing their method in Mathematica started investigating its convergence properties and noticed rather disturbing large oscillations, in the values obtained as the number of time steps increases. These oscillations make it quite difficult to achieve accurate results without applying additional techniques. For example, the last part of  problem 1 of the first problem set asks to "Try to obtain the value of the option with three correct digits." Without using additional methods this is quite hard. For example, using Table to observe oscillations after 300 steps we still can't decide whether the third digit should be 3 or 44:

Table[AmericanPut[60, 60, i, 0.4, 5/12, 0.1], {i, 300, 310}]

When we use 1000 time steps we can assure ourselves that 4 is the right answer:

Table[
      AmericanPut[60, 60, i, 0.4, 5/12, 0.1], {i, 1000, 1010}]

Shaw observed that the most accurate values of the option are obtained by taking a time step so that the strike price lies in between two nodes. It is easy to derive the formulas for the straddle values for different binomial models. For example for the CRR model Shaw gives the formula:

straddleCRR[m_,  σ_,  T_,  S_,  K_]  :=  
      Round[(m + 1/2)^2 σ^2 T/Log[K/S]^2]

and for the JR model

straddleJR[m_,  σ_,  T_,  S_,  K_,  r_]  :=  
      Round[((m  +  1/2)^2*σ^2*
              T)/(Log[K/S]  -  (r - σ^2/2))^2]

Note that in the above example CRR straddle values can't be used since the denominator will become 0 while JR straddle values are:

Table[
      straddleJR[i
, 0.4,  5/12,  60, 60,  0.1], {i, 1, 
        10}]

Problem 7

Shaw does not implement or discuss the well known method to improve the accuracy of results known as "control variate" technique.

Problem 8

Implement the control variate technique for valuing options in a separate package which loads in the OptionPricingWithTrees package (or your own tree based package) and the BlackScholes package. See if you can get more accurate values for the problem above with shorter trees. Also, implement  two other methods to improve accuracy. One consists of taking the arithmetic mean (average) of the results obtained with successive numbers of time steps.  The other takes the geometric mean of successive iterations in the manner described by Benninga on page 264. Compare the accuracy of these two methods and the control variate technique on selected examples. .