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