Suppose there is an asset with the initial spot price \(S_0\) and it may go up by \(u\) to \(S_{1}=uS_0\) with probability \(p\) or by \(d\) to \(S_{1}=dS_0\) with probability \(1-p\) during the time interval of \(dt\). Consider the European Call option and the payoffs as \(C_u = [0,uS_0-K]^+\) and \(C_d = [0,dS_0-K]^+\).
We want to form a (replicating, risk-free) portfolio which has the same value regardless of the outcome (i.e. \(u\) or \(d\)). So, we add \(\Delta\) amount of the underlying asset to get the following equality.
\[-C_u + \Delta u S_0 = -C_d + \Delta d S_0\]
\(\Delta\) value can be calculated as \(\Delta = \dfrac{C_u - C_d}{(u-d)S_0}\).
This risk-free portfolio should return risk-free interest rate.
\[-C_u + \Delta u S_0 = e^{r.dt}(-C + \Delta S_0)\]
Substitute \(\Delta\) and rearrange to find \(C\) in the following way.
\[ C = \dfrac{C_u - C_d}{(u-d)S_0}S_0 + e^{-r.dt}(C_u - \dfrac{C_u - C_d}{(u-d)S_0}uS_0) \]
\[ C = \dfrac{C_u - C_d}{(u-d)} + e^{-r.dt}(C_u - \dfrac{C_u - C_d}{(u-d)}u) \]
\[ C = e^{-r.dt}(e^{r.dt}\dfrac{C_u}{(u-d)} - e^{r.dt}\dfrac{C_d}{(u-d)} + \dfrac{uC_u}{u-d} - \dfrac{dC_u}{u-d} - \dfrac{uC_u}{(u-d)} + \dfrac{uC_d}{(u-d)}) \]
\[ C = e^{-r.dt}(e^{r.dt}\dfrac{C_u}{(u-d)} - \dfrac{dC_u}{u-d} + \dfrac{uC_d}{(u-d)} - e^{r.dt}\dfrac{C_d}{(u-d)}) \]
\[ C = e^{-r.dt}(\dfrac{e^{r.dt}-d}{u-d}C_u + \dfrac{u - e^{r.dt}}{(u-d)}C_d) \]
Define \(p = \dfrac{e^{r.dt}-d}{u-d}\) and the value of the option becomes \(C = e^{-r.dt}[pC_u + (1-p)C_d]\). \(p\) can be interpreted as the risk neutral ‘up’ probability1.
An n-step binomial tree is slightly different. At the terminal (maturity) nodes payoff of the option is calculated.
\[C_{n,j} = [0,S_{n,j} - K]^+\]
where \(S_{n,j}\) is the final asset price at the state \(j\). Then each step back is calculated as the discounted expectation value of their next steps.
\[C_{i,j} = e^{-r.dt}(pC_{i+1,j+1} + (1-p)C_{i+1,j})\]
Repeat until the initial node and the initial node gives the price of the European Call option. R code is given below.
binomial_tree_EOPT<-function(S_0,K,r,T_in_days,vol,n,callput="call",return_tree=FALSE){
#Binomial tree code for the European Call
#n is the number of steps
#dt is the step size
dt <- (T_in_days/252)/n
#calculate u and d
u_val <- exp(vol*sqrt(dt))
d_val <- 1/u_val
#calculate p_u
p_star <- (exp(r*dt) - d_val) / (u_val - d_val)
#Create option value matrix
V_mat <- matrix(ncol=n+1,nrow=n+1)
#Calculate payoffs
if(callput == "call"){
V_mat[n+1,] <- pmax(S_0*d_val^n * cumprod(c(1,rep(u_val^2,n))) - K,0)
}else{
V_mat[n+1,] <- pmax(K - S_0*d_val^n * cumprod(c(1,rep(u_val^2,n))),0)
}
#Calculate discounted expectations in recursion
for(i in n:1){
V_mat[i,1:i] <- exp(-r*dt)*(p_star*V_mat[i+1,2:(i+1)]+(1-p_star)*V_mat[i+1,-((i+1):(n+1))])
}
if(return_tree){
#Return option value tree if return tree is TRUE
return(V_mat)
}else{
#Return option price if return tree is FALSE
return(V_mat[1,1])
}
}
#Price of the European Call
binomial_tree_EOPT(S_0=100,K=100,r=0.06,T_in_days=252,vol=0.166,n=3,callput="call",return_tree=FALSE)## [1] 10.18245#Price of the European Put
binomial_tree_EOPT(S_0=100,K=100,r=0.06,T_in_days=252,vol=0.166,n=3,callput="put",return_tree=FALSE)## [1] 4.358908American options are different than European options in a single but significant way: It is also possible to exercise the option and get the payoff anytime before expiration whereas it is possible to exercise only at the expiration for European options. While index options are European, most (if not all) individual stock options are American type options. If the asset is not paying dividends, value of the American Call is equivalent to the value of the European Call. Briefly it is never advisable to exercise early for American Calls. 2. For American Put, it might be.
Payoff function the of put option is for European: \(P_{European} = [0,K-S_T]^+\) where \(T\) is the expiration time. For American Put, the payoff function is \(P_{American} = [0,K-S_t]^+\) where \(0 < t < T\).
Binomial setting for the American Put is slightly different than the European part.
binomial_tree_AP<-function(S_0,K,r,T_in_days,vol,n,return_tree=FALSE){
#Binomial tree code for the American Put
#n is the number of steps
#dt is the step size
dt <- (T_in_days/252)/n
#calculate u and d
u_val <- exp(vol*sqrt(dt))
d_val <- 1/u_val
#calculate p_u
p_star <- (exp(r*dt) - d_val) / (u_val - d_val)
#Create asset price matrix
S_mat <- matrix(ncol=n+1,nrow=n+1)
S_mat[1,1] <- S_0
#Calculate terminal stock values
S_mat[n+1,] <- S_0*d_val^n * cumprod(c(1,rep(u_val^2,n)))
#Create option value matrix
V_mat <- matrix(ncol=n+1,nrow=n+1)
#Calculate payoffs
V_mat[n+1,] <- pmax(K - S_mat[n+1,],0)
#Calculate discounted expectations in recursion
for(i in n:1){
S_mat[i,1:i] <- S_mat[i+1,1:i]*u_val
V_mat[i,1:i] <- pmax(exp(-r*dt)*(p_star*V_mat[i+1,2:(i+1)]+(1-p_star)*V_mat[i+1,-((i+1):(n+1))]),K-S_mat[i,1:i])
}
if(return_tree){
#Return option value tree if return tree is TRUE
return(V_mat)
}else{
#Return option price if return tree is FALSE
return(V_mat[1,1])
}
}
binomial_tree_AP(S_0=100,K=100,r=0.06,T_in_days=252,vol=0.166,n=3,return_tree=FALSE)## [1] 4.692452Values \(u\) and \(d\) can be defined with either Cox-Ross-Rubinstein (CRR) or Jarrow and Rudd (JR) methods.
CRR method
\[u = e^{(r-1/2\sigma^2)\Delta t + \sigma \sqrt{\Delta t}}\] \[d = e^{(r-1/2\sigma^2)\Delta t - \sigma \sqrt{\Delta t}}\] \[p = 1/2\]
JR method
\[u = 1/d = e^{\sigma \sqrt{\Delta t}}\] \[p = 1/2 + \dfrac{r-1/2\sigma^2}{2\sigma}\sqrt{\Delta t }\]
These methods will only work if time step \(dt\) is small. To work on larger step sizes, we need another way.
Under Black-Scholes assumptions asset price movement is governed by Geometric Brownian Motion (GBM). GBM process is as follows.
\[dS = rSdt + \sigma S dz\]
Under GBM log of the asset price is normally distributed. Define \(x = ln(S)\) and apply Ito’s formula.
\[dx = \dfrac{1}{S}dS -\dfrac{1}/{S^2}dS^2 \]
Replace \(dS\)
\[dx = \dfrac{1}{S}(rSdt + \sigma S dz) -\dfrac{1}{S^2}(\sigma^2 S^2 dt) \]
\[dx = (r - 1/2\sigma^2)dt + \sigma dz \]
Define \(v = r - 1/2\sigma^2\) and \(dx = v dt + \sigma dx\). Differently from the multiplicative binomial model, at the next stage of \(x\) it can either be \(x + \Delta x_u\) or \(x + \Delta x_d\) where \(x_u\) and \(x_d\) are defined as the up and down step sizes. The mean and variance values of \(\Delta x\) can be calculated as follows.
\[E[\Delta x] = p\Delta x_u + (1-p)\Delta x_d = v\Delta t\] \[E[\Delta x^2] = p\Delta x_u^2 + (1-p)\Delta x_d^2 = \sigma^2\Delta t + v^2\Delta t^2\]
As with CRR and JR methods we can define equal probabilities (\(p=1/2\)) or equal jump sizes (\(\Delta x_u = - \Delta x_d\)). At equal probabilities we get the following.
\[\Delta x_u = \dfrac{1}{2}v\Delta t + \dfrac{1}{2}\sqrt{4\sigma^2 \Delta t - 3v^2 \Delta t^2}\] \[\Delta x_d = \dfrac{1}{2}v\Delta t - \dfrac{1}{2}\sqrt{4\sigma^2 \Delta t - 3v^2 \Delta t^2}\]
At equal jump sizes we get the following.
\[\Delta x = \sqrt{\sigma^2\Delta t + v^2 \Delta t^2}\] \[p = \dfrac{1}{2}(1 + \dfrac{v\Delta t}{\Delta x})\]
Reportedly equal jump sizes for the additive model is more accurate than the other three.3
#Binomial General Additive European Call
additiveBinomialEC<-function(K=100,T=1,S=100,sig=0.2,r=0.06,N=200){
#Equal jump size version
dt=T/N
nu=r-0.5*sig^2
dxu=sqrt(sig^2*dt+(nu*dt)^2)
dxd=-dxu
pu=0.5+0.5*(nu*dt/dxu)
disc=exp(-r*dt)
dpu=disc*pu
dpd=disc*(1-pu)
edxud=exp(dxu-dxd)
edxd=exp(dxd)
St<-S*exp(N*dxd)
for(j in 2:(N+1)){
St[j]<-St[j-1]*edxud
}
C<-0
for(j in 1:(N+1)){
C[j]=max(0,St[j]-K)
}
for(i in (N-1):0){
for(j in 0:i){
C[j+1]=dpu*C[j+2]+dpd*C[j+1]
St[j+1]=St[j+1]/edxd
}
}
return(C[1])
}#Binomial General Additive American Put
additiveBinomialAP<-function(K=100,T=1,S=100,sig=0.2,r=0.06,N=200){
#Equal jump size version
dt=T/N
nu=r-0.5*sig^2
dxu=sqrt(sig^2*dt+(nu*dt)^2)
dxd=-dxu
pu=0.5+0.5*(nu*dt/dxu)
disc=exp(-r*dt)
dpu=disc*pu
dpd=disc*(1-pu)
edxud=exp(dxu-dxd)
edxd=exp(dxd)
St<-S*exp(N*dxd)
for(j in 2:(N+1)){
St[j]<-St[j-1]*edxud
}
C<-0
for(j in 1:(N+1)){
C[j]=max(0,K-St[j])
}
for(i in (N-1):0){
for(j in 0:i){
C[j+1]=dpu*C[j+2]+dpd*C[j+1]
St[j+1]=St[j+1]/edxd
C[j+1]=max(C[j+1],K-St[j+1])
}
}
return(C[1])
}Hedging is as important as calculating the price of options. There are different types of hedges, also called as Greeks. Most common and basic hedge is Delta (\(\Delta\)). Delta hedge is simply the rate of change in the option price with respect to the change in the underlying asset price. Gamma (\(\Gamma\)) is the rate of change in the Delta with respect to the change in the underlying asset price. It is a second order Greek.
\[\Delta = (C_u - C_d)/(S_u-S_d)\] \[\Gamma = (\Delta_u - \Delta_d)/(S_u-S_d)\]
BinomialTreeEOPTDeltaGamma<-function(S_0=100,K=100,r=0.06,T_in_days=252,vol=0.2,n=100,callput="call"){
#Binomial tree code for the European Call
#n is the number of steps
#dt is the step size
dt <- (T_in_days/252)/n
#calculate u and d
u_val <- exp(vol*sqrt(dt))
d_val <- 1/u_val
#calculate p_u
p_star <- (exp(r*dt) - d_val) / (u_val - d_val)
#Create option value matrix
S_mat <- matrix(ncol=n+1,nrow=n+1)
S_mat[1,1] <- S_0
for(i in 2:(n+1)){
S_mat[i,1:i] <- S_0*d_val^(i-1) * cumprod(c(1,rep(u_val^2,(i-1))))
}
V_mat <- matrix(ncol=n+1,nrow=n+1)
#Calculate payoffs
if(callput == "call"){
V_mat[n+1,] <- pmax(S_mat[n+1,] - K,0)
}else{
V_mat[n+1,] <- pmax(K - S_mat[n+1,],0)
}
#Calculate discounted expectations in recursion
for(i in n:1){
V_mat[i,1:i] <- exp(-r*dt)*(p_star*V_mat[i+1,2:(i+1)]+(1-p_star)*V_mat[i+1,-((i+1):(n+1))])
}
#Delta Matrix
D_mat <- matrix(ncol=n+1,nrow=n+1)
D_mat[n+1,] <- ifelse(V_mat[n+1,]>0,1,0)
for(i in 1:n){
D_mat[i,1:i]<-(V_mat[i+1,2:(i+1)]-V_mat[i+1,1:i])/((S_mat[i+1,2:(i+1)]-S_mat[i+1,1:i]))
}
#Calculation of Gamma is different than lecture notes
G_mat <- matrix(ncol=n+1,nrow=n+1)
for(i in 1:(n-1)){
G_mat[i,1:i]<-(D_mat[i+1,2:(i+1)]-D_mat[i+1,1:i])/((S_mat[i+1,2:(i+1)]-S_mat[i+1,1:i]))
}
return(list(Price=V_mat[1,1],Delta=D_mat[1,1],Stock_tree=S_mat,Value_tree=V_mat,Delta_tree=D_mat,Gamma_tree=G_mat))
}There are options (or contract combinations) which depend on the outcome of two or more different (but correlated) assets. We can use a multidimensional binomial tree to price them.
\[dS_1(t) = (r - \delta_1)S_1(t)dt + \sigma_1S_1(t)dz_1\] \[dS_2(t) = (r - \delta_2)S_2(t)dt + \sigma_2S_2(t)dz_2\]
where the correlation between the asssets is denoted as \(dz_1dz_2=\rho dt\) and \(\delta\) is the continuous dividend yield. Check the lecture notes for the mathematical derivations of the code below.
American Spread Call is a type of option which gives the owner the right to earn the difference between two asset prices above a predefined strike price at any time till the expiration.
americanSpreadCallBinomial<-function(S1=100,S2=100,K=1,T=1,sigma1=0.2,sigma2=0.3,rho=0.5,r=0.06,N=3,div1=0.03,div2=0.04){
dt <- T/N
mu1<-r - div1 - 0.5*sigma1^2
mu2<-r - div2 - 0.5*sigma2^2
dx1<-sigma1*sqrt(dt)
dx2<-sigma2*sqrt(dt)
disc<-exp(-r*dt)
puu<-(dx1*dx2+(dx2*mu1+dx1*mu2+rho*sigma1*sigma2)*dt)/(4*dx1*dx2)*disc
pud<-(dx1*dx2+(dx2*mu1-dx1*mu2-rho*sigma1*sigma2)*dt)/(4*dx1*dx2)*disc
pdu<-(dx1*dx2+(-dx2*mu1+dx1*mu2-rho*sigma1*sigma2)*dt)/(4*dx1*dx2)*disc
pdd<-(dx1*dx2+(-dx2*mu1-dx1*mu2+rho*sigma1*sigma2)*dt)/(4*dx1*dx2)*disc
edx1<-exp(dx1)
edx2<-exp(dx2)
S1t<-S1*exp(-N*dx1)
S2t<-S2*exp(-N*dx2)
for(j in 2:(2*N+1)){
S1t[j]=S1t[j-1]*edx1
S2t[j]=S2t[j-1]*edx2
}
C<-matrix(0,ncol=2*N+1,nrow=2*N+1)
for(j in seq(-N,N,by=2)){
for(k in seq(-N,N,by=2)){
C[j+N+1,k+N+1] <- max(0,S1t[j+N+1]-S2t[k+N+1] - K)
}
}
for(i in (N-1):0){
for(j in seq(-i,i,by=2)){
for(k in seq(-i,i,by=2)){
C[j+N+1,k+N+1]<-
max(pdd*C[j+N,k+N]+pud*C[j+N+2,k+N]+pdu*C[j+N,k+N+2]+puu*C[j+N+2,k+N+2],
S1t[j+N+1]- S2t[k+N+1]-K)
}
}
}
return(C[N+1,N+1])
}Binomial trees are both practical and popular. Then, why stop there? Let’s build a trinomial tree.4
In additional to the probabilities of binomial tree (\(p_u,p_d=1-p_u\)), we now have a probability of staying at the same level (\(p_m\)). So, \(p_m + p_u + p_d = 1\). Changes in the equations are as follows. Choice of \(\Delta x\) is recommended to be \(\sigma\sqrt{3\Delta t}\).
\[E[\Delta x] = p_u(\Delta x) + p_m 0 + p_d(-\Delta x) = \nu \Delta t\] \[E[\Delta x^2] = p_u(\Delta x^2) + p_m 0 + p_d(-\Delta x^2) = \sigma^2\Delta t + \nu^2\Delta t^2\]
The system of equations above will result in the following properties.
\[p_u = \dfrac{1}{2}(\dfrac{\sigma^2\Delta t + \nu^2\Delta t^2}{\Delta x^2}+\dfrac{\nu\Delta t}{\Delta x})\] \[p_d = \dfrac{1}{2}(\dfrac{\sigma^2\Delta t + \nu^2\Delta t^2}{\Delta x^2}-\dfrac{\nu\Delta t}{\Delta x})\] \[p_m = 1 - p_u - p_d = \dfrac{\sigma^2\Delta t + \nu^2\Delta t^2}{\Delta x^2}\]
Rest is similar to Binomial pricing method. R code is as follows.
##Trinomial Trees European Call
trinomialEC<-function(S0=100,K=100,T=1,sig=0.2,r=0.06,div=0.03,N=3){
dt <- T/N
dx <- sig*sqrt(3*dt)
nu <- r - div - 0.5*sig^2
edx <- exp(dx)
pu <- 0.5*((sig^2*dt+nu^2*dt^2)/(dx^2)+nu*dt/dx)
pd <- 0.5*((sig^2*dt+nu^2*dt^2)/(dx^2)-nu*dt/dx)
pm <- 1 - pu - pd
disc <- exp(-r*dt)
S_mat <- matrix(ncol=2*N+1,nrow=N+1)
S_mat[1,N+1]<-S0
for(i in 1:N){
S_mat[i+1,(N+1-i):(N+1+i)]<- S0*(exp(((-i:i)*dx)))
}
V_mat <- matrix(ncol=2*N+1,nrow=N+1)
V_mat[N+1,] <- pmax(S_mat[N+1,]-K,0)
for(i in N:1){
V_mat[i,(N-i+2):(i+N)]<-disc*(pu*V_mat[i+1,(N-i+3):(i+N+1)] + pm*V_mat[i+1,(N-i+2):(i+N)] + pd*V_mat[i+1,(N-i+1):(i+N-1)])
}
return(list(Price=V_mat[1,N+1],V_mat=V_mat,S_mat=S_mat))
}Binomial and trinomial trees allow for 1 additional state at each time step. For instance, in a 3-step binomial tree there are 4 final states of option prices. Therefore, in order to increase the accuracy of the method there should be more time steps and decreased \(\Delta t\) so we have more states of option prices. If we can free the number of price states from the time steps we can get more granular therefore more accurate results. Suppose number of state jumps are \(N_j\).
But upper and lower boundaries of option values cannot be computed therefore boundary conditions should be defined. For European Call boundary conditions can be defined as follows.
\[\dfrac{\partial C}{\partial S} = 1\] for large \(S\). \[\dfrac{\partial C}{\partial S} = 0\] for small \(S\).
Discretize the values to grid and the result is as follows.
\[\dfrac{C_{i,N_j}-C_{i,N_j-1}}{S_{i,N_j}-S_{i,N_j-1}} = 1\] for large \(S\). \[\dfrac{C_{i,N_j}-C_{i,N_j-1}}{S_{i,N_j}-S_{i,N_j-1}} = 0\] for small \(S\).
Main idea behind finite difference methods is to discretize Black-Scholes PDE5.
\[-\dfrac{\partial C}{\partial t} = \dfrac{1}{2}S^2\sigma^2\dfrac{\partial^2 C}{\partial S^2} + (r - \delta)S\dfrac{\partial C}{\partial S} - rC\]
We can simplify it by using \(x = ln(S)\).
\[-\dfrac{\partial C}{\partial t} = \dfrac{1}{2}\sigma^2\dfrac{\partial^2 C}{\partial x^2} + \nu\dfrac{\partial C}{\partial x} - rC\]
This equation can be discretized as follows.
\[-\dfrac{C_{i+1,j}-C_{i,j}}{\Delta t} = \dfrac{1}{2}\sigma^2\dfrac{C_{i+1,j+1}-2C_{i+1,j}+C_{i+1,j-1}}{\Delta x^2} + \nu\dfrac{C_{i+1,j+1}-C_{i+1,j-1}}{2\Delta x}-rC_{i+1,j}\]
which we can rewrite as
\[C_{i,j} = p_uC_{i+1,j+1} + p_mC_{i+1,j} + p_dC_{i+1,j-1}\]
\[p_u = \Delta t (\dfrac{\sigma^2}{2\Delta x^2}+\dfrac{\nu}{2\Delta x})\] \[p_m = 1 - \Delta t \dfrac{\sigma^2}{2\Delta x^2} - r\Delta t\] \[p_d = \Delta t (\dfrac{\sigma^2}{2\Delta x^2}-\dfrac{\nu}{2\Delta x})\]
Hence the discount parameter at \(p_m\).
##Explicit Finite Difference Method European Call
explicitFDM_EC<-function(S0=100,K=100,T=1,sig=0.2,r=0.06,div=0.03,N=3,Nj=3,dx=0){
if(Nj < N){
print("Nj cannot be smaller than N. Equating Nj to N.")
Nj <- N
}
dt <- T/N
if(dx == 0){
dx <- sig*sqrt(3*dt)
}
nu <- r - div - 0.5*sig^2
edx <- exp(dx)
pu <- 0.5*dt*(sig^2/dx^2+nu/dx)
pd <- 0.5*dt*(sig^2/dx^2-nu/dx)
pm <- 1 - pu - pd - r*dt
St<-S0*exp((-Nj:Nj)*dx)
V_mat <- matrix(ncol=2*Nj+1,nrow=N+1)
V_mat[N+1,] <- pmax(0,St - K)
for(i in N:1){
V_mat[i,2:(2*Nj)]<-pu*V_mat[i+1,3:(2*Nj+1)] + pm*V_mat[i+1,2:(2*Nj)] + pd*V_mat[i+1,1:(2*Nj-1)]
V_mat[i,1]<- V_mat[i,2]
V_mat[i,(2*Nj+1)]<-V_mat[i,(2*Nj)] + (St[(2*Nj+1)]-St[(2*Nj)])
}
return(list(Price=V_mat[1,Nj+1],V_mat=V_mat))
}Here is the code for American Put.
##Explicit Finite Difference Method American Put
explicitFDM_AP<-function(S0=100,K=100,T=1,sig=0.2,r=0.06,div=0.03,N=3,Nj=3,dx=0){
dt <- T/N
if(dx == 0){
dx <- sig*sqrt(3*dt)
}
nu <- r - div - 0.5*sig^2
edx <- exp(dx)
pu <- 0.5*dt*(sig^2/dx^2+nu/dx)
pd <- 0.5*dt*(sig^2/dx^2-nu/dx)
pm <- 1 - pu - pd - r*dt
St<-S0*exp((-Nj:Nj)*dx)
V_mat <- matrix(ncol=2*Nj+1,nrow=N+1)
V_mat[N+1,] <- pmax(0, K - St)
for(i in N:1){
V_mat[i,2:(2*Nj)]<-pu*V_mat[i+1,3:(2*Nj+1)] + pm*V_mat[i+1,2:(2*Nj)] + pd*V_mat[i+1,1:(2*Nj-1)]
V_mat[i,1]<- V_mat[i,2]
V_mat[i,(2*Nj+1)]<-V_mat[i,(2*Nj)] + (St[(2*Nj+1)]-St[(2*Nj)])
V_mat[i,]<-pmax(V_mat[i,],K-St)
}
return(list(Price=V_mat[1,Nj+1],V_mat=V_mat))
}Some notes about the explicit FDM.
Implicit FDM is the more stable version of the explicit FDM with additional computational complexity. It discretizes Black-Scholes PDE with space derivatives with time step \(i\) instead of \(i+1\).
\[-\dfrac{C_{i+1,j}-C_{i,j}}{\Delta t} = \dfrac{1}{2}\sigma^2\dfrac{C_{i,j+1}-2C_{i,j}+C_{i,j-1}}{\Delta x^2} + \nu\dfrac{C_{i,j+1}-C_{i,j-1}}{2\Delta x}-rC_{i,j}\]
which we can rewrite as
\[C_{i+1,j} = p_uC_{i,j+1} + p_mC_{i,j} + p_dC_{i,j-1}\]
\[p_u = - \Delta t (\dfrac{\sigma^2}{2\Delta x^2}+\dfrac{\nu}{2\Delta x})\] \[p_m = 1 + \Delta t \dfrac{\sigma^2}{2\Delta x^2} + r\Delta t\] \[p_d = - \Delta t (\dfrac{\sigma^2}{2\Delta x^2}-\dfrac{\nu}{2\Delta x})\]
This time option values cannot be calculated independently, therefore should be solved as a system of equations with boundary conditions.
\[C_{i,N_j}-C_{i,N_j-1} = \lambda_U\] for large \(S\). \[C_{i,-N_j+1}-C_{i,-N_j} = \lambda_L\] for small \(S\).
where \(\lambda_U\) and \(\lambda_L\) depends on the type of the function. For the European Call it is \(\lambda_U = S_{i,N_j}-S_{i,N_j-1}\) and \(\lambda_L = 0\).
The set of equations is called tridiagonal.
\[ \begin{vmatrix} 1 & -1 & 0 & \cdots & \cdots & \cdots & 0 \\ p_u & p_m & p_d & 0 & \cdots & \cdots & 0 \\ 0 & p_u & p_m & p_d & 0 & \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & \cdots & 0 & p_u & p_m & p_d & 0 \\ 0 & \cdots & \cdots & 0 & p_u & p_m & p_d \\ 0 & \cdots & \cdots & \cdots & 0 & 1 & -1 \\ \end{vmatrix} \begin{vmatrix} C_{i,N_j} \\ C_{i,N_{j}-1} \\ C_{i,N_{j}-2} \\ \cdots \\ C_{i,-N_{j}+2} \\ C_{i,-N_{j}+1} \\ C_{i,-N_{j}} \\ \end{vmatrix} = \begin{vmatrix} \lambda_U \\ C_{i+1,N_{j}-1} \\ C_{i+1,N_{j}-2} \\ \cdots \\ C_{i+1,-N_{j}+2} \\ C_{i+1,-N_{j}+1} \\ \lambda_L \\ \end{vmatrix} \]
##Explicit Finite Difference Method American Put
implicitFDM_AP<-function(S0=100,K=100,T=1,sig=0.2,r=0.06,div=0.03,N=3,Nj=3,dx=0.2){
dt <- T/N
if(dx == 0){
dx <- sig*sqrt(3*dt)
}
nu <- r - div - 0.5*sig^2
edx <- exp(dx)
pu <- -0.5*dt*(sig^2/dx^2+nu/dx)
pd <- -0.5*dt*(sig^2/dx^2-nu/dx)
pm <- 1 + dt*(sig^2/dx^2) + r*dt
St<-S0*exp((-Nj:Nj)*dx)
V_mat <- matrix(ncol=2*Nj+1,nrow=N+1)
V_mat[N+1,] <- pmax(0, K - St)
lambdaL <- -1*(St[2] - St[1])
lambdaU <- 0
amat<-matrix(0,ncol=2*Nj+1,nrow=2*Nj+1)
amat[1,]<-c(1,-1,rep(0,2*Nj-1))
for(i in 2:(2*Nj)){
amat[i,]<-c(rep(0,i-2),pu,pm,pd,rep(0,2*Nj-i))
}
amat[2*Nj+1,]<-c(rep(0,2*Nj-1),1,-1)
solve(amat,vmp)
for(i in N:1){
vmp<-rev(c(lambdaL,unlist(V_mat[i+1,-c(1,2*Nj+1)]),lambdaU))
V_mat[i,]<-pmax(rev(solve(amat,vmp)),K-St)
}
return(list(Price=V_mat[1,Nj+1],V_mat=V_mat))
}For FDM methods on trinomial trees the Delta and Gamma of the option can be approximated as follows.
\[\Delta = \dfrac{\partial C}{\partial S} \approx \dfrac{C_{0,j+1} - C_{0,j-1}}{S_{0,j+1} - S_{0,j-1}}\]
\[\Gamma = \dfrac{\partial^2 C}{\partial S^2} \approx \dfrac{\left(\dfrac{C_{0,j+1} - C_{0,j-1}}{S_{0,j+1} - S_{0,j-1}}\right)-\left(\dfrac{C_{0,j+1} - C_{0,j-1}}{S_{0,j+1} - S_{0,j-1}}\right)}{1/2(S_{0,j+1}-S_{0,j-1})}\]
Risk-free rate \(r\) should be choosen such that \(d < e^{r.dt} < u\) to avoid arbitrage admission.↩
Check Merton’s ‘Theory of Rational Option Pricing’ for the proof.↩
See the textbook (Clewlow and Stickland)↩
For those who are curious about quadrinomial tree; yes there is research about that too. See http://www.tandfonline.com/doi/abs/10.1080/13504860701596745↩
Partial differential equation.↩