Evaluating Financial Options Using Continuous Cellular Automata

Dr Michael Kelly

Stuart Graduate School of Business, IIT

Introduction

The traditional use of cellular automata in the study of scientific topics has been to specify simple rules and see whether the resulting CA structures behave like some known real world structure that the scientist wishes to model. To this end Wolfram specified a simple model for financial market price fluctuations in Chapter 8 of NKS. While exhibiting the characteristic fractal property of market prices this model did not allow an analysis of actual financial instruments. The evaluation of such instruments, especially financial derivatives, constitutes the main task of financial quants employed by banks and financial institutions. What this paper aims to do is show that the appropriate modification of the continuous CA transformation rules makes it possible to implement efficient code for the evaluation of financial options by modelling the trinomial lattice structures which correspond to the partial differential equations which describe all option behaviour. We show that there a number of possible models using inbuilt functions NestList, ListCorrelate and CellularAutomaton which also evaluate options by mimicking the stochastic evolution of underlying stock prices as cellular automata on a trinomial grid.

Basic Ideas

Consider the problem of trying to determine the path of a stock option as it evolves through time. Amongst the simplest models used by financiers is a trinomial structure which hypothesises that the next value of the stock can be in one of three possible future states: either up, down or at the same level. This is the so called trinomial model which envisages three possible future scenarios for each current state. This can be represented diagramtically by Figure 1 which has the horizontal axis representing time and the vertical axis for the equity under consideration. This three state model can be thought of as a continuous cellular automaton moving horizontally through time, where each cell represents the current state of the stock, equity or option and future states are obtained in the next step from the previous state. The up, down and same levels can be obtained from the prior or initial state by shifting, that is rotating to the left the initial state either zero, one or two steps. The transition probabilities to the up, down and same levels are given by p, q and r respectively so that a one time step trinomial operation is described by the following function:

Off[General::"spell", General::"spell1"]

\(\(OneTimeStepEuroTri\  = \ 
        Compile[{{initial, \ _Real, \ 1}, {p, \ _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {r, \ _Real, \ 0}}, 
	
          p*Drop[initial, \(-2\)]\  + \ 
            q*Drop[RotateLeft[initial, 2], \ \(-2\)] + 
            r*Drop[RotateLeft[initial], \ \(-2\)]];\)\)

The combined operation of left shifting and multiplying by transition probabilities constitute a function which describes the action of a cellular automaton. What is not explicit in such trinomial models is the amount by which the stock goes up, represented by the function up[.], the amount by which it goes down, represented by the function down[.] and finally the amount by which the stock remains at the same level, represented by the function stay[.]. All of these functions are dependent upon market parameters which have to be specified initially. These include the current value of the stock, represented by the symbol S, the strike or exercise price K, which is the amount for which the option writer agrees to buy or sell the stock at some future date, called the maturity, T. Other parameters include the current interest rate r, the stock dividend yield q and the annualised standard deviation of the stock or its volatility σ. Finally since this is a discrete process we also have the number of times the time increment dt divides into the time interval T/dt = n, which constitute the number of time steps in the model.

Neighbouring Symmetry Constraints

The problem with all such models is in the specification of the amounts up[.], down[.] and stay[.] by which the stochastic process representing the equity deviates at ech incremental time step. There are two generally accepted models for this named after the people who first specified them, namely the Jarrow-Rudd and Cox-Ross-Rubinstein models. Both models are determined by equating the  values for the mean and variance of the continuous equity process with its corresponding first two moments in the discrete trinomial representation. The problem with such an approach is that there are three parameters to specify, but only two equations. This leaves one of the parameters as a variable, which accounts for the difference between the models. A better approach is described by William Shaw in his book "Modelling Financial Derivatives with Mathematica", chapters 18 to 20. Starting with the general power solution to the Black-Scholes model

\(TraditionalForm\`V = 
      Z\ \((S\/K)\)\^n\ \(e\^\(\(-r\)\ \((T - 
                  t)\)\)\) e\^\(\((\((r - q)\)\ n + \(\(n(n - 1)\)\ σ\^2\)\/2)\)\ \((T - t)\)\)\)

It can be shown that for appropriate choices of the parameters a and b, the relation between the up (u), down (d) and same level (s) amounts, along with the up, down and remain the same probabilities (p, q and r) are linked with the number of time steps n in the following relation:

p u^n  +  q d^n  +  r s^n  =  
      a^n b^((n - 1) n/2)

We first need to define appropriate functions for all of the above parameters before checking that the neighbouring symmetries apply.

Discount[n_, r_, T_] := N[Exp[-r T/n]];
afunc[n_, r_, q_, T_] := N[Exp[(r-q) T/n]];
bfunc[n_, σ_, T_] := N[Exp[σ^2 T/n]];

U[ b_]  :=  
        1/2 (b*(b + 1) + 
              Sqrt[b^4  +  2 b^3  +  b^2  -  4]);
up[a_,  b_]  :=  a U[b]
down[a_,  b_]  :=  a /U[b];
stay[a_,  b_]  :=  a;
pfunc[
          b_]  :=  (b - 1)/(1 - U[b])/(1/U[b]  -  
              U[b]);
qfunc[
          b_]  :=  (b - 1)/(1 - 1/U[b])/(U[b]  -  
              1/U[b]);
rfunc[b_]  :=  1  -  pfunc[b] - qfunc[b];

We will now define a function to check the supersymmetric structures defined by equation (2):

CheckSymmetry[u_,  d_, s_,   p_,  q_,  r_,  n_,  a_,  
        b_]  :=  Factor[
        Simplify[
          p u^n  +  q d^n  +  r s^n -  
            a^n b^((n - 1) n/2)]]

Table[
      CheckSymmetry[up[a, b], down[a, b], stay[a, b], pfunc[b], qfunc[b], 
        rfunc[b], 3, a, b], {n, -10, 10}]

The above clearly shows that supersymmetric identity of equation (2) is satisfied for a wide range of values of n

Basic European Options

Perhaps the simplest type of options are European options which can only be exercised at the maturity date T. This fixed date for the cashing in of the option allows a description of the probability distribution of the option from that of the underlying asset it is based upon. This was the essential insight of Black and Scholes and it can be shown that the trinomial models converge in value to that of the theoretical Black-Scholes mode for sufficiently large n. 
We will first give the trinomial code for the european call and put options and then show that we can also replicate these values from an appropriately defined continuous cellular automata.

\(\(OneTimeStepEuroTri\  = \ 
        Compile[{{initial, \ _Real, \ 1}, {p, \ _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {r, \ _Real, \ 0}}, 
	
          p*Drop[initial, \(-2\)]\  + \ 
            q*Drop[RotateLeft[initial, 2], \ \(-2\)] + 
            r*Drop[RotateLeft[initial], \ \(-2\)]];\)\)

\(\(ManyTimeStepEuroTri\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, \ 
              0}, \ {r, \ _Real, \ 0}, {n, \ _Integer, \ 0}}, 

          Module[{new\  = \ initial, \ k = 1, \ old = initial}, 

            For[k\  = \ 1, \ 
              k <= n, \ \(k++\), \ 
\((old\  = \ new; \ 
                new\  = \ OneTimeStepEuroTri[old, p, q, r])\)]; 

            Join[new, old]]];\)\)

EuroTrinomialCall[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[S s^node u^(n-node)-K,0], {node, 0,n, 1}],
		 Table[Max[S d^node s^(n-node)-K,0], {node, 1,n, 1}]];
result = ManyTimeStepEuroTri[payoffin, P, Q, R, n];
value = result[[1]];
delta = (result[[2]] - result[[4]])/((u-d)*S);
gamma = ((result[[2]] - result[[3]])/(u - s) -
(result[[3]] - result[[4]])/(s - d))*2/(S^2 (u-d));
theta = (result[[3]]-(s-1)*S*delta- result[[1]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

EuroTrinomialPut[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[K-S s^node u^(n-node),0], {node, 0,n, 1}],
		 Table[Max[K-S d^node s^(n-node),0], {node, 1,n, 1}]];
result = ManyTimeStepEuroTri[payoffin, P, Q, R, n];
value = result[[1]];
delta = (result[[2]] - result[[4]])/((u-d)*S);
gamma = ((result[[2]] - result[[3]])/(u - s) -
(result[[3]] - result[[4]])/(s - d))*2/(S^2 (u-d));
theta = (result[[3]]-(s-1)*S*delta- result[[1]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

The points to notice about the above code are that it is decidely procedural and that both option functions call the same compiled general procedure ManyTimeStepEuroTri which in turn repeatedly calls the single step function OneTimeStepEuroTri. If we first load in a package designed to calculate the traditional known theoretical values we can then see how well the trinomial tree schemes compare.

Needs["Derivatives`BlackScholes`"];

In the following calculations we determine the prices of european options and their sensitivity functions, known as the greeks, for the following market values: The stock price is S = $9, the strike price is K = $10, the risk free interst rate r = 10%, the dividend yield q = 3%, volatility σ = 20% and the time remaining tothe expiry date is T = 1 year. Here we choose n = 500 time steps.

\({BlackScholesCall[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], \ 
      BlackScholesCallDelta[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], 
      BlackScholesCallGamma[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], 
	
      BlackScholesCallTheta[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesCallRho[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesCallVega[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1]}\)

EuroTrinomialCall[9, 10,  500,  0.1,  0.03,  0.2,  1]

\({BlackScholesPut[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], \ 
      BlackScholesPutDelta[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], 
      BlackScholesPutGamma[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 1], 
	
      BlackScholesPutTheta[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesPutRho[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesPutVega[9, \ 10, \ 0.2, \ 0.1, \ 0.03, \ 
        1]}\)

EuroTrinomialPut[9, 10,  500,  0.1,  0.03,  0.2,  1]

Clearly the trinomial model is giving accurate results up to the fourth decimal place.

Continuous CAs for European Options

The typical structure of a continuous cellular automaton is described in chapter 3 of Wolfram's book "A New Kind Of Science" and consists of two main components: an evolutionary Step involving the mapping of some function onto the cells of an initial condition and another evolutionary List function which repeatedly applies the first function using the built-in command NestList. We shall also need a few supplementary functions padLeftRight and paddedPayoffs, whose main purpose is to fill out the tree structure into a rectangular array suitable for repeated application of the EvolveStep procedure.

padLeftRight[list_List, n_Integer] := 
        PadLeft[list, 
          n, {0}, (n - Length[list])/2];
paddedPayoffs[payoff_List, 
        payoffs_List] := padLeftRight[#, Length[payoff]] &  /@  
        payoffs

There are a number of ways we can construct the appropriate cellular automata. Either we can specify the structure of the Step and Evolve List functions ourselves, or we can utilise the inbuilt CellularAutomaton[.] function. In the first case we will replicate the behaviour of the OneTimeStepEuroTri function and build a compiled procedure that shifts the previous array of cells depending upon whether we have moved into the up, down states or remained in the same state. Suppose that we are in the stae called initial, then at the next step we will be at the same price level with probability r we will be at the level down with probability q or finally we will be at the upper level with probability p. If we are at the same price level then the prices remain unchanged, whereas if we go down one level then the second value now occupies the leading position and is therefore shifted to the left, but if we go one level up then the first element is forced into the second position. In this way the expected future value of the next step of the process is given by:  p*RotateRight[initial] +r*initial+q*RotateLeft[initial]. In order that first and last values at either end do not get incorrectly transposed, we remove them successively from each step by applying the procedure Take[.., {2,-2}]. Finally we repeatedly apply the single step procedure using NestList

CompCAEuropeanOptionEvolveStep  = 
        Compile[{{initial,  _Real,  1}, {p,  _Real,  0},  {q,  _Real, 
              0},  {r,  _Real,  0}}, 
Take[
            p*RotateRight[initial]  + r*initial + 
              q*RotateLeft[initial], {2, -2}]];
CAEuropeanOptionCompEvolveList[initial_List, p_, q_, r_, t_Integer] := 
        NestList[CompCAEuropeanOptionEvolveStep[#, p, q, r] &, initial, 
          t];

Of course both procedures can be collapsed into one line of code:

CompCAEuropeanOptionEvolveList  = 
        Compile[{{initial,  _Real,  1}, {p,  _Real,  0},  {q,  _Real, 
              0},  {r,  _Real,  0}, {t, _Integer, 
              0}}, 
NestList[
            p*RotateRight[#] + r*# + q*RotateLeft[#] &, initial, 
            t]];

We can also avail ourselves of the inbuilt function ListCorrelate (or ListConvolve) which carry out discrete convolution of the kernel of probabilities {p,q,r} with the state of the system. This has the same effect as one step of the previous procedure but can be more optimised for large arrays. Again NestList is applied for repeated application of the basic convolution as we move through time.

CompCAEuroEvolveStep  = 
        Compile[{{initial,  _Real,  1}, {p,  _Real,  0},  {q,  _Real, 
              0},  {r,  _Real,  0}}, 
ListCorrelate[{p, 
              r, q}, initial]];
CAEuroCompEvolveList[initial_List, p_, q_, r_, t_Integer] := 
        NestList[CompCAEuroEvolveStep[#, p, q, r] &, initial, 
          t];

As before both procedures can be collapsed into one line of code. But with the minor adjustment required by the fact that compiled functions only work on arrays of full rank and which therefore require the further application of the padLeftRight[.] function.

CompCAEuroEvolveList  = 
        Compile[{{initial,  _Real,  1}, {p,  _Real,  0},  {q,  _Real, 
              0},  {r,  _Real,  0}, {t, _Integer, 
              0}}, 
NestList[
            padLeftRight[ListCorrelate[{p, r, q}, #], Length[#]] &, initial, 
            t], {{padLeftRight[_], _Real, 1}}];

Finally we consider the use of the inbuilt function CellularAutomaton[.]. This undoubtedly constitutes the simplest use of the CA methodology and is also highly optimised. Again we are using the specification that requires a function in the first argument that convolves the probability vector {p,q,r} with the current state of the system.

EuroCAEvolveStep[initial_List, p_, q_, r_] := 
        CellularAutomaton[{# . {p, r, q} &, {}, 1}, initial, 1];

The use of the function CellularAutomaton[.] removes the need for further use of the NestList operator. However we do not need to carry all the results as we move to the final step of the automaton and we therefore specify the inclusion of only the last two time steps and cells ranging from positions t-1 to t+1.

EuroCAEvolveList[initial_List, p_, q_, r_, t_Integer] := 
        CellularAutomaton[{# . {p, r, q} &, {}, 1}, initial, 
          t, {-2, {t - 1, t + 1}}];

In order to create the european option programs we merely have to take the existing trinomial function pattern and replace the function ManyTimeStepEuroTri[.] wherever we see it with the new CA functions CompCAEuropeanOptionEvolveList, CompCAEuroEvolveList or EuroCAEvolveList. The following functions not only determine the european option values but with very little further calculation, all of the sensitivity functions or greeks also.

CompCAEuropeanOptionCall[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[S s^node u^(n-node)-K,0], {node, 0,n, 1}],
		 Table[Max[S d^node s^(n-node)-K,0], {node, 1,n, 1}]];
result = CompCAEuropeanOptionEvolveList[payoffin, P, Q, R, n][[{-2,-1},{n,n+1,n+2}]];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

CompCAEuropeanOptionPut[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[K-S s^node u^(n-node),0], {node, 0,n, 1}],
		 Table[Max[K-S d^node s^(n-node),0], {node, 1,n, 1}]];
result = CompCAEuropeanOptionEvolveList[payoffin, P, Q, R, n][[{-2,-1},{n,n+1,n+2}]];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

CompCAEuroEvolveCall[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[S s^node u^(n-node)-K,0], {node, 0,n, 1}],
		 Table[Max[S d^node s^(n-node)-K,0], {node, 1,n, 1}]];
result = CompCAEuroEvolveList[payoffin, P, Q, R, n][[{-2,-1},{n,n+1,n+2}]];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

CompCAEuroEvolvePut[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[K-S s^node u^(n-node),0], {node, 0,n, 1}],
		 Table[Max[K-S d^node s^(n-node),0], {node, 1,n, 1}]];
result = CompCAEuroEvolveList[payoffin, P, Q, R, n][[{-2,-1},{n,n+1,n+2}]];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

EuroCAEvolveCall[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[S s^node u^(n-node)-K,0], {node, 0,n, 1}],
		 Table[Max[S d^node s^(n-node)-K,0], {node, 1,n, 1}]];
result = EuroCAEvolveList[payoffin, P, Q, R, n];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

EuroCAEvolvePut[S_, K_, n_, r_, q_, σ_, T_] :=
Module[{a = afunc[n,r,q,T],b = bfunc[n, σ, T],
disc = Discount[n,r,T],u, d, s, P, Q, R, result,payoffin, value, delta, gamma, theta, rho, vega},
u = up[a,b];
d = down[a,b];
s = stay[a,b];
P = disc*pfunc[b];
Q = disc*qfunc[b];
R = disc*rfunc[b];
payoffin =Join[
		 Table[Max[K-S s^node u^(n-node),0], {node, 0,n, 1}],
		 Table[Max[K-S d^node s^(n-node),0], {node, 1,n, 1}]];
result = EuroCAEvolveList[payoffin, P, Q, R, n];
value = result[[-1,2]];
delta = (result[[1,1]] - result[[1,3]])/((u-d)*S);
gamma = ((result[[1,1]] - result[[1,2]])/(u - s) -
(result[[1,2]] - result[[1,3]])/(s - d))*2/(S^2 (u-d));
theta = (result[[1,2]]-(s-1)*S*delta- result[[-1,2]])*n/T;
rho = -T*(value - S*delta);
vega = T*S*S*gamma*σ;
{value, delta, gamma, theta, rho, vega}]

In the following calculations we determine the prices of european options and their sensitivity functions, known as the greeks, for the following market values: The stock price is S = $100, the strike price is K = $102, the risk free interst rate r = 10%, the dividend yield q = 3%, volatility σ = 30% and the time remaining tothe expiry date is T = 1 year. Here we choose n = 500 time steps.

\({BlackScholesCall[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], \ 
      BlackScholesCallDelta[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], 
      BlackScholesCallGamma[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], 
	
      BlackScholesCallTheta[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesCallRho[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesCallVega[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1]}\)

CompCAEuropeanOptionCall[100, 102,  500,  0.1,  0.03,  0.3,  
      1]

CompCAEuroEvolveCall[100, 102,  500,  0.1,  0.03,  0.3,  
      1]

EuroCAEvolveCall[100, 102,  500,  0.1,  0.03,  0.3,  1]

\({BlackScholesPut[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], \ 
      BlackScholesPutDelta[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], 
      BlackScholesPutGamma[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 1], 
	
      BlackScholesPutTheta[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesPutRho[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1], 
		BlackScholesPutVega[100, \ 102, \ 0.3, \ 0.1, \ 0.03, \ 
        1]}\)

CompCAEuropeanOptionPut[100, 102,  500,  0.1,  0.03,  0.3,  
      1]

CompCAEuroEvolvePut[100, 102,  500,  0.1,  0.03,  0.3,  
      1]

EuroCAEvolvePut[100, 102,  500,  0.1,  0.03,  0.3,  1]

American Option Theory

American options are intrinsically more difficult to analyse because of their property of early exercise, giving the holder of the option the right to exercise her option before the expiry date if she believes that it is unlikely that the underlying stock is either going to increase in the case of an american call or decrease in the case of an american put. The main emphasis of the solution of the Black-Scholes formula is that the stock has a lognormal distibution which menas that the expected future value at one time point - the maturity date can be reasonably described by a normal cumulative distribution function. But if the option can be exercised at any future time then there are an infinite number of such eventualities to calculate and the original argument breaks down. There are three alternative approaches however: a backward recursive binomial or trinomial model which can give approximate results for very large iterations, a quadratic approximation technique due to Baroni-Adesi and Whaley or lastly the determination of the moving free boundary, which is an associated problem concerning the structure of the time dependent curve which desribes those moments when exercise is optimal and the value of the option is best realised and placed in a bank account rather than left to depreciate with time. The recursive trinomial approach is described by the functions SSFastTriAmerCallFun[.] and SSFastTriAmerPutFun[.], but involves the use of procedural For[.] loops. It will be shown in the following section how these procedural constructs can be replaced by continuous cellular automata.

We will first give the trinomial code for the american call and put options and then show that we can also replicate these values from an appropriately defined continuous cellular automata.

OneTimeStepAmerFunTri = 
        Compile[{{initial,  _Real,  1}, {payoff, _Real, 1}, {p,  _Real,  
              0},  {q,  _Real, 0},  {r,  _Real,  
              0}}, 
MapThread[
            Max, {payoff, 
              p*Drop[initial, -2]  + 
                q*Drop[RotateLeft[initial, 2],  -2] + 
                r*Drop[RotateLeft[initial],  -2]}]];

\(\(ManyTimeStepAmerCallFunTri\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, 
              0}, \ {r, \ _Real, 0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 
              0}, \ {d, \ _Real, 0}, \ {s, \ _Real, 0}, {n, \ _Integer, \ 
              0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k = 1, \ R = 0, \ old = initial, \ 
              payoff = Drop[initial, \(-2\)], 
              finpay\  = \ {Max[S - K, 0]}}, 

            For[k\  = \ 1, 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = 
                  Join[Table[
                      Max[S\ s^R\ u^\((n - k - R)\) - K, 0], \ {R, \ 0, 
                        n - k, 1}], 
                    Table[Max[S\ d^R\ s^\((n - k - R)\) - K, 0], \ {R, 1, 
                        n - k, 1}]]; 
old\  = \ 
                  new; 
new\  = \ 
                  OneTimeStepAmerFunTri[old, payoff, \ p, q, r])\)]; 

            old\  = \ new; 

            new\  = \ \ OneTimeStepAmerFunTri[old, finpay, p, q, r]; 

            Join[new, old]]];\)\)

\(\(ManyTimeStepAmerPutFunTri\  = \ 

        Compile[{{initial, \ _Real, \ 1}, \ {p, _Real, \ 0}, \ {q, \ _Real, 
              0}, \ {r, \ _Real, 0}, \ {S, \ _Real, \ 0}, \ {u, \ _Real, \ 
              0}, \ {d, \ _Real, 0}, \ {s, \ _Real, 0}, {n, \ _Integer, \ 
              0}, \ {K, \ _Real, \ 0}}, 

          Module[{new\  = \ initial, \ k = 1, \ R = 0, \ old = initial, \ 
              payoff = Drop[initial, \(-2\)], 
              finpay\  = \ {Max[K - S, 0]}}, 

            For[k\  = \ 1, 
              k <= n - 1, \ \(k++\), \ 
\((payoff\  = 
                  Join[Table[
                      Max[K - S\ s^R\ u^\((n - k - R)\), 0], \ {R, \ 0, 
                        n - k, 1}], 
                    Table[Max[K - S\ d^R\ s^\((n - k - R)\), 0], \ {R, 1, 
                        n - k, 1}]]; 
old\  = \ 
                  new; 
new\  = \ 
                  OneTimeStepAmerFunTri[old, payoff, \ p, q, r])\)]; 

            old\  = \ new; 

            new\  = \ \ OneTimeStepAmerFunTri[old, finpay, p, q, r]; 

            Join[new, old]]];\)\)

\(SSFastTriAmerCallFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, payoffin}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
payoffin = 
          Join[Table[
              Max[S\ s^node\ u^\((n - node)\) - K, 0], {node, 0, n, 1}], 
            Table[Max[S\ d^node\ s^\((n - node)\) - K, 0], {node, 1, n, 
                1}]]; 
result = 
          ManyTimeStepAmerCallFunTri[payoffin, P, Q, R, S, u, d, s, n, 
            K]; 
value = 
          result[\([1]\)]; 
delta = \((result[\([2]\)] - 
                result[\([4]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([2]\)] - result[\([3]\)])\)/\((u - 
                      s)\) - \((result[\([3]\)] - result[\([4]\)])\)/\((s - 
                      d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([3]\)] - \((s - 1)\)*S*delta - result[\([1]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)

\(SSFastTriAmerPutFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, payoffin}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
payoffin = 
          Join[Table[
              Max[K - S\ s^node\ u^\((n - node)\), 0], {node, 0, n, 1}], 
            Table[Max[K - S\ d^node\ s^\((n - node)\), 0], {node, 1, n, 
                1}]]; 
result = 
          ManyTimeStepAmerPutFunTri[payoffin, P, Q, R, S, u, d, s, n, 
            K]; 
value = 
          result[\([1]\)]; 
delta = \((result[\([2]\)] - 
                result[\([4]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([2]\)] - result[\([3]\)])\)/\((u - 
                      s)\) - \((result[\([3]\)] - result[\([4]\)])\)/\((s - 
                      d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([3]\)] - \((s - 1)\)*S*delta - result[\([1]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)

The points to notice about the above code are that it is decidely procedural and that both option functions call either the compiled general procedure ManyTimeStepAmerCallFunTri[.] or ManyTimeStepAmerPutFunTri[.], both of which in turn repeatedly calls the single step function OneTimeStepAmerFunTri[.].

In the following calculations we determine the prices of american options and their sensitivity functions, known respectively as the delta, gamma, theta, rho and vega greeks, for the following market values: The stock price is S = $100, the strike price is K = $102, the risk free interst rate r = 10%, the dividend yield q = 0, volatility σ = 30% and the time remaining tothe expiry date is T = 1 year. Here we choose n = 500 time steps.

SSFastTriAmerCallFun[100, 102, 500, 0.1, 0, 0.3, 1]

SSFastTriAmerPutFun[100, 102, 500, 0.1, 0, 0.3, 1]

If we now load in a package designed to calculate the above described theoretical values we can then see how well the trinomial tree schemes compare. The quadratic approximation is given by the function AmericanPutMBAW[.],while the moving boundary approach due to Carr,Jarrow and Myneni is desribed by the functions AmericanPutCJM[.],AmericanPutDeltaCJM[.] and AmericanPutGammaCJM[.].These last three describe the put option value,as well as its delta and gamma greek sensitivities.

Needs["Derivatives`AnalyticAmerican`"];

?AmericanPutCJM

?AmericanPutMBAW

\({AmericanPutMBAW[100, 102, \ 0.3, \ 0.1, \ 0, \ 1], \ 
      AmericanPutCJM[100, 102, \ 0.3, \ 0.1, \ 1], 
      AmericanPutDeltaCJM[100, 102, \ 0.3, \ 0.1, \ 1], 
	
      AmericanPutGammaCJM[100, 102, \ 0.3, \ 0.1, \ 1]}\)

Clearly the trinomial model is giving accurate results up to the third decimal place for the Carr,Jarrow and Myneni model, but differs significantly from the Baroni-Adesi and Whaley quadratic formula, which is to be expected as it is only a second order approximation.

SSFastTriAmerPutFun[100, 102, 500, 0.1, 0, 0.3, 1]

Continuous CAs for American Options

Unlike the european options we have to not only specify the initial payoff function but also the payoff structure which reflects the possible exercise values at different nodes in the trinomial lattice. These are not the final values of the lattice but are updated when compared with the expected values determined by a recursion technique. In this way the american options are conditional cellular automata, with an initial set of values reflecting the environment dictated by the overall payoff structure of path dependent options such as american derivatives. These payoffs will be later inserted in the code for the final continuous cellular american options.

With[{u = 1.1, d = 0.9, s = 1.002, S = 100, K = 102, n = 20}, 
        callpayoffs  = 
          Table[Join[
              Table[Max[S s^node u^(n - k - node) - K, 0],  {node,  
                  0, n - k, 1}], 
              Table[Max[S d^node s^(n - k - node) - K, 0],  {node, 1, 
                  n - k, 1}]], {k, 0, n}]];

With[{u = 1.1, d = 0.9, s = 1.002, S = 100, K = 102, n = 20}, 
        putpayoffs  = 
          Table[Join[
              Table[Max[K - S s^node u^(n - k - node), 0],  {node,  
                  0, n - k, 1}], 
              Table[Max[K - S d^node s^(n - k - node), 0],  {node, 1, 
                  n - k, 1}]], {k, 0, n}]];

As with the european options we hall need to use the padLeftRight[.] function to describe individual payoff values and paddedPayoffs[.] to denote the construction of payoff structures appropriate to american options.   The main purpose of these supplementary functions padLeftRight and paddedPayoffs  is to fill out the tree structure into a rectangular array suitable for repeated application of the EvolveStep procedure.

padLeftRight[list_List, n_Integer] := 
        PadLeft[list, 
          n, {0}, (n - Length[list])/2];
paddedPayoffs[payoff_List, 
        payoffs_List] := padLeftRight[#, Length[payoff]] &  /@  
        payoffs

As before there are a number of ways we can construct the appropriate cellular automata. Either we can specify the structure of the Step and Evolve List functions ourselves, or we can utilise the inbuilt CellularAutomaton[.] function. In the first case we will replicate the behaviour of the OneTimeStepAmerFunTri function and build a compiled procedure that shifts the previous array of cells depending upon whether we have moved into the up, down states or remained in the same state. Suppose that we are in the stae called initial, then at the next step we will be at the same price level with probability r we will be at the level down with probability q or finally we will be at the upper level with probability p. If we are at the same price level then the prices remain unchanged, whereas if we go down one level then the second value now occupies the leading position and is therefore shifted to the left, but if we go one level up then the first element is forced into the second position. In this way the expected future value of the next step of the process is given by:  p*RotateRight[initial] +r*initial+q*RotateLeft[initial]. In order that first and last values at either end do not get incorrectly transposed, we remove them successively from each step by applying the procedure Take[.., {2,-2}]. However american options have the additional attribute of early exercise which translates into the payoff at any particular time being the maximum of either the expected future value or the actual exercised value. This necessitates mapping the Max function over the exercise payoff values and comparing them with the expected values. This is accomplished with the MapThread function and then repeatedly applied with the FoldList command.

CompCAAmericanOptionEvolveStep  = 
        Compile[{{initial,  _Real,  1}, {payoff, _Real, 1}, {p,  _Real,  
              0},  {q,  _Real, 0},  {r,  _Real,  
              0}}, 
MapThread[
            Max, {payoff, 
              Take[p*RotateRight[initial]  + r*initial + 
                  q*RotateLeft[initial], {2, -2}]}]];
CAAmericanOptionCompEvolveList[initial_List, payoffs_List, p_, q_, 
          r_] := 
FoldList[
          CompCAAmericanOptionEvolveStep[#1, #2, p, q, r] &, initial, 
          payoffs];

Of course both procedures can be collapsed into one line of code:

CompCAAmericanOptionEvolveList = 
        Compile[{{initial,  _Real,  1}, {paddedPayoffs, _Real, 
              2}, {p,  _Real,  0},  {q,  _Real, 0},  {r,  _Real,  0}}, 
          FoldList[
            MapThread[
                Max, {p*RotateRight[#1] + r*#1 + q*RotateLeft[#1], #2}] &, 
            initial, paddedPayoffs]];

We can also avail ourselves of the inbuilt function ListCorrelate (or ListConvolve) which carries out discrete convolution of the kernel of probabilities {p,q,r} with the state of the system. This has the same effect as one step of the previous procedure but can be more optimised for large arrays. Again FoldList is applied for repeated application of the basic convolution as we move through time.

CompCAAmerEvolveStep  = 
        Compile[{{initial,  _Real,  1}, {payoff, _Real, 1}, {p,  _Real,  
              0},  {q,  _Real, 0},  {r,  _Real,  
              0}}, 
MapThread[
            Max, {payoff, ListCorrelate[{p, r, q}, initial]}]];
CAAmerCompEvolveList[initial_List, payoffs_List, p_, q_, 
          r_] := 
FoldList[
          CompCAAmerEvolveStep[#1, #2, p, q, r] &, initial, 
          payoffs];

As before both procedures can be collapsed into one line of code. But with the minor adjustment required by the fact that compiled functions only work on arrays of full rank and which therefore require the further application of the padLeftRight[.] function. Note also that the time dependent payoff structure is also padded out to match an array of full rank.

CompCAAmerEvolveList = 
        Compile[{{initial,  _Real,  1}, {paddedPayoffs, _Real, 
              2}, {p,  _Real,  0},  {q,  _Real, 0},  {r,  _Real,  0}}, 
          FoldList[
            MapThread[
                Max, {padLeftRight[ListCorrelate[{p, r, q}, #1], 
                    Length[#1]], #2}] &, initial, 
            paddedPayoffs], {{padLeftRight[_], _Real, 1}}];

Finally we consider the use of the inbuilt function CellularAutomaton[.]. This undoubtedly constitutes the simplest use of the CA methodology and is also highly optimised. Again we are using the specification that requires a function in the first argument that convolves the probability vector {p,q,r} with the current state of the system. But unlike the situation encountered earlier for european options, we cannot avoid using the FoldList command as the lattice values have to be updated relative to the time dependent payoff structure.

AmerCAEvolveStep[initial_List, payoff_List, p_, q_, r_] := 
      MapThread[
        Max, {padLeftRight[payoff, 
            Length[initial]], CellularAutomaton[{# . {p, r, q} &, {}, 1}, 
              initial, 1][[-1]]}]
AmerCAEvolveList[initial_List, payoffs_List, p_, q_, r_] := 
      FoldList[AmerCAEvolveStep[#1, #2, p, q, r] &, initial, 
        payoffs]

Again both procedures can be collapsed into one line of code using pure functions.

AmerCAEvolveList[initial_List, payoffs_List, p_, q_, r_] := 
      FoldList[Function[{alist, blist}, 
          MapThread[
            Max, {padLeftRight[blist, 
                Length[alist]], CellularAutomaton[{# . {p, r, q} &, {}, 1}, 
                  alist, 1][[-1]]}] ], initial, 
        payoffs]

In order to create the american option programs we merely have to take the existing trinomial function pattern and replace the functions ManyTimeStepAmerCallFunTri[.] and ManyTimeStepAmerPutFunTri[.] wherever we see them with the new CA functions CompCAAmericanOptionEvolveList, CompCAAmerEvolveList or AmerCAEvolveList in conjunction with the callpayoff and putpayoff matrix structure. The following functions not only determine the american option values but with very little further calculation, all of the sensitivity functions or greeks also.

\(\(\(\ \)\(CompCAAmericanOptionEvolveCallFun[S_, K_, n_, r_, 
        q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, callpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
callpayoffs\  = 
          Table[Join[
              Table[Max[S\ s^node\ u^\((n - k - node)\) - K, 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[S\ d^node\ s^\((n - k - node)\) - K, 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(CompCAAmericanOptionEvolveList[callpayoffs[\([1]\)], 
              paddedPayoffs[callpayoffs[\([1]\)], Rest[callpayoffs]], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)\)\)

\(CompCAAmericanOptionEvolvePutFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, putpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
putpayoffs\  = 
          Table[Join[
              Table[Max[K - S\ s^node\ u^\((n - k - node)\), 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[K - S\ d^node\ s^\((n - k - node)\), 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(CompCAAmericanOptionEvolveList[putpayoffs[\([1]\)], 
              paddedPayoffs[putpayoffs[\([1]\)], Rest[putpayoffs]], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)

\(\(\(\ \)\(CompCAAmerEvolveCallFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, callpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
callpayoffs\  = 
          Table[Join[
              Table[Max[S\ s^node\ u^\((n - k - node)\) - K, 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[S\ d^node\ s^\((n - k - node)\) - K, 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(CompCAAmerEvolveList[
              callpayoffs[\([1]\)], 
              paddedPayoffs[callpayoffs[\([1]\)], Rest[callpayoffs]], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)\)\)

\(CompCAAmerEvolvePutFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, putpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
putpayoffs\  = 
          Table[Join[
              Table[Max[K - S\ s^node\ u^\((n - k - node)\), 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[K - S\ d^node\ s^\((n - k - node)\), 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(CompCAAmerEvolveList[
              putpayoffs[\([1]\)], 
              paddedPayoffs[putpayoffs[\([1]\)], Rest[putpayoffs]], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)

\(\(\(\ \)\(AmerCAEvolveCallFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, callpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
callpayoffs\  = 
          Table[Join[
              Table[Max[S\ s^node\ u^\((n - k - node)\) - K, 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[S\ d^node\ s^\((n - k - node)\) - K, 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(AmerCAEvolveList[
              callpayoffs[\([1]\)], Rest[callpayoffs], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)\)\)

\(AmerCAEvolvePutFun[S_, K_, n_, r_, q_, σ_, T_] := 
      Module[{a = afunc[n, r, q, T], b = bfunc[n, σ, T], 
          disc = Discount[n, r, T], P, Q, R, value, delta, gamma, theta, rho, 
          vega, result, putpayoffs, node}, 
u = 
          up[a, b]; 
\ \ \ \ d = down[a, b]; 
\ \ \ \ s = 
          stay[a, b]; 
P = 
          disc\ pfunc[b]; 
Q = 
          disc\ qfunc[b]; 
R = 
          disc\ rfunc[b]; 
putpayoffs\  = 
          Table[Join[
              Table[Max[K - S\ s^node\ u^\((n - k - node)\), 0], \ {node, \ 
                  0, n - k, 1}], 
              Table[Max[K - S\ d^node\ s^\((n - k - node)\), 0], \ {node, 1, 
                  n - k, 1}]], {k, 0, 
              n}]; 
result\  = \ \(AmerCAEvolveList[
              putpayoffs[\([1]\)], Rest[putpayoffs], P, Q, 
              R]\)[\([{\(-2\), \(-1\)}, {n, n + 1, 
                n + 2}]\)]; 
value = 
          result[\([\(-1\), 
              2]\)]; 
delta = \((result[\([1, 1]\)] - 
                result[\([1, 3]\)])\)/\((\((u - d)\)*S)\); 
        gamma = \((\((result[\([1, 1]\)] - result[\([1, 2]\)])\)/\((u - 
                      s)\) - \((result[\([1, 2]\)] - 
                      result[\([1, 3]\)])\)/\((s - d)\))\)*2/\((S^2\ \((u - 
                      d)\))\); 
theta = \((result[\([1, 
                    2]\)] - \((s - 1)\)*S*delta - result[\([\(-1\), 2]\)])\)*
            n/T; 
rho = \(-T\)*\((value - 
                S*delta)\); 
vega = 
          T*S*S*gamma*σ; 
{value, delta, gamma, 
          theta, rho, vega}]\)

On[General::"spell", General::"spell1"]

In the following calculations we determine the prices of american options and their sensitivity functions, known as the greeks, for the following market values: The stock price is S = $100, the strike price is K = $102, the risk free interst rate r = 10%, the dividend yield q = 3%, volatility σ = 30% and the time remaining tothe expiry date is T = 1 year. Here we choose n = 30 time steps.

SSFastTriAmerCallFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

CompCAAmericanOptionEvolveCallFun[100, 102,  30,  0.1,  0.03,  
      0.3,  1]

CompCAAmerEvolveCallFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

AmerCAEvolveCallFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

SSFastTriAmerPutFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

CompCAAmericanOptionEvolvePutFun[100, 102,  30,  0.1,  0.03,  
      0.3,  1]

CompCAAmerEvolvePutFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

AmerCAEvolvePutFun[100, 102,  30,  0.1,  0.03,  0.3,  
      1]

Graphics

European

ListPlot3D[
      Table[EuroTrinomialCall[S, 102,  30,  0.1,  0.03,  σ,  
            1][[1]], {S, 80. , 120. , 5. }, {σ, 0.2, 
          0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Euro Call"}, 
      PlotLabel -> "Euro Trinomial Call"]

ListPlot3D[
      Table[CompCAEuropeanOptionCall[S, 102,  30,  0.1,  
            0.03,  σ,  1][[1]], {S, 80. , 120. , 
          5. }, {σ, 0.2, 0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Euro Call"}, 
      PlotLabel -> "CA Euro Call"]

ListPlot3D[
      Table[CompCAEuroEvolveCall[S, 102,  30,  0.1,  0.03,  σ,  
            1][[1]], {S, 80. , 120. , 5. }, {σ, 0.2, 
          0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Euro Call"}, 
      PlotLabel -> "CA Evolve Euro Call"]

ListPlot3D[
      Table[EuroCAEvolveCall[S, 102,  30,  0.1,  0.03,  σ,  
            1][[1]], {S, 80. , 120. , 5. }, {σ, 0.2, 
          0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Euro Call"}, 
      PlotLabel -> "CA Evolve Euro Call"]

American

ListPlot3D[
      Table[SSFastTriAmerCallFun[S, 102,  30,  0.1,  0.03,  σ,  
            1][[1]], {S, 80. , 120. , 5. }, {σ, 0.2, 
          0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Amer Call"}, 
      PlotLabel -> "Amer Trinomial Call"]

ListPlot3D[
      Table[CompCAAmericanOptionEvolveCallFun[S, 102,  30,  0.1,  
            0.03,  σ,  1][[1]], {S, 80. , 120. , 
          5. }, {σ, 0.2, 0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Amer Call"}, 
      PlotLabel -> "Amer CA Call"]

ListPlot3D[
      Table[CompCAAmerEvolveCallFun[S, 102,  30,  0.1,  
            0.03,  σ,  1][[1]], {S, 80. , 120. , 
          5. }, {σ, 0.2, 0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Amer Call"}, 
      PlotLabel -> "Amer CA Evolve Call"]

ListPlot3D[
      Table[AmerCAEvolveCallFun[S, 102,  30,  0.1,  0.03,  σ,  
            1][[1]], {S, 80. , 120. , 5. }, {σ, 0.2, 
          0.4,  .05}], 
      AxesLabel -> {"Stock", "volatility", "Amer Call"}, 
      PlotLabel -> "Amer CA Evolve Call"]