Comments on  Shaw's Chapter 20.

Fast Option Pricing using trees

Benninga's algorithm for pricing American otpions is elegant but not efficient. This is due mostly to his using a recursive method. There is of course no need to do that. A non-recursive algortihm is given on page 190 of WHD,  in the form of  "pseudo-code" for valuing an American option.  This pseudo-code is written in a language similar to the popular procedural programming language C.

The algorithm for valuing an American option in WHD is basically the same as the one used in Chapter 20 of Shaw's book, although the implementation is somewhat different. Shaw divides his program into two parts. He first defines a function which perfomrs a single step of the "backward induction".  His second funciton basically repeats this inductive step the required number of times. 
He gives two different codes for his single step function. The first is C-like code very similar to the one given in WHD.

\(\(OneTimeStepAmer[initial_, \ payoff_, p_, 
          q_] := 
		Module[{wold\  = \ initial, 
            wnew\  = \ Rest[initial], \ len\  = \ Length[initial] - 1, \ 
            m}, 
			For[m = 1, \ 
            m <= \ len, \ \(m++\), 
				wnew[\([m]\)]\  = 
              Max[payoff[\([m]\)], \ \ p*wold[\([m]\)] + 
                  q*wold[\([m + 1]\)]]]; 
			wnew];\)\)

Shaw also gives a functional version of the single step code.

\(\(OneTimeStepAmerFun[initial_, \ payoff_, p_, q_] := 
		MapThread[
          Max, \ {payoff, \ 
            p*Drop[initial, \(-1\)]\  + \ 
              q*Drop[RotateLeft[initial], \ \(-1\)]}];\)\)

Of course there are lots of different ways to program this. Here are a few other examples:

OneTimeStepAmerFun1[initial_,  payoff_, p_, q_] := 
      Max /@ Transpose[ {payoff,  
            p*Drop[initial, -1]  +  
              q*Drop[RotateLeft[initial],  -1]}]

Here is another way to do this, whihc uses the built-in function ListCorrelate.

OneTimeStepAmerFun2[initial_,  payoff_, p_, q_] := 
      Max /@ Transpose[ {payoff, ListCorrelate[{p, q}, initial]}]

One we compare performance we see that all the functional programs outperform the procedural one.

Do[OneTimeStepAmer[{1, 1, 1},  {10,  0},  1/6,  1/100], {1000}] // 
      Timing

Do[OneTimeStepAmerFun[{1, 1, 1},  {10,  0},  1/6,  1/100], {1000}] // 
      Timing

Do[OneTimeStepAmerFun1[{1, 1, 1},  {10,  0},  1/6,  
          1/100], {1000}] // Timing

Do[OneTimeStepAmerFun2[{1, 1, 1},  {10,  0},  1/6,  
          1/100], {1000}] // Timing

One of the strong points of Mathematica as a tool for numericla compuations is the ability to compile certain programs. When a program is compiled it is converted to a pseudo-assemly like language that usually allows it to run faster. Here is the result of compiling the procedural and functional OntTimeStep functions:

\(CompiledOneTimeStepAmer\  = \ 
      Compile[{{initial, \ _Real, \ 1}, \ {payoff, \ _Real, \ 
              1}, \ {p, \ _Real, \ 0}, \ {q, \ _Real, \ 
              0}}, 
		Module[{wold\  = \ initial, 
              wnew\  = \ Rest[initial], \ len\  = \ Length[initial] - 1, \ 
              m}, 
			For[m = 1, \ 
              m <= \ len, \ \(m++\), 
				wnew[\([m]\)]\  = 
                Max[payoff[\([m]\)], \ \ p*wold[\([m]\)] + 
                    q*wold[\([m + 1]\)]]]; 
			wnew]] // 
        InputForm\)

\(\(\(CompiledOneTimeStepAmerFun\  = \ 
      Compile[{{initial, \ _Real, \ 1}, \ {payoff, \ _Real, \ 
              1}, \ {p, \ _Real, \ 0}, \ {q, \ _Real, \ 0}}, 
		MapThread[
            Max, \ {payoff, \ 
              p*Drop[initial, \(-1\)]\  + \ 
                q*Drop[RotateLeft[initial], \ \(-1\)]}]] // 
        InputForm\)\(
\)
    \)\)

You can see that oth functions are now much faster than before and the difference between them has basically disappeared. In general functional Mathematica programs run faster than procedural ones, but if both can be compiled they perfom roughly as fast.

Do[CompiledOneTimeStepAmer[{1, 1, 1},  {10,  0},  1/6,  
          1/100], {1000}] // Timing

Do[CompiledOneTimeStepAmerFun[{1, 1, 1},  {10,  0},  1/6,  
          1/100], {1000}] // Timing

The main problem wiht compiling is that quite many funcitons (and programs) fail to compile and sometimes when a program does compile if fails to become faster. But still compiling is a powerful technique, of great importnace in computationally intensive numerical programming.

The main difference between the AmericanPut and AmericanCall functions in Chapter 20 of Shaw's book and the one in my package OptionPricingWithTrees is that the first one is basically procedural (even the function that Shaw calls functional is in part procedural, while my code is essentially functional. However, becasue both codes are compiled there is no significant speed difference between them.)

This is how Shaw defines the function to compute the value of an AmericanPut.  Note that in order to compute the Greek parameters the function is somwhat more complicated than it woudl be if oe only wanted to compute the option value.

ManyTimeStepAmerPut = 
        Compile[{{initial, _Real, 1}, {p, _Real, 0}, {q, _Real, 
              0}, {S, _Real, 0}, {u, _Real, 0}, {d, _Real, 0}, {n, _Integer, 
              0}, {K, _Real, 0}}, 
          Module[{new = initial, k, r, old = initial, vold = initial, payoff, 
              finpay = {Max[K - S, 0]}}, 
            For[k = 1, 
              k ≤ 
                n - 1, k++, (payoff = 
                  Table[Max[K - S d^r u^(n - k - r), 0], {r, 0, 
                      n - k}]; 
old = new; 
                new = OneTimeStepAmer[old, payoff, p, 
                    q])]; 
vold = old; 
            old = new; 
new = 
              OneTimeStepAmer[old, finpay, p, q]; 
{Join[
                new, old], vold}]];

CRRFastAmerPut[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{u = CRRupBinApp[n,r,q,σ,T], d = CRRdownBinApp[n,r,q,σ,T],
a = afunc[n,r,q,T],
disc = Discount[n,r,T],
P,Q,TreeOp,result,payoffin},
P = Pfunc[u,d,a,disc];
Q = Qfunc[u,d,a,disc];
payoffin = Table[Max[K-S d^node u^(n-node),0], {node, 0, n, 1}];
result = ManyTimeStepAmerPut[payoffin,P,Q,S,u,d,n,K][[1,1]];
]

Here is how the same thing is done in my OptionPricingWithTrees. For simplicity I have only included the code to compute the option value. THe actual function in the package also computes the Greek parameters.

CompiledAmerPut = 
        Compile[{{initial, _Real, 1}, {p, _Real, 0}, {q, _Real, 
              0}, {S, _Real, 0}, {u, _Real, 0}, {d, _Real, 0}, {K, _Real, 
              0}}, Last[
            Nest[Max /@ 
                  Transpose[ {Table[
                        K - Max[S d^r u^(Length[#] - 2 - r), 0], {r, 0, 
                          Length[#] - 2}],  
                      p*Drop[#, -1]  +  
                        q*Drop[RotateLeft[#],  -1]}] &, initial, 
              Length[initial] - 1]]
];

FastAmericanPut[K_, S_, n_, σ_, T_, Rf_] := 
      Module[{u = up[n, σ, T], d = down[n, σ, T], 
          r = R[n, Rf, T], p, q}, p = P[u, d, r]; 
        q = Q[u, d, r]; 
initial = 
          Table[Max[K - S d^node u^(n - node), 0], {node, 0, 
              n}]; 
CompiledAmerPut[initial, p, q, S, u, d, 
          K]]

The main difference between the function is that Nest is used instead of the For loop and there are no assignments (=) in CompiledAmericanPut.  This is in line with the functional programming style., As in the case of a single step, the uncompiled functional program is considerable faster than the procedural one but the difference largely disappears as a result of compiling.

Problem 8

Try to implement the code  from page 190 of WHD  in Mathematica. Compare the correctness and efficiency of your code with that of Benninga's code. Try to compile your code.  Compare its performance with un-compiled version.

Problem 9

On page 396 of chapter 20 of Shaw's book there compiled code for a trinomial scheme for valuing European options. Write a similar compiled CRR scheme for valuing American call and put options.

Problem 10

Try to modify the code for a trinomial tree to use the adaptive mesh model described on page 406 of Hull's book.

Problem 11 (very difficult, I think)

The slowest function in Benninga article is the AsianOption function designed to value path dependent options by means of binary trees and the functions AmerAvgCall and AmerAvgPut, , EurAvgCall and EurAvgPut, which depend on it. Most of the slowness of this function is due to the exponential complexity of the algorithm. Try to improve the speed of the functions AmerAvgCall and AmerAvgPut, EurAvgCall and EurAvgPut, by compiling them or parts of them.