I am a big fan of trees. It is a very nice way to see how financial pricing works, for derivatives. An with a matrix-based language (R for instance), it is extremely simple to compute almost everything. Even multiple assets options. Let us see how it works. But first, I have to assume that everyone knows about trees (or at least is familiar), and about risk neutral probabilities, and with standard financial derivatives. Just in case, I can upload some old slides of the first course on asset pricing we gave a few years ago at École Polytechnique.

Let us get back on the code to price a (European) call option with trees.The idea is simple. We have to fix the number of periods. Let us start with only one (as described in the slides above). The stock has price and can go either up, and then have price or go down, and have price . And the fundamental theorem of asset pricing says that we do not really care about probabilities of going up, or down. Assuming that we can buy or sell that stock, and that a risk free asset is available on the market, it is possible to price any contingent financial product, like a financial option. Since we know the final value of the option when the stock goes either up, or down, it is possible to replicate the payoff of that option using the stock and the risk free asset. And we can prove that the price of the option is simply

where the probability is the so-called risk neutral probability

So, we’ve done it with only one single period, but it is possible to extend it to multiperiods. The idea is to keep that multiplicative representation of possible values of the stock, and to get a recombinant tree. At step 2, the stock can take only three different values: when up twice, when down twice, or went up and down (or the reverse, but we don’t care, this is the point of recombining). If we write things down, then we can prove that

for some probability. But we do not really care about those closed formula, the goal is to write an algorithm which compute the tree, and return the price of a call option (say). But before starting, we have to make a connection between that model with up and down prices, and the parameter of the Black-Scholes diffusion, for the stock price. The idea is to identify the first and the second moment, i.e.

(where, under the risk neutral probability, the trend is the risk free rate) and

The code might look like that

n=5; T=1; r=0.05; sigma=.4;S=50;K=50
price=function(n){
u.n=exp(sigma*sqrt(T/n));
d.n=1/u.n
p.n=(exp(r*T/n)-d.n)/(u.n-d.n)
SJ=matrix(0,n+1,n+1)
SJ[1,1]=S
for(i in(2:(n+1)))
{for(j in(1:i)){SJ[i,j]=S*u.n^(i-j)*d.n^(j-1)}}
OPT=matrix(0,n+1,n+1)
OPT[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
for(i in(n:1))
{for(j in(1:i)){OPT[i,j]=exp(-r*T/n)*(OPT[i+1,j]*p.n+
(1-p.n)*OPT[i+1,j+1])}}
return(OPT[1,1])
}

We can plot the evolution of the price, as a function of the number of time periods.

N=10:400
V=Vectorize(price)(N)
plot(N,V,type="l")

Note that we can compare with the Black-Scholes price of this call option, given by

where

and

d1=1/(sigma*sqrt(T))*(log(S/K)+(r+sigma^2/2)*T)
d2=d1-sigma*sqrt(T)
BS=S*pnorm(d1)-K*exp(-r*T)*pnorm(d2)
abline(h=BS,lty=2,col="red")

The code is clearly not optimal, but at least, we see what’s going on. For instance, we do not need a matrix when we calculate using backward recursions the price of the option. We can just keep a single vector. But this matrix is nice, because we can use in to price American options.

But the great thing with trees, is that we can extend them to model comovments of two underlying prices. For instance, with the code below, we compare the price of an American put option, and the price of European put option.

price.american=function(n,opt="put"){
u.n=exp(sigma*sqrt(T/n)); d.n=1/u.n
p.n=(exp(r*T/n)-d.n)/(u.n-d.n)
SJ=matrix(0,n+1,n+1)
SJ[1,1]=S
for(i in(2:(n+1)))
{for(j in(1:i)) {SJ[i,j]=S*u.n^(i-j)*d.n^(j-1)}}
OPTe=matrix(0,n+1,n+1)
OPTa=matrix(0,n+1,n+1)
if(opt=="call"){
OPTa[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
OPTe[n+1,]=(SJ[n+1,]-K)*(SJ[n+1,]>K)
}
if(opt=="put"){
OPTa[n+1,]=(K-SJ[n+1,])*(SJ[n+1,]<K)
OPTe[n+1,]=(K-SJ[n+1,])*(SJ[n+1,]<K)
}
for(i in(n:1))
{
for(j in(1:i))
{if(opt=="call"){
OPTa[i,j]=max((SJ[i,j]-K)*(SJ[i,j]>K),
exp(-r*T/n)*(OPTa[i+1,j]*p.n+
(1-p.n)*OPTa[i+1,j+1]))}
if(opt=="put"){
OPTa[i,j]=max((K-SJ[i,j])*(K>SJ[i,j]),
exp(-r*T/n)*(OPTa[i+1,j]*p.n+
(1-p.n)*OPTa[i+1,j+1]))}
 
OPTe[i,j]=exp(-r*T/n)*(OPTe[i+1,j]*p.n+
(1-p.n)*OPTe[i+1,j+1])}}
priceop=c(OPTe[1,1],OPTa[1,1])
names(priceop)=c("E","A")
return(priceop)}

It is possible to compare those price, obtained on trees, with prices given by closed (approximated) formulas.

> d1=1/(sigma*sqrt(T))*(log(S/K)+(r+sigma^2/2)*T)
> d2=d1-sigma*sqrt(T)
> (BS=-S*pnorm(-d1)+K*exp(-r*T)*pnorm(-d2)  )
[1] 6.572947
> N=10:200
> M=Vectorize(price.american)(N)
> plot(N,M[1,],type='l',col='blue',ylim=range(M))
> lines(N,M[2,],type='l',col='red')
> abline(h=BS,lty=2,col='blue')
> library(fOptions)
> (am=BAWAmericanApproxOption(TypeFlag =
+ "p", S = S,X = K, Time = T, r = r,
+ b = r, sigma =sigma)@price)
[1] 6.840335
> abline(h=am,lty=2,col='red')

Another great thing with trees, is that it becomes possible to plot to region where it is optimal to exercise our right to sell the stock.

Let us move now to a model with two assets, as suggested by Rubinstein (1994). First, observe that a discretization of two independent Brownian motions will be based on two independent random walk, taking values

i.e. both went up (NW), both went dowm (SE), and one went up while the other went down (either NE or SW).

With independent and symmetric random walks, the probabilities will be respectively 1/4. Here, the idea is to consider possible correlation, i.e.

And again, the idea is then to identify the first two moments. This gives us the following system of equations for the four respective (risk neutral) probabilities 

For those willing to do the maths, please do. The answer should be

and for the last one

The code here looks like that

price.spead=function(n){
T=1; r=0.05; K=0
S1=105
S2=100
sigma1=0.4
sigma2=0.3
rho=0.5
u1.n=exp(sigma1*sqrt(T/n)); d1.n=1/u1.n
u2.n=exp(sigma2*sqrt(T/n)); d2.n=1/u2.n
 
v1=r-sigma1^2/2; v2=r-sigma2^2/2
puu.n=(1+rho+sqrt(T/n)*(v1/sigma1+v2/sigma2))/4
pud.n=(1-rho+sqrt(T/n)*(v1/sigma1-v2/sigma2))/4
pdu.n=(1-rho+sqrt(T/n)*(-v1/sigma1+v2/sigma2))/4
pdd.n=(1+rho+sqrt(T/n)*(-v1/sigma1-v2/sigma2))/4
k=0:n
un=matrix(1,n+1,1)
SJ= (S1 * d1.n^k * u1.n^(n-k-1)) %*% t(un) -
un %*%t(S2 * d2.n^k * u2.n^(n-k-1))
OPT=(SJ)*(SJ>K)
for(k in(n:1))
{
OPT0=matrix(0,k,k)
for(i in(1:k))
{
for(j in(1:k))
{OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)}}
OPT=OPT0}
return(OPT[1,1])}

If we look at the details, consider two periods, like on the figure above. The are nine values for the spread,

> n=2
> SJ
[,1]      [,2]       [,3]
[1,]  32.02217  84.86869 119.443578
[2,] -47.84652   5.00000  39.574891
[3,] -93.20959 -40.36308  -5.788184

and the payoff of the option is here

> OPT
[,1]     [,2]      [,3]
[1,] 32.02217 84.86869 119.44358
[2,]  0.00000  5.00000  39.57489
[3,]  0.00000  0.00000   0.00000

So if we go backward of one step, we have the following square of values

> k=n
> OPT0<-matrix(0,k,k)
> for(i in(1:k))
+ {
+   for(j in(1:k))
+   {
+     OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
+ OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)
+ }
+ }
> OPT0
[,1]      [,2]
[1,] 22.2741190 58.421275
[2,]  0.5305465  5.977683

The idea is then to move backward once more,

> OPT=OPT0
> OPT0<-matrix(0,k,k)
> for(i in(1:k))
+ {
+   for(j in(1:k))
+   {
+     OPT0[i,j]=(OPT[i,j]*puu.n+OPT[i+1,j]*pdu.n+
+ OPT[i,j+1]*pud.n+OPT[i+1,j+1]*pdd.n)*exp(-r*T/n)
+ }
+ }
> OPT0
[,1]
[1,] 16.44106

Here calculations are much longer,

> price.spead(250)
[1]  15.66496

and again, it is possible to use standard approximations to compare that price with a more standard one,

> (sp=SpreadApproxOption(TypeFlag =
+ "c", S1 = 105, S2 = 100, X = 0,
+ Time = 1, r = .05, sigma1 = .4,
+ sigma2 = .3, rho = .5)@price)
[1]  15.65077