This notebook implements binomial option pricing models for european and american style options.
See also https://github.com/JuliaQuant/FinancialDerivatives.jl for a more ambitious derivatives package (not used here).
In [1]:
using Dates
include("jlFiles/printmat.jl")
include("jlFiles/printTable.jl")
Out[1]:
In [2]:
using Plots
#pyplot(size=(600,400))
gr(size=(480,320))
default(fmt = :svg)
In a binomial model the underlying can change from $S$ today to either $Su$ or $Sd$ in the next period (which is $h$ years from now).
Let $f_u$ and $f_d$ be the values of the derivative in the up- and down-state next period. Then, the value today is
$f =e^{-yh}\left[ pf_{u}+\left( 1-p\right) f_{d}\right] \: \text{ with } \: p=\frac{e^{yh}-d}{u-d}$
For a call option that expires next period, the payoffs in the two states are
$ f_u = \max(Su-K,0) $
$ f_d = \max(Sd-K,0) $
In [3]:
S = 10 #underlying price today
u = 1.1
d = 0.95
K = 10 #strike price: try different values
y = 0 #interest rate
h = 1/12 #time to expiration
fu = max(S*u-K,0) #value of call option in up node
fd = max(S*d-K,0) #in down node
p = (exp(y*h) - d)/(u-d)
C = exp(-y*h)*(p*fu+(1-p)*fd)
printblue("Pricing of a call option with strike $K, one time step:")
xx = [fu,fd,p,C]
printTable(xx,[""],["Payoff up","Payoff down","p","call price"])
The CRR approach to construct $(u,d)$ for a small time step of length $h$ is
$u=e^{\sigma\sqrt{h}} \: \text{ and }\: d=e^{-\sigma\sqrt{h}}$,
where $\sigma$ is the annualized standard deviation of the underlying. Notice that $p$ depends on the choice of $(u,d)$
$ p=\frac{e^{yh}-d}{u-d}. $
In [4]:
m = 0.5 #time to expiration (in years)
y = 0.05 #interest rate (annualized)
σ = 0.2 #annualized std of underlying asset
n = 50 #number of time steps of size h
h = m/n #time step size (in years)
u = exp(σ*sqrt(h))
d = exp(-σ*sqrt(h))
p = (exp(y*h) - d)/(u-d) #recalibrate p, it depends on u and d
printblue("CRR parameters in tree when m=$m, y=$y, σ=$σ and n=$n:")
xx = [h,u,d,p,exp(y*h)]
printTable(xx,[""],["h","u","d","p","exp(yh)"])
printblue("Checking if u > exp(y*h) > d: $(u > exp(y*h) > d)")
In [5]:
x = [ones(i)*i for i = 1:3] #vector or vectors (of different lengths)
printblue("Illustrating a vector of vectors:\n")
for i = 1:length(x)
printblue("i=$i:")
printmat(x[i])
end
We create a tree by starting at the current spot price $S$. The future nodes are then created by multiplying by $u$ or $d$
$\text{step 0}: (S)$
$\text{step 1}: (Su,Sd)$
$\text{step 2}: (Suu,Sud,Sdd)$ (since $Sud=Sdu$)
...
Clearly, the vectors are of different lengths at the different time steps.
For step 2, the code in the next cell creates the vector $(Suu,Sud,Sdd)$ by first multiplying the step 1 vector by the up factor ($u(Su,Sd)$) and then attaching (as the last element) the last element of the step 1 vector times the down factor ($dSd$).
Notice that the first index of a traditional vector is 1, so STree[1] corresponds to time step 0, and STree[n+1] to time step $n$.
In [6]:
"""
BuildSTree(S,n,u,d)
Build binomial tree, starting at `S` and having `n` steps with up move `u` and down move `d`
# Output
- `STree:: Vector or vectors`: each (sub-)vector is for a time step
"""
function BuildSTree(S,n,u,d)
STree = [fill(NaN,i) for i = 1:n+1] #vector of vectors (of different lengths)
STree[1][1] = S #step 0 is in STree[1], element 1
for i = 2:n+1 #move forward in time
STree[i][1:end-1] = u*STree[i-1] #up move from STree[i-1][1:end]
STree[i][end] = d*STree[i-1][end] #down move from STree[i-1][end]
end
return STree
end
Out[6]:
In [7]:
S = 42.0 #current price of underlying asset
STree = BuildSTree(S,n,u,d)
printblue("printing the first 3 time steps in the tree")
for i = 1:3
printblue("i=$i:")
printmat(STree[i])
end
In [8]:
p1 = plot( legend = false,
xlim = (-0.1,2.5),
ylim = (40,44),
xticks = 0:2,
title = "The first few nodes in the tree",
xlabel = "time step",
ylabel = "value of asset",
annotatation = (0,43.8,text("Structure of STree[timestep][node]",8,:left)) )
for i = 1:3, j = 1:length(STree[i])
scatter!( [i-1],[STree[i][j]], #[] to make a 1-element vector
annotation = (i-1+0.1,STree[i][j],text("STree[$i][$j]",8,:left)) )
end
display(p1) #needed since plot() does not have data points to plot
In [9]:
p1 = plot( legend = false,
xlim = (-0.1,2.5),
ylim = (40,44),
xticks = 0:2,
title = "The first few nodes in STree",
xlabel = "time step",
ylabel = "value of asset" )
for i = 1:3, j = 1:length(STree[i])
annotate!(i-1,STree[i][j],text(string(round(STree[i][j],digits=2)),8))
end
display(p1)
In [10]:
p1 = scatter( [0.0],STree[1],
markercolor = :blue,
xlim = (-2,n+1.5),
ylim = (0,115),
legend = false,
title = "All nodes in STree",
xlabel = "time step",
ylabel = "value of asset" )
for i = 1:n+1
scatter!(fill(i-1,i),STree[i],color=:blue)
end
display(p1)
Let $f_{ij}$ be the option price at time step $i$ when the underlying price is $S_{ij}$. (We use $S_{ij}$ as a shorthand notation, to avoid writing things like $Sudd$.) In the code below, $S_{i1}$ is the highest node in time step $i$, $S_{i2}$ is the second highest, etc.
For a European call option, the call price at the last time step $n$ is $f_{nj} = \max(0,S_{nj}-K)$.
Instead, the put price is $f_{nj} = \max(0,K-S_{nj})$.
For all earlier time steps, the value is
$f_{ij} = e^{-yh}[p f_{i+1,j} + (1-p) f_{i+1,j+1}]$
For instance,
$f_{i1} = e^{-yh}[p f_{i+1,1} + (1-p) f_{i+1,2}]$
$f_{i2} = e^{-yh}[p f_{i+1,2} + (1-p) f_{i+1,3}]$
In [11]:
"""
EuOptionPrice(STree,K,y,h,p,isPut=false)
Calculate price of European option from binomial model
# Output
- `Value:: Vector of vectors`: option values at different nodes, same structure as STree
"""
function EuOptionPrice(STree,K,y,h,p,isPut=false) #price of European option
Value = deepcopy(STree) #tree for derivative, to fill
n = length(STree) - 1 #number of steps in STree
if isPut
Value[n+1] = max.(0,K.-STree[n+1]) #put, at last time node
else
Value[n+1] = max.(0,STree[n+1].-K) #call, at last time node
end
for i = n:-1:1 #move backward in time
Value[i] = exp(-y*h)*(p*Value[i+1][1:end-1] + (1-p)*Value[i+1][2:end])
end #p*up + (1-p)*down, discount
return Value
end
Out[11]:
In [12]:
K = 42.0 #strike price
Pe = EuOptionPrice(STree,K,y,h,p,true) #Pe[1] is a 1x1 vector with the put price (node 0)
Ce = EuOptionPrice(STree,K,y,h,p,false) #use Pe[1][1] or Pe[1][]
Ce_parity = Pe[1][] + S - exp(-m*y)*K #put-call parity, Pe[1][] is a scalar
printlnPs("European put and call prices at K=$K and S=$S: ",[Pe[1] Ce[1]])
printlnPs("\nCall price according to put-call parity: ",Ce_parity)
The option values are calculated as for the European option, except that that the option value is
$f_{ij} = \max(\text{value if exercised now},\text{continuation value})$
The continuation value has the same form as in the European case, and thus assumes that the option has not been exercised before the next period.
The value of exercising now is $S_{ij}-K$ for a call and $K-S_{ij}$ for a put.
In [13]:
"""
AmOptionPrice(STree,K,y,h,p,isPut=false)
Calculate price of American option from binomial model
# Output
- `Value:: Vector of vectors`: option values at different nodes, same structure as STree
- `Exerc::` Vector of vectors`: true if early exercise at the node, same structure as STree
"""
function AmOptionPrice(STree,K,y,h,p,isPut=false) #price of American option
Value = deepcopy(STree) #tree for derivative, to fill
n = length(STree) - 1
Exerc = [falses(i) for i = 1:length(STree)]
if isPut
Value[n+1] = max.(0,K.-STree[n+1]) #put, at last time node
else
Value[n+1] = max.(0,STree[n+1].-K) #call, at last time node
end
Exerc[n+1] = Value[n+1] .> 0 #exercise
for i = n:-1:1 #move backward in time
fa = exp(-y*h)*(p*Value[i+1][1:end-1] + (1-p)*Value[i+1][2:end])
if isPut
Value[i] = max.(K.-STree[i],fa) #put
else
Value[i] = max.(STree[i].-K,fa) #call
end
Exerc[i] = Value[i] .> fa #early exercise
end
return Value, Exerc
end
Out[13]:
In [14]:
K = 42.0 #strike price
(Pa,Exerc) = AmOptionPrice(STree,K,y,h,p,true)
Ca, = AmOptionPrice(STree,K,y,h,p,false)
printlnPs("American put and call prices at K=$K and S=$S: ",[Pa[1] Ca[1]])
printlnPs("\nCompare with European put and call prices: ",[Pe[1] Ce[1]])
In [15]:
p1 = plot( xlim = (-0.5,11),
ylim = (33,53),
xticks = 0:10,
legend = false,
title = "American put values at different nodes",
xlabel = "time step",
ylabel = "value of asset",
annotation = (0,50,text("Some nodes may not be\nreached due to early exercise",8,:left)) )
for i = 1:11, j = 1:length(STree[i])
annotate!(i-1,STree[i][j],text(string(round(Pa[i][j],digits=2)),8))
end
display(p1)
In [16]:
p1 = plot( xlim = (-2,n+1),
ylim = (0,115),
legend = false,
title = "Exercise nodes",
xlabel = "time step",
ylabel = "value of asset",
annotation = (0,100,text("Some nodes may not be\nreached due to early exercise",8,:left)) )
for i = 1:n+1, j = 1:length(STree[i])
if Exerc[i][j]
scatter!([i-1],[STree[i][j]],markercolor=:blue)
end
end
display(p1)
In [ ]: