A New Binomial Model for Discrete Lookback Options using Mathematica

Dr Michael Kelly and Dr James Brown

Lookback options are a special case of exotic derivatives where the value of the derivative depends upon some function of all of the previous underlying asset(s) prices. In particular it attempts to mimic the optimal situation where one can purchase stocks at the minimum price and sell stocks at the maximum price, ie buy at the low and sell at the high. There are two types of call and put Lookback options, depending upon whether the strike or the final asset price is replaced by the optimal asset price. Using the standard notation where S is the asset and K is the strike price that one pays for the asset at some future date, then the Floating strike or Standard Lookback call has the value (Subscript[S, T]) - 
        Min((Subscript[S, t]),  0 ≤ t ≤ T), whereas the Fixed strike Lookback call or call on the Max has the value Max(0, 
        Max((Subscript[S, t]),  0 ≤ t ≤ T) - K). Similarly the Floating strike put or Standard Lookback put has the value Max((Subscript[S, t]),  0 ≤ t ≤ T) - 
        (Subscript[S, T]), whereas the Fixed strike Lookback put or put on the Min has the value Max(0, 
        K - Min((Subscript[S, t]),  0 ≤ t ≤ T)). Standard Lookback options were first evaluated by Goldman, Sosin and Gatto (1979) and probabilistic models for all types of Lookback options were later calculated by Conze and Viswanathan (1991). All of these models were determined within the continuous framework, where the value of the underlying security is modeled as changing continuously with time. However such prices were found by Babbs (2000) and others to significantly overprice the true value of the market Lookbacks, because the actual Lookback options were determined from asset prices valued at discrete intervals of a fraction of a day. This paper derives two new methods for determining discrete Lookback options by using Mathematica's powerful symbolic capabilities to find simplified expressions for the optimal asset prices in each aseet path in a binomial tree model.

Initialisations

The implementation of descriptive functions and graphics in the paper will require pre-defined functions contained in Mathematica Add-On packages from the contexts of Statistics and Graphics.

Needs["Statistics`"];  Needs["Graphics`Graphics`"];  
    Needs["Graphics`MultipleListPlot`"];
Off[General::"spell1"]; Off[General::"spell"];

Previous Lookback Option models

Conze-Viswanathan Lookback Option models

Conze and Viswanathan (1991) used probabilistic theorems about the joint distributions of stochastic processes with their maxima and minima taken from Harrison (1985) to evaluate the probabilities for maxima (Subscript[M, (Subscript[t, 1])]^((Subscript[t, 2]))) and minima (Subscript[m, (Subscript[t, 1])]^((Subscript[t, 2]))) of stock prices (Subscript[S, t]) in terms of the normal cumulative distribution function 𝒩. These resulted in formulae for the Call on the Max and Put on the Min which are coded by Shaw (1998). In order to describe these results we first define auxiliary parameters based upon the usual Black-Scholes market parameter (Subscript[d, T]) and its extensions \(TraditionalForm\`\(d\_T\%'\)\) and \(TraditionalForm\`\(d\_T\%''\)\) to the cases when the stike price K is replaced by the price maxima and minima respectively.

\(TraditionalForm\`\(\(
\)\(d\_T = \ \((ln(
                    S\_0/K) + \((r - 
                        q + σ\^2/
                          2)\)\ T)\)/\((σ\ \@T)\)

          \(d\_T\%'\) = \ \((ln(
                    S\_0/M\_\(T\_0\)\%0) + \((r - 
                        q + σ\^2/
                          2)\)\ T)\)/\((σ\ \@T)\)

          \(d\_T\%''\) = \ \((ln(
                    S\_0/m\_\(T\_0\)\%0) + \((r - 
                        q + σ\^2/
                          2)\)\ T)\)/\((σ\ \@T)\)

          
          M\_\(t\_1\)\%\(t\_2\) = 
            max(S\_t, \ t ∈ \([t\_1, \ t\_2]\))

          
          m\_\(t\_1\)\%\(t\_2\) = 
            min(S\_t, \ t ∈ \([t\_1, \ t\_2]\))\)\)\)

The formulae for the Call on the Max, C and the Put on the Min, P are described below

C=(Piecewise[{{\(\(S\_0\) \(𝒩(
                            d\_T)\) - \(e\^\(\(-r\)\ T\)\) 
                          K\ \(𝒩(
                            d\_T - σ \@ T)\) + \(e\^\(\(-r\)\ T\)\) \(σ\^2\/\(2  
                                r\)\) \(\(S\_0\)(\(e\^\(r\ T\)\) \(𝒩(d\_T)\) - \(\((S\_0\/K)\)\^\(-\(\(2  
                                        r\)\/σ\^2\)\)\) \(𝒩(d\_T - \(\(2  
                                        r\)\/σ\) \@T)\))\)\), K≥(Subscript[M, (Subscript[T, 0])]^0) }, {\(\(e\^\(\(-r\)\ T\)\)(
                          M\_\(T\_0\)\%0 - 
                            K) + \(S\_0\) \(𝒩(\(d\_T\%'\))\) - \(e\^\(\(-r\)\ T\)\) 
                          M\_\(T\_0\)\%0\ \(𝒩(\(d\_T\%'\) - σ \@ T)\) + \(e\^\(\(-r\)\ T\)\) \(σ\^2\/\(2  
                                r\)\) \(S\_0\) \((\(e\^\(r\ T\)\) \(𝒩(\(d\_T\%'\))\) - \(\((S\_0\/M\_\(T\_0\)\%0)\)\^\(-\(\(2  
                                        r\)\/σ\^2\)\)\) \(𝒩(\(d\_T\%'\) - \(\(2  r\)\/σ\) \@T)\))\)\), K < 
                        (Subscript[M, (Subscript[T, 0])]^0)}}])

P=(Piecewise[{{\(\(-S\_0\) \(𝒩(\(-d\_T\))\) + \(e\^\(\(-r\)\ T\)\) 
                          K\ \(𝒩(\(-d\_T\) + σ \@ T)\) + \(e\^\(\(-r\)\ T\)\) \(σ\^2\/\(2  
                                r\)\) \(\(S\_0\)(\(-e\^\(r\ T\)\) \(𝒩(\(-d\_T\))\) + \(\((S\_0\/K)\)\^\(-\(\(2  
                                        r\)\/σ\^2\)\)\) \(𝒩(\(-d\_T\) + \(\(2  
                                        r\)\/σ\) \@T)\))\)\), K≤(Subscript[m, (Subscript[T, 0])]^0) }, {\(\(e\^\(\(-r\)\ T\)\)(
                          K - m\_\(T\_0\)\%0) - \(S\_0\) \(𝒩(\(-\(d\_T\%''\)\))\) + \(e\^\(\(-r\)\ T\)\) 
                          m\_\(T\_0\)\%0\ \(𝒩(\(-\(d\_T\%''\)\) + σ \@ T)\) + \(e\^\(\(-r\)\ T\)\) \(σ\^2\/\(2 
                                 r\)\) \(\(S\_0\)(\(-e\^\(r\ T\)\) \(𝒩(\(d\_T\%''\))\) + \(\((S\_0\/m\_\(T\_0\)\%0)\)\^\(-\(\(2  r\)\/σ\^2\)\)\) \(𝒩(\(d\_T\%''\) - \(\(2  
                                        r\)\/σ\) \@T)\))\)\), K > 
                        (Subscript[m, (Subscript[T, 0])]^0)}}])

\(\(\(\(d1[σ_, r_, q_, x_, y_, 
            t_] := \((Log[
                  x/y] + \((r - 
                      q + σ\^2/2)\)\ t)\)/\((σ\ \@t)\);\)

      \(NormCDF[z_] := \((1 + Erf[z/\@2])\)/2\ ;\)\)\(\ \)\)\)
done[p_, pmax_, sd_, r_, q_, t_] := d1[sd, r, q, p, pmax, t];
dtwo[p_, pmin_, sd_, r_, q_, t_] := d1[sd, r, q, p, pmin, t];

CVPutOnMin[p_, pmin_, k_, r_, q_, sd_, t_] :=
If[k <= pmin,
-p*NormCDF[-d1[sd, r, q, p, k, t]] + 
Exp[-r*t]*k*NormCDF[-d1[sd, r, q, p, k, t] + sd*Sqrt[t]] + 
Exp[-r*t]*(sd^2/(2*r))*p*(
(p/k)^(-2r/sd^2)*NormCDF[-d1[sd, r, q, p, k, t] + 2*r*Sqrt[t]/sd] - 
Exp[r*t]*NormCDF[-d1[sd, r, q, p, k, t]]),
Exp[-r*t]*(k - pmin) - 
p*NormCDF[-dtwo[p, pmin, sd, r, q, t]] + 
Exp[-r*t]*pmin*NormCDF[-dtwo[p, pmin, sd, r, q, t] + sd*Sqrt[t]] + 
Exp[-r*t]*(sd^2/(2*r))*p*(
(p/pmin)^(-2r/sd^2)*NormCDF[-dtwo[p, pmin, sd, r, q, t] + 2*r*Sqrt[t]/sd] - 
Exp[r*t]*NormCDF[-dtwo[p, pmin, sd, r, q, t]])]

CVCallOnMax[p_, pmax_, k_, r_, q_, sd_, t_] :=
If[k >= pmax, 
p*NormCDF[d1[sd, r, q, p, k, t]] - 
Exp[-r*t]*k*NormCDF[d1[sd, r, q, p, k, t] - sd*Sqrt[t]] + 
Exp[-r*t]*(sd^2/(2*r))*p*(
-(p/k)^(-2*r/sd^2)*NormCDF[d1[sd, r, q, p, k, t] - 2*r*Sqrt[t]/sd] + 
Exp[r*t]*NormCDF[d1[sd, r, q, p, k, t]]),
Exp[-r*t]*(pmax - k) + 
p*NormCDF[done[p, pmax, sd, r, q, t]] - 
Exp[-r*t]*pmax*NormCDF[done[p, pmax, sd, r, q, t] - sd*Sqrt[t]] + 
Exp[-r*t]*(sd^2/(2*r))*p*(-(p/pmax)^(-2*r/sd^2)*NormCDF[done[p, pmax, sd, r, q, t]
- 2*r*Sqrt[t]/sd] + 
Exp[r*t]*NormCDF[done[p, pmax, sd, r, q, t]])]

We will use the above code to compare with our discrete valued functions. In doing so it will be noted, as with Babbs (2000) that the continuous functions of Conze and Viswanathan overestimate the true value of lookback options.

Stochastic and Binomial Paths

Using the NormalDistribution package we can generate values for a Gaussian distribution with the function RandomArray. This approach can then be used to create random paths for equities by using the lognormal distribution structure of the updating process between successive values of the equities, where Z ~ 𝒩(0,1) is the standard normal random variable.

\(TraditionalForm\`S\_\(t + Δt\) = 
            S\_t\ e\^\(\((r - q - σ\^2\/2)\)\ Δt + Z\ σ\ \@Δt\)\)

Equation (4) allows the determination of new stochastic values (Subscript[S, t + Δt]) from the previous ones (Subscript[S, t]) by observing that the elements in the path are lognormally distributed with mean (r - q - ((σ^2)/2)) Δt and standard devistion σ (Sqrt[Δ])t where the time increment Δt = t/(navg - 1) is determined by the number of points dispersed along each path up to the final date t. The function coarsepath utilises equation (4) and the random number generator to create 4 n paths for an equity with risk free interest rate r, dividend yield q, volatility σ, starting from an initial value of s and updated at a succession of navg points up to time t. An example is given in the plot below:

coarsepath[s_, r_, q_, σ_, t_, n_, navg_] :=
Module[
{mpath = N[(r-q-σ^2/2)*t/(navg - 1)], 
spath = N[σ Sqrt[t/(navg - 1)]]},
Flatten[
 Table[
  Transpose[
   NestList[# Exp[RandomArray[NormalDistribution[mpath,spath],4]]&, {s, s, s, s}, navg - 1]], 
  {i, n}],
 1]
];

example = 
        coarsepath[100, 0.1, 0, 0.2, 1, 5, 30];
With[{len  = Length[example]},  
        MultipleListPlot[example, SymbolShape -> None, 
          PlotJoined -> True,  
          PlotStyle -> Table[GrayLevel[i/len], {i, 0, len - 1}], 
          AxesOrigin -> {0, Automatic}, 
          PlotRange -> All ] ];

The plot of stock paths in Fig.1 above illustrates the stochastic nature of the underlying equity. What we need is a tree-like pattern of asset growth whose variance spreads out linearly with T that will expand in the same manner as Fig. 1. The first two steps in the asset and option paths are displayed in Fig. 2 below, while a more detailed tree pattern approximating that shown above is given in Fig. 3.

networkPlot[points_, labels_, lines_, r_, opts___] := 
      Show[Graphics[
{Map[Line, 
              Map[{points[[#[[1]]]], points[[#[[2]]]]} &, 
                lines]], {GrayLevel[1], 
              Map[Disk[#, Offset[{r, r}]] &, points]}, 
            Map[Circle[#, Offset[{r, r}]] &, points], 
Map[
              Text[#[[1]], #[[2]]] &, Transpose[{labels, points}]]}, 
          opts]]

The function BinomialConnectedOptionPlot is similar to BinomialConnectedStockPlot but is concerned with showing how the corresponding option prices change with different values of S and is illustrated in Fig. 2 below.

BinomialConnectedStockPlot[S_, u_, d_, n_Integer] := 
        With[{binPoints = 
              Flatten[Table[{i, j}, {i, 0, n}, {j, n - i, n + i, 2}], 1],  
            binLabels = 
              Flatten[Table[
                  S (d^((n - j + i)/2)) (u^((j - n + i)/2)), {i, 
                    0, n}, {j, n - i, n + i, 2}]], 
            binLines = 
              Flatten[Table[{{i (i + 1)/2 + 
                        k, (i + 1) (i + 2)/2 + k}, {i (i + 1)/2 + 
                        k, (i + 1) (i + 2)/2 + k + 1}}, {i, 0, 
                    n - 1}, {k, 1, i + 1}], 2]}, 
          networkPlot[binPoints, binLabels, binLines, 20, 
            PlotRange -> {{-0.5, n + 0.5}, {-1, 2  n + 1}}, 
            TextStyle -> {FontSize -> 11}, 
            DisplayFunction -> Identity]];

BinomialConnectedOptionPlot[V_, u_, d_, n_Integer] := 
        With[{binPoints = 
              Flatten[Table[{i, j}, {i, 0, n}, {j, n - i, n + i, 2}], 1],  
            binLabels = 
              Flatten[Table[
                  Subscript[V,  
                    (d^((n - j + i)/2)) (u^((j - n + i)/2))], {i, 
                    0, n}, {j, n - i, n + i, 2}]], 
            binLines = 
              Flatten[Table[{{i (i + 1)/2 + 
                        k, (i + 1) (i + 2)/2 + k}, {i (i + 1)/2 + 
                        k, (i + 1) (i + 2)/2 + k + 1}}, {i, 0, 
                    n - 1}, {k, 1, i + 1}], 2]}, 
          networkPlot[binPoints, binLabels, binLines, 20, 
            PlotRange -> {{-0.5, n + 0.5}, {-1, 2  n + 1}}, 
            TextStyle -> {FontSize -> 11}, 
            DisplayFunction -> Identity]];

Consider the situation where the are only two possible states of the financial universe for our instruments: at the next moment of time, eithe asset S go up by some fixed amount u to S*u or go down by an amount d to the value S*d. At the same time the option V which is a function of S eithe goes up to (Subscript[V, u]) or down to (Subscript[V, d]) (note that these are subscripts, not multiples of u and d). To see what is happening consider Fig. 2 below which shows the state of our ideal market after two steps into the future.

Show[
        GraphicsArray[{BinomialConnectedStockPlot[S, u, d, 2], 
            BinomialConnectedOptionPlot[V, u, d, 2]}]];

Show[BinomialConnectedStockPlot[S, u, d, 4], 
        DisplayFunction -> $DisplayFunction];

Fig. 3 is a graphical description of a binomial model where the up and down increments in each step are reciprocal to each other. This is the well known Cox, Ross and Rubinstein (1979) model, generally referred to by its acronym CRR.

There already exist algorithms by Goldman, Sasin and Gatto (1979) and Conze and Viswanathan (1991) for the evaluation of continuous Lookback options which assume that the underlying equity can be measured continually. However the reality is that most Lookbacks are only determined daily and that consequently the continuous values significantly exceed (Babbs, 2000) the real values. Attempts to use discrete models have floundered because they are slow. But in the following sections we shall show that by identifying the hidden structure for optimal prices in binomial CRR trees it is possible to find accurate values which are just as fast as the analytical continous formulae.

Hull and Babb's method - Compiled binomial tree

Hull (2003) reviewed the difficulties associated with the calculation of lookback options with binomial trees and examined one method for improving the computation for a floating lookback put.  The essential difficulty in all lookback options is that not only must the price be computed at each node but the maximum or minimum price achieved so far must also be computed.  The result is a rapid increase in the computation time with tree size.  For the floating lookback put, Hull constructs a tree based not on the price of the stock, but on the ratio of the highest stock price achieved so far to the current stock price.  This is effectively a change of numeraire to the current stock price, which reduces the payoff (in the new numeraire) to the excess of the ratio over 1.  Although the geometry of the new tree is not the same as a conventional binomial tree, for example the ratio starts at 1 and can never decrease, its parameters are closely related to those for the standard CRR binomial tree.  The method does not make an allowance for a previous maximum price in excess of the price at the time of valuation; it also does not address the closely related call on the max option.

Babbs gives an efficient method for the valuation of standard lookback options. The trick is to use the stock as a numeraire. The tree is still recombining if CRR parameters are used and it can be easily modified to accommodate early exercise. This method is proposed by Reiner and Babbs independently. Hull also presents this method in his book, but he only shows a simple case where we are pricing a new option, so the current max/min price is equal to the spot price. Babbs extends his method to allow max/min prices that are different from the current spot price.

Babbs (2000) develops a method similar to Hull, but also investigates discrete sampling of the maximum price, rather than only looking at continuous sampling.  Babbs adopts as his single stochastic variable the logarithm of the ratio maxS/S of the highest price achieved so far to the stock price, essentially the same change of numeraire as Hull.  The restriction on the possible value of the new variable is that it must be greater than zero, but rather than reconstructing the tree to the special case of starting at zero, a reflecting barrier is inserted at zero. This recreates Hull's structure in the special case, but also extends the method to include a prior maximum.  Babbs then extends the tree to replicate discrete sampling by only allowing reflections at the sampling times, not continuously along the tree.

The following pieces of code are taken from Shaw (1998) for the description of the essential parameters in a CRR binomial model. They give respectively for each binomial step the up and down reciprocal coefficients, the drift factor (afunc), the discounted probabilities for up and down movements (Pfunc and Qfunc) and the discount factor. These shall be repeatedly used to describe the binomial tree pictured in Fig. 1.

CRRupBinApp[n_, r_, q_, σ_, T_] := 
        Exp[σ Sqrt[T/n]];
CRRdownBinApp[n_, r_, q_, σ_, T_] := 
        Exp[-σ Sqrt[T/n]];
afunc[n_,  r_,  q_,  T_]  :=  
        N[Exp[(r - q) T/n]];
Pfunc[up_, down_, a_, disc_] := 
        N[disc*(a - down)/(up - down)];
Qfunc[up_, down_, a_, disc_] := 
        N[disc - Pfunc[up, down, a, disc]];
Discount[n_, r_, T_] := N[Exp[-r T/n]];

The above auxiliary functions allow a CRR model for the evaluation of the Lookback put as determined by Babbs (2000). In the code maxS is actually the standardised maximum value divided by the current price maxS/S, so that its current value at time 0 would be 1. If the current ratio is maxS then at the next moment in time ie level+1 it will be either maxS-1 or the maximum of u maxS in the case of a downward jump with probability Q d or max(d maxS, 1) in the case of an upwards jump with probability P u . This is described in the recursive definition given below.

BabbsAmericanLookbackPut[S_, M0_, r_, q_, σ_, T_, 
          n_]:=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, value, delta, gamma, 
              theta},
P = Pfunc[u, d, a, disc];
Q = Qfunc[u, d, a, disc];
TreeOp[maxS_,level_]:=TreeOp[maxS, level]=If[level == n,maxS - 1,Max[P u TreeOp[Max[d maxS, 1],level + 1]+Q d TreeOp[u maxS, level + 1],maxS - 1]];
value=M0 TreeOp[1,0];
Clear[TreeOp];
value];

BabbsAmericanLookbackPut[50, 50, 0.2, 0.0, 0.25, 1.0, 100]

The above code can then be compiled to yield a version that can determine results which can be used in practise to calculate European discrete Lookback options. Here m0 refers to the minimum historical value achieved so far and M0 refers to the maximum historical value achieved so far.

BabbsLookbackCallCompiled = 
        Compile[{{S, _Real}, {m0, _Real}, {r, _Real}, {q, _Real}, {σ, _Real}, {T, _Real}, {n, _Integer}}, 
Module[{dt, d, p, 
              payoff, Mt, k, f, i, j, wnew, wold}, 
dt = 
              T/n; 
d = -σ*
                Sqrt[dt]; 
p = (Exp[(q - r)*dt] - 
                    Exp[d])/(Exp[-d] - 
                    Exp[d]); 
Mt = 
              Log[m0/S]; 
k = 
              IntegerPart[Mt/d]; 
f = 
              FractionalPart[Mt/d]; 
payoff = 
              Table[1 - 
                  Exp[If[i ≥ n - k, Mt - (n - 2*i)*d, 
                      If[Mod[n - i - k, 2] == 0, (i + f)*
                          d, (i + 1 - f)*d]]], {i, 0, 
                  n}]; 
wnew = 
              payoff; 
Do[
              wold = wnew; 
Do[
                wnew[[i + 1]] = 
                  Exp[-q*
                        dt]*(p*
                          If[i ≤ j - 1 - k && i > 0, 
                            wold[[i]], 
                            wold[[i + 
                                1]]] + (1 - p)*
                          wold[[i + 
                              2]]), {i, 0, 
                  j - 1}], 
{j, n, 
                1, -1}]; 
wnew[[1]]*S]];

BabbsLookbackPutCompiled = 
        Compile[{{S, _Real}, {M0, _Real}, {r, _Real}, {q, _Real}, {σ, _Real}, {T, _Real}, {n, _Integer}}, 
Module[{dt, u, p, 
              payoff, Nt, k, f, i, j, wnew, wold}, 
dt = 
              T/n; 
u = σ*
                Sqrt[dt]; 
p = (Exp[(q - r)*dt] - 
                    Exp[-u])/(Exp[u] - 
                    Exp[-u]); 
Nt = 
              Log[M0/S]; 
k = 
              IntegerPart[Nt/u]; 
f = 
              FractionalPart[Nt/u]; 
payoff = 
              Table[Exp[
                    If[i ≤ k, +(n - 2*i)*u, 
                      If[Mod[i - k, 2] == 0, (n - i + f)*
                          u, (n - i + 1 - f)*u]]] - 1, {i, 0, 
                  n}]; 
wnew = 
              payoff; 
Do[
              wold = wnew; 
Do[
                wnew[[i + 1]] = 
                  Exp[-q*
                        dt]*(p*
                          wold[[i + 
                              1]] + (1 - p)*
                          If[i ≥ k && i < j - 1, 
                            wold[[i + 
                                3]], 
                            wold[[i + 
                                2]]]), {i, 0, 
                  j - 1}], 
{j, n, 
                1, -1}]; 
wnew[[1]]*S]];

Note also that the European option values are smaller than the corresponding American option prices, as expected from theory.

BabbsLookbackPutCompiled[50, 50, 0.2, 0.0, 0.25, 1.0, 100] // 
      Timing

The following table shows that the Babbs calculation approaches the continuous values of Hull and Conze and Viswanathan from below.

ListPlot[
        Table[{n, 
            BabbsLookbackPutCompiled[50, 50, 0.2, 0.0, 0.25, 1.0, n]}, {n, 
            100, 500, 50}], PlotJoined -> True];

Adapted CRR Binomial model

Before developing our own model we will first construct a model which uses the CRR binomial structure in the recursive manner introduced by Benninga and Wiener (1997) and applied to barrier options in chapter 19 of Shaw's (2000) text. The price afterj down jumps at time level i is (Subscript[S, i, j]) = (Subscript[S, 0]) (d^j) (u^(i - j)). Let (Subscript[M, i, j]) denote the maximum price up to node (i,j). One of the possible values of (Subscript[M, i, j]) is (Subscript[S, i, j]) itself. Another possible value is S (d^(j - 1)) (u^(i - j)), which corresponds to a down move from the previous time level. Using this relationship recursively, we can list all the possible values of (Subscript[M, i, j]) as S (d^0) (u^(i - j)), 
      S (d^1) (u^(i - j)),  ... , S (d^j) (u^(i - j)). And of course (Subscript[M, i, j]) shouble satisfy (Subscript[M, i, j]) ≥ (Subscript[M, 0, 0]).

(Subscript[S, i, j]) = (Subscript[S, 0]) (d^j) (u^(i - j))
(Subscript[M, i, j]) = max((Subscript[S, i, j]) , 1≤j≤i≤n)

At the next time level, the price will either move up to (Subscript[S, i + 1, j]) = S (d^j) (u^(i + 1 - j)) or move down to (Subscript[S, i + 1, j + 1]) = 
        S (d^(j + 1)) (u^(i - j)). If it is an up move, the maximum price will be (Subscript[M, i + 1, j]) = 
        max((Subscript[S, i + 1, j]), (Subscript[M, i, j])). If it is a down move, the maximun price will be unchanged at (Subscript[M, i + 1, j + 1]) = (Subscript[M, i, j]). The option value f() at (i, j) given (Subscript[M, i, j]) is therefore

f ((Subscript[S, i, j]) | (Subscript[M, i, j]))=(e^(-r Δt)) (p f(
                        (Subscript[S, i + 1, j]) | (Subscript[M, i + 1, j])) + (1 - 
                          p) f(
                        (Subscript[S, i + 1, j + 1]) | 
                          (Subscript[M, i + 1, j + 1]))  )
\(\(=\)\(\(e\^\(\(-r\)\ Δt\)\) 
                  p\ \(f(S\ \(d\^j\) u\^\(i + 1 - j\) | 
                      max(\ S\ d\^j\ u\^\(i + 1 - j\), 
                        M\_\(i, j\))\ )\) + \(\(e\^\(\(-r\)\ Δt\)\)(1 - p)\)\ \(f(S\ \(d\^\(j + 1\)\) u\^\(i - j\) | M\_\(i, j\))\)\)\)

At maturity, the option value at node (N, j) given (Subscript[M, N, j]) is the payoff f((Subscript[S, N, j]) | (Subscript[M, N, j])) which is now given by:

f((Subscript[S, N, j]) | (Subscript[M, N, j])) = 
            max((Subscript[M, N, j]) - K, 0)

Unfortunately applying the above recursive technique we get unacceptably slow code. Here is a typical but slow binomial model to determine the European call on the Max Lookback which uses a recursive procedure adapted from the CRR model. As before the symbols have the following meanings: S is the current equity price, M0 is the previous maximum value of the equity before the initiation of the option, K is the strike price for which the buyer of the call option can purchase the underlying equity at the future maturity date T. The market parameters are the interest rate r, the dividend yield q and the volatility or standard devistion of the equity returns, σ. Finally the size of the binomial tree is determined by n, the number of steps in the tree that are required to span the time to maturity T.

CRRAppEuropeanCallOnMax[S_, M0_, K_, r_, q_, σ_, T_, n_] := 
        Module[{dt, a, u, d, P, Q, TreeOp, value}, 
dt = 
            T/n; 
a = 
            Exp[(r - q)*dt]; 
u = 
            Exp[σ*Sqrt[dt]]; 
d = 
            Exp[-σ*
                Sqrt[dt]]; 
P = (a - d)/(u - 
                  d); 
Q = 1 - P; 
P = 
            P*Exp[-r*dt]; 
Q = 
            Q*Exp[-r*dt]; 
TreeOp[node_, level_, 
              M_] := TreeOp[node, level, M] = 
              If[level == n, 
                Max[M - K, 
                  0], {P, 
                    Q} . {TreeOp[node, level + 1, 
                      Max[S*d^node*u^(level + 1 - node), M]], 
                    TreeOp[node + 1, level + 1, 
                      M]}]; 
value = 
            TreeOp[0, 0, M0]; 
Clear[
            TreeOp]; 
value];

This is the compiled procedural version of the above recursive model

CRRAppEuropeanCallOnMaxCompiled = 
        Compile[{{S, _Real, 0}, {M0, _Real, 0}, {K, _Real, 0}, {r, _Real, 
              0}, {q, _Real, 0}, {σ, _Real, 0}, {T, _Real, 
              0}, {n, _Integer, 0}}, 
Module[{dt, a, u, d, 
              P, Q, k, kmax, j, jmax, Smax, Mold, Mnew, fold, fnew, i, posu, 
              posd}, 
dt = T/n; 
a = 
              Exp[(r - q)*dt]; 
u = 
              Exp[σ*Sqrt[dt]]; 
d = 
              Exp[-σ*
                  Sqrt[dt]]; 
P = (a - d)/(u - 
                    d); 
Q = 
              1 - P; 
P = 
              P*Exp[-r*dt]; 
Q = 
              Q*Exp[-r*dt]; 
 (* 
              all possible max values and maturity and payoff *) 
jmax = 
              n; 
kmax = (jmax + 1 + Mod[jmax + 1, 2])/
                2; 
Mnew = 
              Table[0.0, {j, 0, jmax}, {k, 1, 
                  kmax}]; 
fnew = 
              Mnew; 
Do[k = 1; 
              Smax = S*u^(jmax - j); 
If[Smax < M0, 
                Mnew[[j + 1, 1]] = 
                  M0, 
While[
                  k ≤ kmax && 
                    Smax ≥ S*u^(n - j)*d^j, 
                  If[Smax > M0, 
                    Mnew[[j + 1, k]] = 
                      Smax; Smax *= d; k++, 
                    Mnew[[j + 1, k]] = 
                      M0; Break[];]]], {j, 0, jmax}]; 
Do[
              If[Mnew[[j + 1, 
                    k]] ≠ 0.0, 
                fnew[[j + 1, k]] = 
                  Max[Mnew[[j + 1, k]] - 
                      K, 0]], {j, 0, jmax}, {k, 1, 
                kmax}]; 
 (* 
              going backward through time *) 
Do[
Mold = Mnew; 
fold = 
                fnew; 
Do[
                Mnew[[j + 1, k]] = 
                  0.0, {j, 0, jmax}, {k, 1, kmax}]; 
fnew = 
                Mnew; 
jmax -= 
                1; 
kmax = (jmax + 1 + 
                      Mod[jmax + 1, 2])/2; 
 (* 
                possible max values *) 
Do[k = 1; 
                Smax = S*u^(jmax - j); 
If[
                  Smax ≤ M0, 
                  Mnew[[j + 1, 1]] = 
                    M0, 
While[
                    k ≤ kmax && 
                      Smax ≥ S*u^(jmax - j)*d^j, 
                    If[Smax > M0, 
                      Mnew[[j + 1, k]] = 
                        Smax; Smax *= d; k++, 
                      Mnew[[j + 1, k]] = 
                        M0; Break[];]]], {j, 0, 
                  jmax}]; 
 (* 
                discounted risk neutral expectation *) 
Do[
                Do[If[Mnew[[j, 
                        k]] ≠ 
                      0, 
 (* 
                      search the table for the corresponding max value at the next time level *) 
posu = 1; 
                    While[Mold[[j, 
                          posu]] ≠ 
                        Max[S*u^(jmax - j + 2)*d^(j - 1), 
                          Mnew[[j, 
                            k]]], posu++]; 
posd = 1; 
                    While[Mold[[j + 1, 
                          posd]] ≠ 
                        Mnew[[j, 
                          k]], posd++]; 
fnew[[j, k]] = 
                      P*fold[[j, 
                            posu]] + 
                        Q*fold[[j + 1, 
                            posd]]], 
{k, 1, kmax}], {j, 1, jmax + 1}], 
{i, 
                n - 1, 0, -1}]; 
fnew[[1, 1]]]];

In the case of an American option, the option value f() at (i, j) given (Subscript[M, i, j]) takes the maximum of the previous european value and the exercise price (Subscript[M, i, j]) - K.
While at maturity, the option value at node (N, j) given (Subscript[M, N, j]) is again the same payoff f((Subscript[S, N, j]) | (Subscript[M, N, j])) as described in the previous equation and this is described by the following code:

\(TraditionalForm\`f(S\_\(i, j\) | M\_\(i, j\)) = 
            max(\(e\^\(\(-r\)\ Δt\)\) 
                  p\ \(f(S\ \(d\^j\) u\^\(i + 1 - j\) | 
                      max(\ S\ d\^j\ u\^\(i + 1 - j\), 
                        M\_\(i, j\))\ )\) + \(\(e\^\(\(-r\)\ Δt\)\)(1 - p)\)\ \(f(S\ \(d\^\(j + 1\)\) u\^\(i - j\) | M\_\(i, j\))\), 
              M\_\(i, j\) - K\ )\)
f((Subscript[S, N, j]) | (Subscript[M, N, j])) = 
            max( (Subscript[M, N, j]) - K, 0)

All of this results in the following code:

CRRAppAmericanCallOnMax[S_, M0_, K_, r_, q_, σ_, T_, n_] := 
        Module[{dt, a, u, d, P, Q, TreeOp, value}, 
dt = 
            T/n; 
a = 
            Exp[(r - q)*dt]; 
u = 
            Exp[σ*Sqrt[dt]]; 
d = 
            Exp[-σ*
                Sqrt[dt]]; 
P = (a - d)/(u - 
                  d); 
Q = 1 - P; 
P = 
            P*Exp[-r*dt]; 
Q = 
            Q*Exp[-r*dt]; 
TreeOp[node_, level_, 
              M_] := TreeOp[node, level, M] = 
              If[level == n, Max[M - K, 0], 
                Max[{P, 
                      Q} . {TreeOp[node, level + 1, 
                        Max[S*d^node*u^(level + 1 - node), M]], 
                      TreeOp[node + 1, level + 1, M]}, 
                  Max[M - K, 0]]]; 
value = 
            TreeOp[0, 0, M0]; 
Clear[
            TreeOp]; 
value];

CRRAmericanLookbackPut[S_, M0_, r_, q_, σ_, T_, n_]:=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, value, delta, gamma, 
              theta},
P = Pfunc[u, d, a, disc];
Q = Qfunc[u, d, a, disc];
TreeOp[maxS_,node_,level_]:=TreeOp[maxS, node, level]=If[level == n,Max[0, 
                        maxS - S d^node u^(level - node)],Max[{P, Q}.{TreeOp[Max[maxS, 
                                      S d^node u^(level + 1 - node)],node,level + 1],TreeOp[maxS, node + 1, 
                                    level + 1]},Max[0, 
                            maxS - S d^node u^(level - node)]]];
value=TreeOp[M0,0,0];
  
Clear[TreeOp];
value];

CRRAmericanLookbackPut[50, 50, 0.2, 0.0, 0.25, 1.0, 50]

CRR Binomial model with modified binomial coefficients

Derivation of Coefficients

A general formula for the paths in a lookback option can be demonstrated in a fashion similar to that for a path independent option.  First, consider a path independent binomial pricing model after n total steps have been taken of which i are downward steps.  The only paths after n-1 steps that contribute to the paths at n steps are those that have the same number or one fewer downward steps, so we have that

paths(n,i) = paths(n-1, i-1) + paths(n-1, i)

But the fundamental recursion relation among the binomial coefficients is exactly the same, i.e.

(({{n}, {i}}))==(({{n - 1}, {i - 1}}))+(({{n - 1}, {i}}))

as is proved in many elementary texts.  Since the binomial coefficient holds for the case n=i=0, we have a proof by induction that

paths(n, i)=(({{n}, {i}}))

A similar result can be derived for path dependent options by looking at the recursion relations in the number of paths passing through each node.  First we introduce an additional parameter m, which is the largest number of net upward price steps achieved so far along the path.  At each node, it can clearly be no larger than total number of upward steps taken along the path so far, n-i.  Likewise, it can be no smaller than the net number of upward steps taken to reach the current node, n-2i, which is equivalent to the restriction that the highest price reached along the path so far must be at least equal to the current price.  Finally, m starts at zero and can never decrease.

As with the path independent case, the only paths after n-1 steps that contribute to the paths at n steps are those that have the same number or one fewer downward steps.  In addition, m must either be the same or increase by one, and it can only have increased from the prior node if it is equal to its minimum value at the current node.  This leads to two recursion relations:

paths(n, i, m) = paths(n-1, i-1, m) + paths(n-1, i, m),        m > n-2i
paths(n, i, m) = paths(n-1, i, m-1) + paths(n-1, i, m),        m = n-2i

The first relation is exactly the same as in the path independent case, corresponding to those movements that do not change the maximum price to date.  The second case corresponds to those paths where the price at the current node is the highest price achieved so far.  There is no contribution from the (n - 1, i - 1) node since after a downward price movement, the current price can never be the highest price achieved so far.  There are now, however, two contributions from the (n - 1, i) node since the prior m could have either been the same (the second term, which is the same in both relations) or the prior m could have been one less, which in any case is the minimum value at that node since if m = n - 2  i, then (m - 1) = (n - 1) - 2  i.

Because the recursion relation is so similar to the simple binomial tree, the general formula for the number of paths in the path dependent case is easily shown to be

paths(n, i, m)=(({{n}, {i + m}}))-(({{n}, {i + m + 1}}))

This modified binomial relation is a generalisation of the usual binomial formula allowing for the inclusion of the path dependent information of the maximum price achieved so far.  Because the general formula only depends on i + m, the two recursion relations yield the same relationship after direct substitution:

(({{n}, {i + m}}))-(({{n}, {i + m + 1}}))=(({{n - 1}, {i + m - 1}}))-(({{n - 1}, {i + m - 1 + 1}}))+(({{n - 1}, {i + m}}))-(({{n - 1}, {i + m + 1}}))

This is seen to be true by the fundamental recursion relation among the binomial coefficients, the two positive terms on the r.h.s. summing to the positive term on the l.h.s. etc.  Since n = i = m = 0 yields one path, the formula is proved by induction.

Coding of Equations

A direct implementation of the above formula (12) including the restrictions on time level i and nodes j results in the following code:

CRRFastCallonMaxApp[S_, M0_, K_, r_, q_, σ_, T_, n_] := 
        Module[{u = Exp[σ Sqrt[T/n]], 
            d = Exp[-σ Sqrt[T/n]], a = N[Exp[(r - q) T/n]], 
            disc = N[Exp[-r T/n]], P, Q, value, i, Payoffs, pq, 
            A}, 
P = 
            N[disc*(a - d)/(u - d)]; 
Q = 
            N[disc - 
                disc*(a - d)/(u - d)]; 
Payoffs = 
            Table[Max[Max[M0, S u^(n - i)] - K, 0], {i, 0, 
                n}]; 
pq = 
            Table[P^(n - i) Q^i, {i, 0, n}]; 
A = 
            SparseArray[{{i_, i_} -> 
                  1, {i_, j_} /; And[j < 2  i, j > i] -> 
                  Binomial[n, n - (j - i)] - 
                    Binomial[n, n - (j - i) + 1]}, {n + 1, 
                n + 1}]; 
value = 
            pq . A . Payoffs; 
value];

Using formula (14) and the restrictions on i and m, the value of the call on the max can now be expressed as a double sum over i and m. This reduces the dimensionality of the problem and results in a faster binomial model to determine the same Lookback option.

CRRFastCallonMaxApp[S_, M0_, K_, r_, q_, σ_, T_, n_]:=
Module[{dt, a, u, d, P, Q, i, m},
dt = T/n;
a = Exp[(r - q)*dt];
u = Exp[σ*Sqrt[dt]];
d = Exp[-σ*Sqrt[dt]];
P = 
                Exp[-r*dt] (a - d)/(u - d);
Q = Exp[-r*dt] - P;
(∑_(i = 0)^n)(P^(n - i)) (Q^i) (∑_(m = Max(0, n - 2 i))^(n - 1))((({{n}, {i + m}}))-(({{n}, {i + m + 1}}))) Max(0, Max(M0, S (u^m)) - K)];

CRRFastCallonMaxApp[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRFastCallonMaxApp[50, 50, 50, 0.05, 0, 0.2, 1, 500] // 
      Timing

CRR Binomial model with polynomial structure

In order to establish all the possible binomial paths that can be generated in n steps, we first construct a list with k down jumps and n-k up jumps and then form all such possible list for all values of k from 0 to n with the command Table[Join[Table[d, {k}], Table[u, {n - k}]], {k, 0, n}]. To list all the possible permutations for each node in a binomial tree we map the function Permutations onto the previously constructed array. Each sublist in this array is then mapped onto with the function FoldList[Times,S,#]& to create the list of all possible asset paths out to time n.

BinPaths[S_, u_, d_, n_Integer] := 
      Map[FoldList[Times, S, #] &, (Permutations  /@  
            Table[Join[Table[d, {k}], Table[u, {n - k}]], {k, 0, 
                n}]), {2}]

BinPaths[S, u, d, 4]

To replace all possible asset paths between any two points by the maximum of those paths we merely replace the function FoldList[Times, S, #] & by Max[FoldList[Times, S, #]] & in the code above

CRRAppEuroTheorMaximumPaths[S_, u_, n_] := 
      Assuming[{{u, S} ∈ Reals , S > 0,  u > 1},  
        Simplify[
          Map[Max[FoldList[Times, S, #]] &, (Permutations  /@  
                Table[Join[Table[1/u, {k}], Table[u, {n - k}]], {k, 0, 
                    n}]), {2}]]]

CRRAppEuroTheorMaximumPaths[S, u, 4]

To determine the expected maximum price at each final node we map Mean onto each collection of path maxima

CRRAppEuroTheorPathMaxima[S_, u_, n_] := 
      Assuming[{{u, S} ∈ Reals , S > 0,  u > 1},  
        Simplify[
          Mean  /@  Map[
              Max[FoldList[Times, S, #]] &, (Permutations  /@  
                  Table[Join[Table[1/u, {k}], Table[u, {n - k}]], {k, 0, 
                      n}]), {2}]]]

CRRAppEuroTheorPathMaxima[S, u, 4]

Finally to calculate the expected Call on Max payoff at each final node we merely replace the function FoldList[Times, S, #] & in the code above by Max[0, 
          Max[Prepend[FoldList[Times, S, #], M0]] - K] &. In this way we account for the prior maximum M0 in every path to each final node at maturity.

CRRAppEuroTheorCallOnMaxPayoffs[S_, M0_, K_, u_, n_] := 
      Assuming[{{u, S} ∈ Reals , S > 0,  u > 1},  
        Simplify[
          Mean  /@  Map[
              Max[0, Max[Prepend[FoldList[Times, S, #], M0]] - 
                    K] &, (Permutations  /@  
                  Table[Join[Table[1/u, {k}], Table[u, {n - k}]], {k, 0, 
                      n}]), {2}]]]

CRRAppEuroTheorCallOnMaxPayoffs[S, M0, K, u, 4]

While the above code correctly determines the maximum asset value and Call on Max payoffs at each final node in a binomial tree, it is slow because we have to calculate all the possible n! permutations of up and down jumps up to time n. However we can also use Mathematica's symbolic manipulation ability to inspect the algebraic structure of each maximum value in an equity path from initiation to maturity. This analysis shows that they are polynomials whose coefficients are recursively determined by earlier values and hence which allow for a more efficient description and faster evaluation. To determine the structure of these polynomials consider the output of the function CRRAppEuroTheorPathMaxima for different values of n, the number of steps to maturity.

CRRAppEuroTheorPathMaxima[S, u, 8]

CRRAppEuroTheorPathMaxima[S, u, 11]

Inspection of the above output can be further aided by using CoefficientList with different values of n

With[{n = 8}, CoefficientList[#, u] &  /@  
        CRRAppEuroTheorPathMaxima[S, u, n]]

Let (Subscript[C, k - i]) represent the ith coefficient of (u^((n - 2  k) + i)) in the kth polynomial (Subscript[P, k]) which yields the maximum asset price in the path to the kth node from the top at the end of a binomial tree, where k = 0 represents the top node. Then inspection of the above output shows that the equations describing the coefficients are given by:

(Subscript[C, k]) = (Subscript[C, n - k]),⌈n/2⌉<k≤n
(Subscript[C, k])=(({{n}, {k}}))-(∑_(i = 0)^(k - 1))(Subscript[C, i]),k ≤ ⌈n/
                              2⌉
(Subscript[P, k]) = (Subscript[P, n - k])/(u^((2  k - n))), ⌊n/2⌋<k≤n
\(P\_k = 
                        u\^\((n - 2  k)\)\ \(∑\+\(i = 0\)\%k C\_\(k - i\)\ u\^i\), k ≤ \ ⌊n/2⌋\)

Furthermore the coefficients are symmetric and the resultant equation describing the path maximum up to the kth node from the top of a binomial tree after n steps after starting from an asset price of S is

(Subscript[M, k]^n) =(Piecewise[{{S(∑_(i = 0)^k)(((u^((n - 2  k))) (Subscript[C, k - 
                                      i]) (u^i))/((({{n}, {k}})))), \(k ≤ \ ⌊n/2⌋\)}, {S(∑_(i = 0)^(n - k))(((u^((n - 2  k))) 
                                (Subscript[C, n - k - i]) (u^i))/((({{n}, {k}})))),  ⌊n/
                              2⌋ < k ≤ n}}])

The above equations are contained in the following code for the function CRRAppFastEuroTheorPathMaxima

CRRAppFastEuroTheorPathMaxima[S_, u_, n_] := 
      Module[{poly, coefft, maxs}, 
coefft[0] = 1; 
        coefft[k_] /; (k <= Ceiling[n/2]) := coefft[k] = 
            Binomial[n, k] - 
              Sum[coefft[i], {i, 0, k - 1}]; 
coefft[
              k_] /; (Ceiling[n/2] < k <= n) := coefft[k] = 
            coefft[n - k]; 
poly[
              k_] /; (k <= Floor[n/2]) = poly[k] = 
            u^(n - 2  k) Sum[coefft[k - i] u^i, {i, 0, k}]; 
        poly[k_] /; (Floor[n/2] < k <= n) = poly[k] = 
            poly[n - k]/u^(2  k - n); 
maxs = 
          S Table[poly[k]/Binomial[n, k], {k, 0, n}]; Clear[poly, coefft]; 
        maxs]

Running the code for this polynomial based function we see that it results in the same outcome as the slower but more straightforward version.

CRRAppFastEuroTheorPathMaxima[S, u, 11]

The code for the more efficient version of the path maxima can then be used to determine a faster version for the Call on Max payoffs. It merely requires replacing all instances of S (u^k) by Max[0, Max[M0, S (u^k)] - K], which is the call payoff that accounts for earlier maxima M0.

CRRAppFastEuroTheorCallOnMaxPayoffs[S_, M0_, K_, u_, n_] := 
      Module[{maxs}, 
maxs  = 
          CRRAppFastEuroTheorPathMaxima[S, u, 
            n]; 
Together[ 
            Expand[maxs]] /. {S u^k_.  -> Max[0, Max[M0, S u^k] - K], 
            S -> Max[0, Max[M0, S ] - K]}]

CRRAppFastEuroTheorCallOnMaxPayoffs[S, M0, K, u, 4] == 
      CRRAppEuroTheorCallOnMaxPayoffs[S, M0, K, u, 4]

This faster code is the same as the previous permutaions code but much faster.

Both CRRAppFastEuroTheorPathMaxima and CRRAppFastEuroTheorCallOnMaxPayoffs can be adapted to account for volatility σ and maturity T rather than the usually unknown upward factor u by observing that in the CRR model u = (e^(σ (Sqrt[T/n])))

CRRAppEuroPathMaxima[S_, σ_, T_, n_] := 
      Module[{dt, u}, 
dt = T/n; 
u = 
          Exp[σ*
              Sqrt[dt]]; 
CRRAppFastEuroTheorPathMaxima[S, 
          u, n]]

CRRAppEuroCallOnMaxPayoffs[S_, M0_, K_, σ_, T_, n_] := 
      Module[{dt, u, S1, M1, K1, u1, n1}, 
dt = 
          T/n; 
u = 
          Exp[σ*
              Sqrt[dt]]; 
CRRAppFastEuroTheorCallOnMaxPayoffs[S1, M1, K1, u1, n] /. 
          Thread[{S1, M1, K1, u1} -> {S, M0, K, u}]]

CRRAppEuroCallOnMaxPayoffs[S, M0, K, σ, T, 4]

Finally CRRAppFastEuroTheorCallOnMaxPayoffs can be used to calculate the Call on Max option value

CRRAppFastEuroCallOnMax[S_, M0_, K_, r_, q_, σ_, T_, 
          n_] := 
Module[{dt, a, u, d, P, Q,  PQ,  
            payoffs}, 
dt = T/n; 
a = 
            Exp[(r - q)*dt]; 
u = 
            Exp[σ*Sqrt[dt]]; 
d = 
            Exp[-σ*Sqrt[dt]]; 
P = 
            Exp[-r*dt] (a - d)/(u - d); 
Q = 
            Exp[-r*dt] - P; 
PQ = 
            Table[Binomial[n, j] Q^j P^(n - j), {j, 0, 
                n}]; 
payoffs = 
            CRRAppFastEuroTheorCallOnMaxPayoffs[S1, M1, K1, u1, n] /. 
              Thread[{S1, M1, K1, u1} -> {S, M0, K, 
                    u}]; 
PQ . payoffs];

CRRAppFastEuroCallOnMax[S, M0, K, r, q, σ, T, 4]

Comparing the results we see that the polynomial structure and adjusted binomial functions are respectively 34 and 110 times as fast as the usual recursive binomial model

CRRAppEuropeanCallOnMax[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRAppFastEuroCallOnMax[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRFastCallonMaxApp[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

Relation to Other Lookback Options

Other Lookback options can be readily valued merely by replacing the payoff function in the original summation, and possibly interchanging P/Q and u/d if the minimum price along the path is required.  However, to use the more efficient forms of computation for the floating lookbacks, the sums must be evaluated somewhat differently.  The floating lookback put is the simplest example as it has a payoff equal to the excess of the maximum price achieved so far over the current price.  In its simplest form, it can be computed as

CRRLookbackPutApp[S_, M0_, r_, q_, σ_, T_, n_]:=
Module[{dt, a, u, d, P, Q, i, m},
dt = T/n;
a = Exp[(r - q)*dt];
u = Exp[σ*Sqrt[dt]];
d = Exp[-σ*Sqrt[dt]];
P = 
                Exp[-r*dt] (a - d)/(u - d);
Q = Exp[-r*dt] - P;
(∑_(i = 0)^n)(P^(n - i)) (Q^i) (∑_(m = Max(0, n - 2 i))^(n - i))((({{n}, {i + m}}))-(({{n}, {i + m + 1}}))) (Max(M0, S (u^m)) - 
                            S (d^i) (u^(n - i)))];

CRRLookbackPutApp[50, 50, 0.05, 0, 0.2, 1, 100] // Timing

It might seem that the more efficient calculation of the sums is not workable since the payoff now depends on both i and m, but that is not the case.  As with the Call on the Max, the payoff can be split into two terms:  (i) the maximum price along the path and (ii) the expected present value of the stock price.  This is almost identical to the call on the max case, except the strike is replaced by the stock price.

CRRLookbackPutAppComp := 
        Compile[{{S, _Real, 0}, {M0, _Real, 0}, {r, _Real, 0}, {q, _Real, 
              0}, {σ, _Real, 0}, {T, _Real, 0}, {n, _Real, 
              0}}, 
Module[{dt, a, u, d, P, Q, j, m, PQ, 
              Payoffs, InnerSums, PQfactors, logPQ, 
              logM0}, 
dt = T/n; 
a = 
              Exp[(r - q)*dt]; 
u = 
              Exp[σ*Sqrt[dt]]; 
d = 
              Exp[-σ*Sqrt[dt]]; 
P = 
              Exp[-r*dt] (a - d)/(u - d); 
Q = 
              Exp[-r*dt] - P; 
logPQ = 
              Log[P/Q]; 
logM0 = 
              Log[M0]; 
PQfactors = 
              Table[logPQ + 
                  Log[(2  j - n - 
                          1) (j + 
                            1)/((2  j - n + 1) (n - j + 
                                1))], {j, 
                  n, (n + Mod[n, 2])/2 + 
                    1, -1}]; 
PQ = 
              Exp[Reverse[
                  FoldList[Plus, n Log[Q], 
                    PQfactors]]]; 
Payoffs = 
              Table[Exp[m*logPQ + Max[logM0,  Log[S] + m Log[u]]],  {m, 0, 
                  n}]; 
InnerSums = 
              If[Mod[n, 2] == 
                  0, 
FoldList[(#1 + Total[#2]) &, 
                  Payoffs[[1]], Partition[Rest[Payoffs], 2]], 
                Rest[FoldList[(#1 + Total[#2]) &, 0, 
                    Partition[Payoffs, 
                      2]]]
]; 
PQ . 
                InnerSums - 
              S Exp[-q T]
]];

CRRLookbackPutAppComp[50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRLookbackPutAppComp[50, 50, 0.05, 0, 0.2, 1, 500] // 
      Timing

CompFastLookbackPut=Compile[{{S, _Real, 0}, {M0, _Real, 0}, {r, _Real, 0}, {q, _Real, 
                0}, {σ, _Real, 0}, {T, _Real, 0}, {n, _Integer, 0}},Module[{u = Exp[σ Sqrt[T/n]],d = Exp[-σ Sqrt[T/n]],a = Exp[(r - q) T/n],disc = Exp[-r T/n],P = 0. ,Q = 0. ,value = 0. ,pmax = 0,node = 0,payoffS,payoffMaxS,pq,pqB,A=Table[If[n - node  ≥  
                                pmax  ≥  n - 2  node, 
                              Binomial[N[n], pmax + node] - 
                                Binomial[N[n], pmax + node + 1], 0],{node, 0, n},{pmax, 0, n}]},
P = disc (a - d)/(u - d);
Q = disc - P;
payoffMaxS = 
                    Table[Max[M0, S u^pmax], {pmax, 0, n}];
payoffS = 
                      Table[S d^node u^(n - node), {node, 0, n}];
pq = 
                    Table[P^(n - node) Q^node, {node, 0, n}];
pqB = 
                      Table[Binomial[N[n], node] 
                          P^(n - node) Q^node, {node, 0, n}];
value=pq . A . payoffMaxS-pqB . payoffS]];

CompFastLookbackPut[50, 50, 0.05, 0, 0.2, 1, 500] // Timing

Speed up of new binomial model

Although much faster than a recursive tree and having a computation time that grows only as (n^2), this is by no means the fastest arrangement for valuing such a sum.  The repeated calculation of the payoffs and the binomial coefficients is a considerable drag on the program, which can most naturally be remedied by replacing the summation variable i with j=i+m.

CRRFastCallonMaxApp2[S_, M0_, K_, r_, q_, σ_, T_, 
          n_]:=
Module[{dt, a, u, d, P, Q, j, m},
dt = T/n;
a = Exp[(r - q)*dt];
u = Exp[σ*Sqrt[dt]];
d = Exp[-σ*Sqrt[dt]];
P = 
                Exp[-r*dt] (a - d)/(u - d);
Q = Exp[-r*dt] - P;
(∑_(j=(1/2) (n+n mod 2))^n)(P^(n - j)) (Q^j) ((({{n}, {j}}))-(({{n}, {j + 1}}))) (∑_(m = 0)^(2 j - n))(P^m) (Q^(-m)) Max(0, Max(M0, S (u^m)) - K)];

CRRFastCallonMaxApp2[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRFastCallonMaxApp2[50, 50, 50, 0.05, 0, 0.2, 1, 500] // 
      Timing

At this point, the calculation has been sped up considerably, but the number of computations still scales as (n^2) because of the double sum.   The inner sum, however, depends only on j through the upper limit on m and the dependence is monotonic, allowing the inner sums to be easily generated recursively.  The resulting vector can then be dotted with a vector of j-dependent terms to generate the final value, effectively eliminating the double sum and yielding a computation time that grows only as n.

CRRFastCallonMaxApp3[S_, M0_, K_, r_, q_, σ_, T_, 
          n_] := 
Module[{dt, a, u, d, P, Q, j, m, PQ, 
            Payoffs, InnerSums}, 
dt = 
            T/n; 
a = 
            Exp[(r - q)*dt]; 
u = 
            Exp[σ*Sqrt[dt]]; 
d = 
            Exp[-σ*Sqrt[dt]]; 
P = 
            Exp[-r*dt] (a - d)/(u - d); 
Q = 
            Exp[-r*dt] - P; 
PQ = 
            Table[P^(n  -  j)*
                Q^j*(Binomial[n,  j]  - 
                    Binomial[n,  j  +  1]), {j,  (n + Mod[n, 2])/2, 
                n}]; 
Payoffs = 
            Table[(P/Q)^m Max[0,  Max[M0,  S*u^m]  -  K],  {m, 0, 
                n}]; 
InnerSums = 
            If[Mod[n, 2] == 
                0, 
FoldList[(#1 + Total[#2]) &, 
                Payoffs[[1]], Partition[Rest[Payoffs], 2]], 
              Rest[FoldList[(#1 + Total[#2]) &, 0, 
                  Partition[Payoffs, 
                    2]]]
]; 
PQ . 
            InnerSums
];

CRRFastCallonMaxApp3[50, 50, 50, 0.05, 0, 0.2, 1, 100] // 
      Timing

CRRFastCallonMaxApp3[50, 50, 50, 0.05, 0, 0.2, 1, 500] // 
      Timing

To fully realize this computation time scaling, it is also necessary to streamline the calculation of the j-dependent terms.  Direct computation of the binomial coefficients is much slower than computing the j-dependent terms recursively, which can be done by computing the ratio of consecutive terms, and then generating them by succesive multiplications starting with j = n.  Finally, by compiling it's possible to get another significant reduction in computation time, but in order for the compiled code to work for large n, the calculation must avoid intermediate results that are larger than machine precision reals.  

Although the adaptations for compiling somewhat obscure the form of the calculation, one is worth noting.  For the compiled form, it is somewhat simpler to evaluate the payoff function  Max[0, Max[M0, S*u^m] - K] as the equivalent  Max[K, M0, S*u^m] - K.  This avoids having to take the log of a sum, but effectively splits the calculation into two pieces:  (i) the expected value of the greater of the strike or the maximum price along the path and (ii) the expected present value of the strike.  This division of the calculation is the natural route to extend the method to standard lookbacks. These methods can be significantly enhanced with compilation to give speeds comparable to the continuous analytical results.

CRRFastCallonMaxAppComp := 
        Compile[{{S, _Real, 0}, {M0, _Real, 0}, {K, _Real, 0}, {r, _Real, 
              0}, {q, _Real, 0}, {σ, _Real, 0}, {T, _Real, 
              0}, {n, _Real, 0}}, 
Module[{dt, a, u, d, P, 
              Q, j, m, PQ, Payoffs, InnerSums, PQfactors, logPQ, 
              logKM0}, 
dt = T/n; 
a = 
              Exp[(r - q)*dt]; 
u = 
              Exp[σ*Sqrt[dt]]; 
d = 
              Exp[-σ*Sqrt[dt]]; 
P = 
              Exp[-r*dt] (a - d)/(u - d); 
Q = 
              Exp[-r*dt] - P; 
logPQ = 
              Log[P/Q]; 
logKM0 = 
              Max[Log[K], Log[M0]]; 
PQfactors = 
              Table[logPQ + 
                  Log[(2  j - n - 
                          1) (j + 
                            1)/((2  j - n + 1) (n - j + 
                                1))], {j, 
                  n, (n + Mod[n, 2])/2 + 
                    1, -1}]; 
PQ = 
              Exp[Reverse[
                  FoldList[Plus, n Log[Q], 
                    PQfactors]]]; 
Payoffs = 
              Table[Exp[m*logPQ + Max[logKM0,  Log[S] + m Log[u]]],  {m, 
                  0, n}]; 
InnerSums = 
              If[Mod[n, 2] == 
                  0, 
FoldList[(#1 + Total[#2]) &, 
                  Payoffs[[1]], Partition[Rest[Payoffs], 2]], 
                Rest[FoldList[(#1 + Total[#2]) &, 0, 
                    Partition[Payoffs, 
                      2]]]
]; 
PQ . 
                InnerSums - 
              K Exp[-r*T]
]];

CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, 1, 500] // 
      Timing

CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, 1, 5000] // 
      Timing

CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, 1, 50000] // 
      Timing

CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, 1, 500000] // 
      Timing

Comparison to prior results

Although there are points of interest in the comparison, continuously valued options for which analytic solutions are available do not reflect the reality that the market is rarely continuously priced.  Comparing the analytic result to the discrete result we find, as others have, a very slow convergence of the discrete to the continuous result.  Using the Conze and Viswanathan continuous formula, we find as did Babbs that the difference between option values for daily and continuous sampling is essentially independent of time to expiry.

Plot[{CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, j/365, 
              Round[j]] - CVCallOnMax[50, 50, 50, 0.05, 0, 0.2, j/365]}, {j, 
          1, 365}, 
        AxesLabel -> {"Days to Expiry",  "Option Value Difference\
"}];

The explanation for this "striking feature" Babbs found is hinted at by the fact that the difference remains constant as the sampling rate remains constant.  When the time to expiry is fixed and the number of time steps is varied, we find that the error is directly proportional to the volatility during one time step.  Since this volatility is directly proportional to the price difference between nodes, it indicates a systematic problem with the valuation of the maximum price in the continuous limit that has not been accounted for.

ListPlot[
        Table[{0.2/Sqrt[j], 
            CRRFastCallonMaxAppComp[50, 50, 50, 0.05, 0, 0.2, 1, j] - 
              CVCallOnMax[50, 50, 50, 0.05, 0, 0.2, 1]}, {j, 10, 1000, 10}], 
        PlotJoined -> True, 
        AxesLabel -> {"One Time Step Volatility",  "Option Value \
Difference"}];

When valuing the price of the stock at a particular node, whether for a vanilla option or a lookback option, the stock price is taken to be S(u^i).  This represents the perfectly reasonable assumption that if we sample the price at that instant and near (i.e. represented by) that node, it's most likely to be the average value.  This leads to oscillations in price about the continuous value when binomial trees are used to value vanilla options (see Shaw, chapter 19 for some examples).  The cause of the oscillations is that the payoff of the option at the tree node nearest the strike price oscillates rapidly as the tree size is changed, unless the tree is of a "magic" size where the strike falls roughly half way in between two nodes.  While the payoff based on the maximum price largely smooths out the oscillation problem, using this standard value fails to account for the fact that we are trying to compare to the continuous limit.

Conclusion

We have compared a number of different methods for the evaluation of Lookback options. The commonly used approach is due to Conze and Viswanathan (1991) continuous formulae which are derived from stochastic theorems about the distributions of processes and their maxima and minima. It has since been observed by Babbs (2000) amongst others that such continuous results assume an infinite sampling rate which far exceeds the actual determination of the underlying asset prices, usually once every day. This has resulted in an overpricing of the observed values. To resolve this difficulty Babbs suggested the introduction of a reflecting barrier which corresponds to the application of a change of measure. This approach introduces a continuity correction error discussed in the last section above. However two alternate approaches using the Cox, Ross and Rubinstein (1979) binomial model are capable of obtaining the correct discrete results without needing to introduce arbitrary barriers. One approach resolves the algebraic structure of the discrete payoffs at each node in the binomial tree, while the other finds a simple extension of the well known binomial formula to the case of the path dependent Lookback option thereby reducing the dimensionality of the problem to the same level as that of a Vanilla option. It is shown that both of these new formulae yield the same results and converge to the Babbs pricing formula. Lastly these different approaches can be sped up enormously by further utilising the extended binomial relation to reduce multiple summations as well as by the method of compilation.

References

Babbs,S., Binomial Valuation of Lookback Options, Journal of Economic Dynamics and Control, 24,1499-1525, 2000.
Benninga,S. and Wiener,Z.,The Binomial Option Pricing Model, Mathematica in Education and Research, 6 (3), 27-34, Summer 1997.
Benninga,S. and Wiener,Z., Binomial Option Pricing, the Black-Scholes Option Pricing Formula and Exotic Options, Mathematica in Education and Research, 6 (4), 11-14, Fall 1997.
Black,F. and Scholes,M., The Pricing of Options and Corporate Liabilities, Journal of Political Economy, 81, 659-673, 1973.
Conze, A. and Viswanathan,R., "Path Dependent Options: The Case of Lookback Options", Journal of Finance,Vol. 46, No 5, pgs 1893-1907, 1991.
Cox,J., Ross,S. and Rubinstein,M., Option Pricing: A Simplified Approach, Journal of Financial Economics, 7, p.229, 1979.
Cox,J. and Ross,S., The Valuation of Options for Alternative Stochastic Processes, Journal of Financial Economics, 3, 145-166, 1976.
Goldman, M.B., Sosin, H.B., Gatto, M.A., 1979a, "Path Dependent Options: buy at the low, sell at the high", Vol. 34, pgs 1111-1127.
Harrison, J.M., Brownian Motion and Stochastic Flow Systems, Wiley New York, 1985.
Hull,J., Options, Futures and Othe Derivatives, Prentice Hall, 6th Edn 2005.
Jarrow,R. and Rudd,A., Option Pricing, Irwin, 1983.
Kelly, M., Evaluation Of Financial Options using Radial Basis Functions in Mathematica, submitted to The Mathematica Journal, June 2006.
MathConsult. UnRisk Pricing Engine. http://www.unriskderivatives.com/
Shaw,W., Modelling Financial Derivatives with Mathematica, Cambridge University Press, 1998.
Trott,M., The Mathematica Guidebook for Numerics, Springer, 2006.