Boundary Conditions in an Implicit Finite-Difference Call Option Solver
Summary
The document presents Python code for pricing a European call or put with an implicit finite-difference grid. It asks why the coefficient matrix is modified at its first and last rows and why boundary grid values at each time step are extrapolated from nearby interior values. The setup includes a terminal payoff, a discretized asset-price axis, and a backward time-stepping loop that solves a tridiagonal linear system.
The code applies adjustments to the matrix coefficients at the asset-price boundaries and uses linear extrapolation to fill those boundary values. However, the document contains no answer explaining the derivation or conditions under which these choices are appropriate. It also gives no convergence study or comparison with an analytical price, so the implementation alone does not establish the accuracy of the boundary treatment.
Key ideas
- An implicit finite-difference method solves a linear system at each backward time step.
- The terminal grid values are set from the option payoff at expiry.
- The code modifies the coefficient matrix at both asset-price boundaries.
- It extrapolates boundary values from nearby grid points, but gives no derivation or validation of that treatment.
Tags
Full text
# Issue in Understanding the Boundary Conditions for European Call Option in Implicit Finite Difference Method
# Issue in Understanding the Boundary Conditions for European Call Option in Implicit Finite Difference Method
I have a working Python code which prices European call option in Implicit Finite Difference setting. However, I am unable to understand the Boundary Conditions implemented on the coefficient matrix as mentioned in the code. Any help regarding the same will be greatly appreciated. Here is the code:
```
import numpy as np
import scipy.linalg
#parameters
S0 = 95
K = 60 #ATM Strike
r = 0.01
T = 0.5
vol = 0.2
Smax = 100
M = 6 # S
N = 6 # t
is_call = True
def Implicit_FDM(S0, K, r, T, vol, Smax, M, N, is_call=True):
M, N = int(M), int(N)
dt = T / float(N)
iVal = np.arange(1, M)
jVal = np.arange(N)
grid = np.zeros(shape=(M+1, N+1))
SValues = np.linspace(0, Smax, M+1)
alpha = 0.5*dt * (r *iVal - vol**2 *iVal**2)
beta = 1+dt * (r + vol**2 * iVal**2)
gamma = -0.5*dt * (r * iVal + vol**2 * iVal**2)
# Coefficient or Tridiagonal Matrix
coeff = np.diag(alpha[1:], -1) + np.diag(beta) + np.diag(gamma[:-1], 1)
# Terminal Condition
if is_call:
grid[:, -1] = np.maximum(SValues - K, 0)
else:
grid[:, -1] = np.maximum(K - SValues, 0)
# Boundary conditions
coeff[0, 0] = coeff[0, 0] + 2*alpha[0]
coeff[0, 1] = coeff[0, 1] - alpha[0]
coeff[-1, -1] = coeff[-1, -1] + 2*gamma[-1]
coeff[-1, -2] = coeff[-1, -2] - gamma[-1]
# Matrix Decomposition
P, L, U = scipy.linalg.lu(coeff)
for j in reversed(jVal):
Ux = scipy.linalg.solve(L, grid[1:-1, j+1])
grid[1:-1, j] = scipy.linalg.solve(U, Ux)
grid[0, j] = 2 * grid[1, j] - grid[2, j]
grid[-1, j] = 2 * grid[-2, j] -grid[-3, j]
return np.interp(S0,SValues,grid[:, 1])
#print the value
print(Implicit_FDM(S0, K, r, T, vol, Smax, M, N, is_call=True))
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.