where the probability is the so-called risk neutral probability

So, we’ve done it here 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: went up twice, went 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 parameter (the so-call risk neutral probability, if it is unique). But we do not really care about those closed formula, the goal is to write an algorithm which computes 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 parameters 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 (or subdivision of the time interval, from now till maturity of the European option),

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 it to price American options. 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 down (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. An if we move one step foreward, we have the following tree.

Here it is still recombining. But the size will increase much faster than in the univariate case. Now, assume that there might be some correlation. Then one can consider the following values, to have a specific correlation,

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 (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

Well, playing with trees is nice, but it might not be optimal for complex products. Next time, we’ll discuss other techniques…