Estoy resolviendo la clásica EDP de Black & Scholes (BS) para una opción europea utilizando diferencias finitas y el esquema implícito. En otras palabras, estoy tratando de resolver
$\displaystyle\frac{\partial V}{\partial t} + \frac 12\sigma^2\frac{\partial^2V}{\partial x^2} + \mu\frac{\partial V}{\partial x} - rV = 0$
con la condición de límite $V(T,x) = \Phi(e^x) = \max\{e^x - K, 0\}$ . El parámetro de deriva es el habitual $\mu = r - y - \frac 12 \sigma^2$ .
El problema es que cuando mi $x$ es demasiado estrecha, mis precios de las opciones se desvían de la BS. En concreto, estoy resolviendo el caso ATM con S0 = K = 2775 y si uso [np.log(2675), np.log(2875)] como los puntos límite, la solución tiene un error enorme (~400). Sin embargo, cuando utilizo una malla más amplia [np.log(1775), np.log(3775)] funciona a las mil maravillas.
Pensé que el esquema implícito sería lo suficientemente estable numéricamente como para no mostrar este tipo de comportamiento, así que ¿debo asumir que hay algo mal en mi código Python?
Mis otros parámetros son r, y, T, sigma = (0.0278, 0.0189, 1, 0.15) .
def buildDs(nx):
D1 = np.diag(np.ones(nx), 1) - np.diag(np.ones(nx), -1)
D2 = np.diag(np.ones(nx), 1) + np.diag(np.ones(nx), -1) - 2 * np.eye(nx + 1)
D1[0, 2] = -1
D1[0, 1] = 4
D1[0, 0] = -3
D1[-1, -1] = 3
D1[-1, -2] = -4
D1[-1, -3] = 1
D2[0, 2] = 1
D2[0, 1] = -2
D2[0, 0] = 1
D2[-1, -1] = 1
D2[-1, -2] = -2
D2[-1, -3] = 1
return csc_matrix(D1), csc_matrix(D2)
def pxFD(Phi, r, mu, sigma,
xs, ts,
nx = 2000, nt = 10000):
dx = (xs[1] - xs[0]) / nx
dt = (ts[1] - ts[0]) / nt
xs = np.linspace(xs[0], xs[1], num = nx + 1, endpoint = True)
ts = np.linspace(ts[0], ts[1], num = nt + 1, endpoint = True)
V = np.zeros((nt + 1, nx + 1))
I = identity(nx + 1)
D1, D2 = buildDs(nx)
L = 1 / 2 * (sigma / dx) ** 2 * D2 + mu / (2 * dx) * D1 - r * I
P = I - dt * L
V[-1] = Phi(xs)
for j in reversed(range(nt)):
V[j] = spsolve(P, V[j + 1])
return V, xs, ts
def pxFDGBM(Phi, r, y, sigma,
Ss, ts,
nx = 2000, nt = 10000):
mu = r - y - 1 / 2 * sigma ** 2
xs = np.log(Ss)
Philog = lambda x: Phi(np.exp(x))
Vs, xs, ts = pxFD(Phi = Philog,
r = r,
mu = mu,
sigma = sigma,
xs = xs,
ts = ts,
nx = nx,
nt = nt)
return Vs, np.exp(xs), ts
