von Neumann boundary in the transformed PDE

Question

I have transformed the BSM PDE $$\frac{\partial V}{\partial t} + \frac{\sigma^2}{2}S^2 \frac{\partial^2 V}{\partial S^2} + rS \frac{\partial V}{\partial S} - rV = 0 $$ to $u(\tau,x) = V(T-\tau,S_{0} e^{x})$ with a change of variables $x = \ln(S/S_{0})$ and $\tau = T -t$ to $$ \frac{\partial u}{\partial \tau} = \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2} + (r - \frac{\sigma^2}{2})\frac{\partial u}{\partial x} - ru $$ For the discretization I use backward difference for the time derivative and central differences for the spatial derivatives. So I have an implicit finite difference scheme of the form $$ A u_{i}^{n+1} = u_{i}^{n} $$ where $i$ is the index for the spatial grid (from 1 to $M = 2B/\Delta x$) and $n$ for the time grid (from 1 to $N = T/\Delta\tau$) and $$ A = \begin{pmatrix} \beta_{1} & \gamma_{1} & 0 & 0 & \cdots & 0\\ \alpha_{2} & \beta_{2} & \gamma_{2} & 0 & \cdots & 0 \\ 0 & \alpha_{3} & \beta_{3} & \gamma_{3} & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots &\ddots &\vdots \\ 0 &\cdots & & & \alpha_M & \beta_{M} \end{pmatrix} $$ and $\beta_i = 1 - k (\frac{\sigma^2}{\Delta x^2} - r)$, $\gamma_i = k(- \frac{\sigma^2}{2 \Delta x^2} - (r - \frac{\sigma^2}{2}) \frac{1}{2\Delta x})$ and $\alpha_i = k( -\frac{\sigma^2}{2 \Delta x^2} + (r - \frac{\sigma^2}{2})\frac{1}{2 \Delta x})$ $\forall 2 \leq i \leq N-1$ From the boundary conditions for a European call $$\lim_{S \rightarrow 0} V(S,t) = 0 \qquad \lim_{x \rightarrow -\infty} u(x,\tau) = 0$$ I set $\beta_1 = 1$, $\gamma_{1} = 0$ and $u_{i}^{n} = 0$ as the lower bound and from $$\lim_{S \rightarrow \infty} V(S,t) = S - K e^{-r(T-t)} \qquad \lim_{x \rightarrow \infty} u(x,\tau) = S_{0} e^{x} - K e^{-r(T-t)}$$ I set $\alpha_M = 0$, $\beta_{M} = 1$ and $u_{M}^{n} = S_{0} e^{x_{M}} - K e^{-r(T-t)}$ to solve this scheme.

Now I want to use von Neumann boundaries and I'm not sure how it's done in the transformed PDE. So far I have tried it by leaving the lower boundary as it is and for the upper boundary I used $$\lim_{S \rightarrow \infty} \frac{\partial V(S,t)}{\partial S} = 1 \qquad \lim_{x \rightarrow \infty} \frac{\partial u(x,\tau)}{\partial x} = S_{0} e^{x} $$ From the central difference of the first spatial derivative I get then $$ \frac{u_{M+1}^{n+1} - u_{M-1}^{n+1}}{2 \Delta x} = e^{x_{M}} S_{0} \\ u_{M+1}^{n+1} = 2 \Delta x e^{x_{M}} S_{0} - u_{M-1}^{n+1}$$ Now I insert this into the last row of the scheme $$\alpha_{M} u_{M-1}^{n+1} + \beta_{M} u_{M} + \gamma_{M} (2 \Delta x e^{x_{M}} S_{0} - u_{M-1}^{n+1}) = u_{M}^{n} \\ (\alpha_{M} - \gamma_{M}) u_{M-1}^{n+1} + \beta_{M} u_{M} = u_{M}^{n} - 2 \Delta x e^{x_{M}} S_{0} \gamma_{M} $$ So the boundary in the scheme becomes $\alpha_{M} = (\alpha_{M} - \gamma_{M})$, $\beta_M = \beta_M$ and $u_{M}^{n} = u_{M}^{n} - 2 \Delta x e^{x_{M}} S_{0} \gamma_{M}$

When I implement this in Matlab with the Dirichlet boundary I get good results, but something must be wrong with the von Neumann implementation because the results for large $S$ have very large errors. Can someone please help me?

Thanks a lot

Answer 1

I suggest this procedure: $$x=\ln\left(\frac{S}{S_0}\right)\,\,\,\ ,\,\,\,\tau=T-t\,\,\,\, ,\,\,\,\,U(x,\tau)=e^{r\tau}V(S_t,t)$$ which leads to a forward PDE with constant coefficients $$U_\tau(x,\tau)=\frac{\sigma^2}{2}U_{xx}(x,\tau)+\mu\,U_{x}(x,\tau)$$ where $$\mu=r-\frac{1}{2}\sigma^2$$ After this change of variables, the terminal condition becomes an initial condition at $\tau = 0$: $$U(x,0)=\max\{S_0\,e^{x_0}-K,0\}$$ In order to apply standard methods of numerical solution of PDEs, such as finite difference or finite element schemes, we need to localize the domain of definition in $x$ to a bounded interval $$x\in(x_{min}\,,\,x_{max})$$ The standard approach consists in choosing the computational domain $(x_{min}\,,\,x_{max})$ sufficiently large” and imposing Dirichlet or Neumann boundary conditions based on the asymptotic of the solution. For instance, in the case of a call option, we have the following asymptotic behavior of the solution: $$V_S(S,t)\simeq\,S-K\,e^{-r(T-t)}\,\, ,\,\,S\to\infty$$ $$V_S(S,t)\simeq0\,\, ,\,\,S\to0$$ Neumann Boundaries: $$u_x(x_{max},\tau)=S_0e^{x_{max}+r\tau}-K$$ $$u_x(x_{min},\tau)=0$$

It was proved that the localization error due to these boundary conditions goes to zero exponential when the size of the computational domain goes to infinity. However, these results do not indicate whether it is better to take Dirichlet or Neumann conditions nor how to choose domain “sufficiently large” in practice. One often suggests to take this interval in relation with the standard deviation of the Brownian motion for the period $[0,T]$. This approach uses the so-called transparent boundary conditions:

• Transparent boundary conditions may be used for any interval $x\in(x_{min}\,,\,x_{max})$ even small, provided it contains the singularities of the payoff.
• They provide exact boundary conditions for the localized PDE (and not only asymptotically exact). There is still a small numerical error due to the discretization of these conditions.
• Transparent boundary conditions are almost as easy to implement as Neumann or Dirichlet conditions. For more details,you can read this

Discretization

$$x=x_{min}+i\Delta\,x\,\, ,\,\,\Delta\,x=\frac{x_{max}-x_{min}}{N-1}$$ $$\tau_n=n\Delta t\,\,\, ,\,\,\Delta t=\frac{T}{M}$$

Iimplementation

We have a usual three-point finite difference scheme: $$a_lU_{i-1}^{n+1}+a_dU_{i}^{n+1}+a_uU_{i+1}^{n+1}=b_lU_{i-1}^{n}+b_dU_{i}^{n}+b_uU_{i}^{n+1}$$

In the case of the $\theta$-scheme, the coefficients are given by \begin{align} & {{a}_{l}}=\theta \,\Delta t\left( -\frac{{{\sigma }^{2}}}{2{{(\Delta x)}^{\,2}}}+\frac{\mu }{2\Delta x} \right) \\ & {{a}_{d}}=1+\theta \Delta t\frac{{{\sigma }^{2}}}{{{(\Delta x)}^{\,2}}} \\ & {{a}_{u}}=\theta \Delta t\left( -\frac{{{\sigma }^{2}}}{2{{(\Delta x)}^{\,2}}}-\frac{\mu }{2\Delta x} \right) \\ & {{b}_{l}}=(\theta -1)\Delta t\left( -\frac{{{\sigma }^{2}}}{2{{(\Delta x)}^{\,2}}}+\frac{\mu }{2\Delta x} \right) \\ & {{b}_{d}}=1+(\theta -1)\theta \Delta t\frac{{{\sigma }^{2}}}{{{(\Delta x)}^{\,2}}} \\ & {{b}_{u}}=(\theta -1)\Delta t\left( -\frac{{{\sigma }^{2}}}{2{{(\Delta x)}^{\,2}}}-\frac{\mu }{2\Delta x} \right) \\ \end{align}

Answer 2

Just bumped into this old question and since it hasn't really been answered and for what it's worth now, the actual Neumann condition in the log coordinate is wrong. Instead of $lim_{x \rightarrow \infty} \frac{\partial u(x,\tau)}{\partial x} = S_{0} e^{x}$, it should read $lim_{x \rightarrow \infty} \frac{\partial u(x,\tau)}{\partial x} = e^{x}$. That was probably the reason for the large errors at large S.