A)

Implement the Explicit Finite Difference method to price both European Call and Put options.

#Explicit Finite Difference method MODIFIED----
  # Explicit Implementation----
Explicit<- function(isCall, K, Tm,S0, r, sig, N, Nj=0, div, dx,returnGreeks=FALSE){
  # Finite Difference Method: i times, 2*i+1 final nodes
  # Precompute constants ----
  if(Nj!=0)
    returnGreeks=FALSE
  dt = Tm/N
  nu = r - div - 0.5 * sig^2
  edx = exp(dx)
  # got the constants formulas from clewlow 3.18, 3.19, 3.20
  pu = 0.5 * dt * ( (sig/dx)^2 + nu/dx )
  pm = 1.0 - dt *   (sig/dx)^2 - r*dt 
  pd = 0.5 * dt * ( (sig/dx)^2 - nu/dx)
  firstRow = 1
  firstCol = 1
  cp = ifelse(isCall, 1, -1)
  if(Nj!=0){
    r = nRows = lastRow = 2*Nj+1
    middleRow = Nj+1
    nCols = lastCol = N+1
    # Intialize asset prices  ----
    V = S = matrix(0, nrow=nRows, ncol=nCols)
    S[middleRow, firstCol] = S0
    S[lastRow,lastCol]= S0*exp(-Nj*dx)
    for(j in (lastRow-1):1){
      S[j,lastCol] = S[j+1,lastCol] * edx
    }
  }
  else{
    Nj=N
    r = nRows = lastRow = 2*Nj+1
    middleRow = s = nCols = lastCol = N+1
    V = S = matrix(0, nrow=nRows, ncol=nCols)
    # Intialize asset prices  ----
    S[middleRow, firstCol] = S0
    for (i in 1:(nCols-1)) {
      for(j in (middleRow-i+1):(middleRow+i-1)) {
        S[j-1, i+1] = S[j, i] * exp(dx)
        S[j ,  i+1] = S[j, i] 
        S[j+1, i+1] = S[j, i] * exp(-dx)
      }
    }
  }
  # Intialize option values at maturity ----
  for (j in 1:lastRow) {
    V[j, lastCol] = max( 0, cp * (S[j, lastCol]-K))
  }
  # Step backwards through the tree ----
  for (i in N:1) {
    for(j in (middleRow+Nj-1):(middleRow-Nj+1)) {
      # This inner for loop is only stepping through the 2 to rowsize-1, to avoid the boundaries
      V[j, i] = pu*V[j-1,i+1] + pm*V[j, i+1] + pd*V[j+1,i+1]
    }
    # Boundary Conditions ----
    stockTerm = ifelse(isCall, S[1, lastCol]-S[2,lastCol],        S[nRows-1,lastCol]-S[nRows,lastCol])
    # The last row contains the discounted value of V[lastRow, lastCol] and since 
    # this is zero for Call, we adopt the below method
    V[lastRow,  i] = V[lastRow-1,  i] + ifelse(isCall, 0, stockTerm)
    # Doing interpolation for Filling the first rows of each column
    V[firstRow, i] = V[firstRow+1, i] + ifelse(isCall, stockTerm, 0)
  }
  # Compute the Greeks ----
  if(returnGreeks && isCall){
    delta = (V[middleRow-1,firstCol+1]-V[middleRow+1,firstCol+1])/
      (S[middleRow-1,firstCol+1]-S[middleRow+1,firstCol+1])
    delta1 =(V[middleRow-1,firstCol+1]-V[middleRow,firstCol+1])/
      (S[middleRow-1,firstCol+1]-S[middleRow,firstCol+1])
    delta2 =(V[middleRow,firstCol+1]-V[middleRow+1,firstCol+1])/
      (S[middleRow,firstCol+1]-S[middleRow+1,firstCol+1])
    gamma = 2*(delta1-delta2)/((S[middleRow-1,firstCol+1]-S[middleRow+1,firstCol+1]))
    theta =((V[middleRow,firstCol+1]-V[middleRow,firstCol])/dt)/252
    return(list(Price=V[middleRow,firstCol],Delta=delta,Gamma=gamma,Theta=theta))
  }
  # Return the price ----
  return(V[middleRow,firstCol])
}
  #Implementation of Explicit----
print(paste("European Call Price=",Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.2, N=3, div=0.03, dx=0.2,returnGreeks = FALSE)))
## [1] "European Call Price= 8.54547955777436"
print(paste("European Put Price=",Explicit(isCall=FALSE, K=100, Tm=1, S0=100, r=0.06, sig=0.2, N=3, div=0.03, dx=0.2,returnGreeks = FALSE)))
## [1] "European Put Price= 5.62168917046175"

B)

Implement the Implicit Finite Difference method to price European Call and Put options.

#Implicit method MODIFIED----
  # Implicit Implementation----
Implicit = function(isCall, K, Tm,S0, r, sig, N, div, dx, Nj=0){
  # Implicit Finite Difference Method: i times, 2*i+1 final nodes
  # Precompute constants ----
  dt = Tm/N
  nu = r - div - 0.5 * sig^2
  edx = exp(dx)
  # got the constants formulas from clewlow 3.33,3.34,3.35
  pu = -0.5 * dt * ( (sig/dx)^2 + nu/dx )
  pm =  1.0 + dt *   (sig/dx)^2 + r*dt 
  pd = -0.5 * dt * ( (sig/dx)^2 - nu/dx)
  firstRow = 1
  if(Nj!=0){
    r = nRows = lastRow = 2*Nj+1
    middleRow = Nj+1
    nCols = lastCol = N+1
  }
  else{
    Nj=N
    r = nRows = lastRow = 2*Nj+1
    middleRow = s = nCols = lastCol = N+1
  }
  firstCol = 1
  cp = ifelse(isCall, 1, -1)
  
  # Intialize asset price, derivative price, primed probabilities  ----
  pp=pmp=V = S = matrix(0, nrow=nRows, ncol=nCols)
  S[middleRow, firstCol] = S0
  S[lastRow,lastCol]= S0*exp(-Nj*dx)
  for(j in (lastRow-1):1){
    S[j,lastCol] = S[j+1,lastCol] * edx
  }
  # Intialize option values at maturity ----
  for (j in firstRow:lastRow) {
    V[j, lastCol] = max( 0, cp * (S[j, lastCol]-K))
  }
  # Compute Derivative Boundary Conditions ----
  # From equation 3.38 and 3.39 in Clewlow
  if(isCall){ 
    lambdaU =(S[1, lastCol] - S[2, lastCol])
    lambdaL = 0
  }else{ #clewlows way
    lambdaU = 0
    lambdaL = -1 * (S[lastRow-1, lastCol] - S[lastRow,lastCol])
  }
  # Step backwards through the lattice ----
  for (i in (lastCol-1):firstCol) {
    h = solveImplicitTridiagonal(V, pu, pm, pd, lambdaL, lambdaU, i)
    pmp[,i] = h$pmp  # collect the pm prime probabilities
    pp [,i] = h$pp   # collect the p prime probabilities
    V = h$V
    # Apply Early Exercise condition ----
  }
  # Return the price ----
  return(list(Price=V[middleRow,firstCol],Probs=round(c(pu=pu, pm=pm, pd=pd),middleRow)))
}
  # Solving the tridiagonal matrix---- 
solveImplicitTridiagonal=function(V, pu, pm, pd, lambdaL, lambdaU, colI)
{
  # Initalize values ----
  firstRow = 1
  secondRow = 2
  thirdRow = 3
  lastRow = nRows = nrow(V)
  lastCol = ncol(V)
  # Substitute boundary condition at j = -Nj into j = -Nj+1 ----
  pp = pmp = numeric(nRows)
  pmp[lastRow-1] = pm + pd
  pp[lastRow-1]  = V[lastRow-1, lastCol] + pd*lambdaL
  
  # Eliminate upper diagonal ----
  for (j in (lastRow-2):(secondRow)) {
    pmp[j] = pm - pu*pd/pmp[j+1]
    pp[j] = V[j, colI+1] - pp[j+1]*pd/pmp[j+1]
  }
  # Use boundary conditions at j = Nj and equation at j=Nj-1 ----
  V[firstRow, colI] = (pp[secondRow] + pmp[secondRow]*lambdaU)/(pu + pmp[secondRow])
  V[secondRow, colI] = V[firstRow,colI] - lambdaU
  # Back-substitution ----
  for(j in thirdRow:lastRow) {
    V[j, colI] =  (pp[j] -pu*V[j-1, colI])/pmp[j]
  }
  V[lastRow, colI] = V[lastRow-1, colI] - lambdaL
  # Return values ----
  list(V=V, pmp=pmp, pp=pp)
}
  # Testing the implicit method ----
print(paste("European Call Price=",Implicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.2, N=3, div=0.03, dx=0.2)$Price))
## [1] "European Call Price= 7.56078104077051"
print(paste("European Put Price=",Implicit(isCall=FALSE, K=100, Tm=1, S0=100, r=0.06, sig=0.2, N=3, div=0.03, dx=0.2)$Price))
## [1] "European Put Price= 4.76799996252376"

C)

For both the Explicit and Implicit Finite Difference schemes estimate the numbers ∆t, ∆x as well as the total number N j of points on the space grid x to obtain a desired error of ε = 0.001.

# Question1.C----
  # Carried out implementation from clewlow pg 65---- 
QuestionC<-function(N=3,sig,Tm,nsd,error){
  repeat{
    dt=Tm/N
    Nj=(sqrt(N)*nsd)/(2*sqrt(3)) -.5
    dx=(nsd*sig*sqrt(Tm))/(2*Nj+1)
    N=N+1
    if(((dx*dx)+dt)<=error){
      print(paste("ANSWER:",dt,Nj,dx,N))
      return(list(dt,Nj,dx,N))
      break
    }
  }
}
  # Testing it for the above conditions 
c=QuestionC(N=3,sig=.2,Tm=1,nsd=6,error = .001)
## [1] "ANSWER: 0.000892857142857143 57.4655069847578 0.0103509833901353 1121"
print(paste("N=",c[4]," Nj=",c[2]," dx=",c[3]))
## [1] "N= 1121  Nj= 57.4655069847578  dx= 0.0103509833901353"

D)

Consider S 0 = 100, K = 100, T = 1 year, σ = 25%, r = 6%, δ = 0.03. Calculate and report the price for European Call and Put using both explicit and implicit FD methods and the number of steps that you calculated in the previous point (part c).

explicit_call=Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))
explicit_put=Explicit(isCall=FALSE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))

implicit_call=Implicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))$Price
implicit_put=Implicit(isCall=FALSE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))$Price
print(paste("Explicit European Call Price=",explicit_call))
## [1] "Explicit European Call Price= 11.012409566251"
print(paste("Explicit European Put Price=",explicit_put))
## [1] "Explicit European Put Price= 8.14417238983029"
print(paste("Implicit European Call Price=",implicit_call))
## [1] "Implicit European Call Price= 11.0098296631494"
print(paste("Implicit European Put Price=",implicit_put))
## [1] "Implicit European Put Price= 8.14181701536698"

E)

Repeat part (c) of this problem but this time get the empirical number of iterations. Specifically, obtain the Black Scholes price for the data in (d), then do an iterative procedure to figure out the ∆x, ∆t, N , and N j to obtain an accuracy of ε = 0.001.

QuestionE<-function(error,N=100,Tm,nsd,isAmerican, isCall, K, S0, r, sig, div, dx,isExplicit=TRUE){
  BS_price = BSM(type=ifelse(isCall, "c", "p"), K=K, t=Tm, S=S0, r=r, sigma=sig,div=div)
  print(BS_price)
  repeat{
    dt=Tm/N
    Nj=(sqrt(N)*nsd)/(2*sqrt(3)) -.5
    dx=(nsd*sig*sqrt(Tm))/(2*Nj+1)
    if(isExplicit)
      fd_price=Explicit(isCall=isCall, K=K, Tm=Tm, 
                              S0=S0, r=r, sig=sig, N=N, Nj=as.integer(Nj), div=div, dx=dx)
    else
      fd_price=Implicit(isCall=isCall, K=K, Tm=Tm, 
                              S0=S0, r=r, sig=sig, N=N, Nj=as.integer(Nj), div=div, dx=dx)$Price
    N=N+100
    if(abs(BS_price-fd_price)<=error){
      print(fd_price)
      print(paste("ANSWER:",dt,Nj,dx,N))
      return(list(dt,Nj,dx,N))
      break
    }
  }
}
BSM<-function(S, K, t, r, sigma,type,div=0){
  d1 <- (log(S/K)+(div+sigma^2/2)*t)/(sigma*sqrt(t))
  d2 <- d1 - sigma * sqrt(t)
  if (type == "c")
    result <- S*exp((div-r)*t)*pnorm(d1) - K*exp(-r*t)*pnorm(d2)
  if (type == "p")
    result <- K*exp(-r*t) * pnorm(-d2) - S*exp((div-r)*t)*pnorm(-d1)
  return(result)
}

E=QuestionE(error=0.001,Tm=1,nsd=6,isAmerican=FALSE,isCall=TRUE,isExplicit = TRUE,
            K=100, S0=100, r=0.06, sig=0.25, div=0.03)
## [1] 11.01308
## [1] 11.01212
## [1] "ANSWER: 0.000588235294117647 70.9142842854285 0.0105021006302101 1800"
print(paste("For Explicit: Nj=",as.integer(E[2])," dx=",E[3]," N=",E[4]))
## [1] "For Explicit: Nj= 70  dx= 0.0105021006302101  N= 1800"
E=QuestionE(error=0.001,Tm=1,nsd=6,isAmerican=FALSE,isCall=TRUE,isExplicit = FALSE,
            K=100, S0=100, r=0.06, sig=0.25, div=0.03)
## [1] 11.01308
## [1] 11.0121
## [1] "ANSWER: 0.000238095238095238 111.749721603218 0.00668153104781061 4300"
print(paste("For Implicit: Nj=",as.integer(E[2])," dx=",E[3]," N=",E[4]))
## [1] "For Implicit: Nj= 111  dx= 0.00668153104781061  N= 4300"

F)

Using the parameters from part (d), plot on the same graph the implicit fi- nite difference probabilities p u , p m , p d as a function of σ, σ ∈ {0.05, 0.1, 0.15, …, 0.6}. Write detailed comments on your observations.

pu = function(sig, dt=1/3, dx=0.2, r=0.06, div=0.03, nu=r-div-0.5*sig^2) {
  0.5 * dt * ( (sig/dx)^2 + nu/dx )
}

pm =  function(sig, dt=1/3, dx=0.2, r=0.06, div=0.03, nu=r-div-0.5*sig^2) {
  1.0 - dt *   (sig/dx)^2 - r*dt 
}
pd =  function(sig, dt=1/3, dx=0.2, r=0.06, div=0.03, nu=r-div-0.5*sig^2) {
  0.5 * dt * ( (sig/dx)^2 - nu/dx)
}

x = seq(from=0.01, to=0.5, by=0.01)
ux = pu(x)
mx = pm(x)
dx = pd(x)
theTitle = expression(paste("Explicit Finite Difference Probabilities", 
                            " as a Function of ", sigma, sep=""))
xlab1 = expression(paste("Volatility (", sigma, ")", sep=""))
xlab2 = expression(paste("K = 100, T = 1.0, S =100, r = 6%, ", delta, " = 3%"))
plot(x, ux, ylim=c(-1,1), type="n", xlab="")
title(theTitle, cex=3)
mtext(text=xlab1, side=1, line=2, cex=1.25)
mtext(text=xlab2, side=1, line=3, cex=1.25)
lines(x, ux, lwd=2, col=2)
lines(x, mx, lwd=2, col=3)
lines(x, dx, lwd=2, col=4)
abline(h=0)
legend(0.4, 0.5, legend=c('pu', 'pm', 'pd'), col=2:4, lty=1, lwd=3)

G)

Implement the Crank-Nicolson Finite Difference method and price both European Call and Put options. Use the same parameters as in part (d) and the same number of steps in the grid. Put the results of the 3 methods (EFD, IFD, CNFD) side by side in a table and write your observations.

#Crank Nicholson Method----
CrankNicholson = function(isAmerican, isCall, K, Tm,S0, r, sig, N, div, dx, Nj=0){
  # Crank Nicholson Finite Difference Method: i times, 2*i+1 final nodes
  # Precompute constants ----
  dt = Tm/N
  nu = r - div - 0.5 * sig^2
  edx = exp(dx)
  pu = -0.25     *dt * ( (sig/dx)^2 + nu/dx )  
  pm =  1.0 + 0.5*dt *   (sig/dx)^2 + 0.5*r*dt 
  pd = -0.25     *dt * ( (sig/dx)^2 - nu/dx)   
  firstRow = 1
  firstCol = 1
  if(Nj!=0){
    r = nRows = lastRow = 2*Nj+1
    middleRow = Nj+1
    nCols = lastCol = N+1
  }
  else{
    Nj=N
    r = nRows = lastRow = 2*Nj+1
    middleRow = s = nCols = lastCol = N+1
  }
  middleRow = nCols = lastCol = Nj+1
  
  cp = ifelse(isCall, 1, -1)
  
  # Intialize asset price, derivative price, primed probabilities  ----
  pp=pmp=V = S = matrix(0, nrow=nRows, ncol=nCols)
  S[middleRow, firstCol] = S0
  S[lastRow,lastCol]= S0*exp(-Nj*dx)
  for(j in (lastRow-1):1){
    S[j,lastCol] = S[j+1,lastCol] * edx
  }
  # Intialize option values at maturity ----
  for (j in firstRow:lastRow) {
    V[j, lastCol] = max( 0, cp * (S[j, lastCol]-K))
  }
  # Compute Derivative Boundary Conditions ----
  if(isCall){
    lambdaU =(S[1, lastCol] - S[2, lastCol])
    lambdaL = 0
  }else{
    lambdaU = 0
    lambdaL = round(-1 * (S[lastRow-1, lastCol] - S[lastRow,lastCol]),2)
  }
  # Step backwards through the lattice ----
  for (i in (lastCol-1):firstCol) {
    h = solveCrankNicholsonTridiagonal(V, pu, pm, pd, lambdaL, lambdaU, i)
    pmp[,i] = round(h$pmp,4)  # collect the pm prime probabilities
    pp [,i] = round(h$pp, 4)  # collect the p prime probabilities
    V = h$V
    # Apply Early Exercise condition for American Options ----
    for(j in lastRow:firstRow) {
      V[j, i] = max(V[j, i], cp * (S[j, lastCol] - K))
    }
  }
  # Return the price ----
  list(Type = paste(ifelse(isCall, "Call", "Put")),Price = V[middleRow,firstCol],
       Probs=round(c(pu=pu, pm=pm, pd=pd), 4), pmp=pmp, pp= pp,
       S=round(S,2), V=round(V,middleRow))
}
  #TridiagonalMatrix solver for crank nicholson----
#it was solveICrankNicholsonTridiagonal
solveCrankNicholsonTridiagonal=function(V, pu, pm, pd, lambdaL, lambdaU, colI)
{
  # Initalize values ----
  firstRow = 1
  secondRow = 2
  thirdRow = 3
  lastRow = nRows = nrow(V)
  lastCol = ncol(V)
  # Substitute boundary condition at j = -Nj into j = -Nj+1 ----
  pp = pmp = numeric(nRows)
  pmp[lastRow-1] = pm + pd
  pp[lastRow-1]  = (- pu   *V[lastRow-2, lastCol] 
                    -(pm-2)*V[lastRow-1, lastCol]
                    - pd   *V[lastRow  , lastCol] + pd*lambdaL)
  # Eliminate upper diagonal ----
  for (j in (lastRow-2):(secondRow)) {
    pmp[j] = pm - pu*pd/pmp[j+1]
    pp[j] = ( - pu   *V[j-1, colI+1] 
              -(pm-2) *V[j  , colI+1]
              - pd    *V[j+1, colI+1] 
              -pp[j+1]*pd/pmp[j+1])
  }               
  # Use boundary conditions at j = Nj and equation at j=Nj-1 ----
  V[firstRow, colI] = (pp[secondRow] + pmp[secondRow]*lambdaU)/(pu + pmp[secondRow])
  V[secondRow, colI] = V[firstRow,colI] - lambdaU
  # Back-substitution ----
  for(j in thirdRow:lastRow) {
    V[j, colI] =  (pp[j] -pu*V[j-1, colI])/pmp[j]
  }
  V[lastRow, colI] = V[lastRow-1, colI] - lambdaL
  # Return values ----
  list(V=V, pmp=pmp, pp=pp)
}



  #Implementation of Crank Nicholson----
cnd_call=CrankNicholson(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))$Price
cnd_put=CrankNicholson(isCall=FALSE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))$Price
# Tabulating all the Finite Difference Methods
expl=c(explicit_call,explicit_put)
impl=c(implicit_call,implicit_put)
cnd=c(cnd_call,cnd_put)
index=c("Call","Put")

finite_difference_df = data.frame(index,expl,impl,cnd)
colnames(finite_difference_df) = c("Type of Options", "Explicit Method", "Implicit Method", "Crank Nicholson")
print(finite_difference_df)
##   Type of Options Explicit Method Implicit Method Crank Nicholson
## 1            Call       11.012410       11.009830       11.011271
## 2             Put        8.144172        8.141817        8.509172

We can see that the Crank Nicholson method converges faster as the convergence rate of Cranknichlson O(dx^2 +(dt/2)^2) is smaller than the convergence rate of Explicit and Implicit O(dx^2 +dt). We can also see that when dx is very small for explicit, it fails to converge to the black shole price and therefore its conditionally stable whereas implicit and crank nicholson converges for a wide range of values of dx and hence we can vary dx and dt in implicit and crank nicholson to tradeoff between speed of convergence and accuracy.

H)

Calculate the hedge sensitivities for the European call option using the Explicit Finite Difference method. You need to calculate Delta, Gamma, Theta, and Vega.

sigma=.25
# explicit_call=Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))
explicit=Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=0.25, N=as.integer(c[4]) , div=0.03, dx=as.double(c[3]),returnGreeks = TRUE)
delta=explicit$Delta
gamma=explicit$Gamma
theta=explicit$Theta
dsigma = .001* sigma
Vega = (Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=sigma+dsigma, N=as.integer(c[4]), div=0.03, dx=as.double(c[3]))-Explicit(isCall=TRUE, K=100, Tm=1, S0=100, r=0.06, sig=sigma-dsigma, N=as.integer(c[4]), div=0.03, dx=as.double(c[3])))/(2*dsigma)/100


print(paste("Delta=",delta,"Gamma=",gamma,"Theta=",theta,"Vega=",Vega))
## [1] "Delta= 0.579141905833163 Gamma= 0.015033207780097 Theta= -0.0229163669908173 Vega= 0.375832253819119"