6

In an intermediary step to completing my undergraduate thesis, I'm attempting to solve the Black-Scholes PDE in R by an explicit finite difference scheme. As of now, I keep getting hung up on simultaneously iterating backwards in time (along the time discretization) and along the stock price discretization. I've implemented the boundary conditions for a European call option and am trying to solve my way backwards according to Phil Goddard's description for the model.

R <- matrix(NA, ncol=length(t), nrow=length(S))

# Value at Maturity:
R[,length(t)] <- Vold <- pmax(S - K,0) 

# Value at S = S_max; V = fwd
R[length(S),] <- V_max <- S_max - K * exp(-r_yr * (t_T - t))

# Value at S = S_0; V = 0
R[1,] = V_min <- t - t

# Run Model
# -----------------------------------------------------------------------------------------

# Define Coefficient-Functions a_j, b_j, and c_j
a_ <- function(j){
  a_j <- (1/2) * dt * (r*j - sigma^2 * j^2)
  return(a_j)
}

b_ <- function(j){
  b_j <- 1 + dt * (sigma^2 * j^2 + r)
  return(b_j)
}

c_ <- function(j){
  c_j <- -(1/2) * dt * (r*j + sigma^2 * j^2)
  return(c_j)
}

# Move backward in time stepping on each stock price within the BCs S = S_0, S = S_max
for (i in length(t):2){
  for (j in (length(S)-1):2){
    R[,(length(i)-1)][j] <- b_(j) * R[,length(i)][j] + c_(j) * R[,length(i)][j+1] + 
                            a_(j) * R[,length(i)][j-1]
  }
}

I've discretized time and stock price such that the user can specify the desired number of mesh points (resulting in "length(t)" and "length(S)" nodes, respectively). Since I specify 3 outer edges of matrix R (2 across time at S_max and S_0 over the columns of R; 1 across initial option value at maturity, the last column of R), I iterate over the length of the columns minus 1 and the length of the rows minus 2.
(please see link to diagram for further description)

Link to BC image, at bottom of page

When I run the model, the boundary conditions are correctly enacted, but I can't move any further back in time... Any guidance on my (likely rookie-mistake-filled) code would be greatly appreciated. Please let me know if there is anything further I can add (code or description) to clarify.

Thanks,

Jake

dunno
  • 205
  • 1
  • 5
  • length(i) should always return 1 because i is simply an integer value taken from length(t):2, right? So all you're doing is recycling the values from R[,1] to recalculate R[,0]. – Evan Johnson Oct 11 '14 at 00:05

0 Answers0