The Question
I am implementing a finite difference scheme for the Heston-Hull-White PDE:
\begin{align} \frac{\partial u}{\partial t} &= \frac{1}{2}s^2v\frac{\partial^2 u}{\partial s^2 } + \frac{1}{2}\sigma_1^2 v \frac{\partial^2 u}{\partial v^2} + \frac{1}{2}\sigma_2^2 \frac{\partial^2 u}{\partial r^2} \\ &+ \rho_{1 2} \sigma_1 s v \frac{\partial^2 u}{\partial s \partial v } + \rho_{1 3} \sigma_2 s \sqrt{v} \frac{\partial^2 u}{\partial s \partial r } + \rho_{2 3} \sigma_1 \sigma_2 \sqrt{v} \frac{\partial^2 u}{\partial v \partial r } \\ &+ r s \frac{\partial u}{\partial s} + \kappa (\eta - v) \frac{\partial u}{\partial v} + a(b(T-t)-r) \frac{\partial u}{\partial r} - ru \text{.} \end{align}
I have some discretization in the $s$, $v$, $r$, and $t$ dimensions:
$${0, s_1, s_2, \ldots, s_{max}}$$ $${0, v_1, v_2, \ldots, v_{max}}$$ $${-r_{max}, r_1, r_2, \ldots, r_{max}}$$ $${0, t_1, t_2, \ldots, t_{max}}$$
I then have some boundary conditions:
\begin{align} u(s, v, r, t) &= 0 \hspace{2em} \text{ whenever } \hspace{2em} s = 0,\\ \frac{\partial u}{\partial s} (s, v, r, t) &= 1 \hspace{2em} \text{ whenever } \hspace{2em} s = s_{max},\\ u(s, v, r, t) &= s \hspace{2em} \text{ whenever } \hspace{2em} v = v_{max}, \\ \frac{\partial u}{\partial r} (s, v, r, t) &= 0 \hspace{2em} \text{ whenever} \hspace{2em} r = \pm r_{max} \end{align}
I am now left with two questions:
- How do I approximate the $\frac{\partial^2 u}{\partial r^2}$ at $r = \pm r_{max}$?
- How do I approximate the $\frac{\partial^2 u}{\partial s^2}$ at $s = s_{max}$?
A Proposed Solution
The paper I am reading is here https://arxiv.org/pdf/1111.4087.pdf . The solution the authors provide (top of page 6) is to use the boundary conditions on the first derivatives above to introduce virtual points. At $s=s_{max}$, for example, the virtual point is based on $u(s_{max-1})$ and $\frac{\partial u}{\partial s}$ at $s_{max}$.
2 - At the boundary $s=s_{max}$, approximate $\frac{\partial^2 u}{\partial s^2}$ using the central FD scheme with the "virtual point" $\Big(2s_{max} - s_{max-1}, u(s_{max-1}) + 2(s_{max} - s_{max-1})\Big)$.
1 - Do something similar at $r = \pm r_{max}$. Presumably, the authors mean that I should use the virtual points $$\Big(2r_{min} - r_{min+1}, u(r_{min+1})\Big) \hspace{2em} \text{ and } \hspace{2em} \Big(2r_{max} - r_{max-1}, u(r_{max-1})\Big)\text{.}$$
I believe that the authors chose this solution because it allowed them to directly apply Thomas' algorithm in their implicit scheme.
However, it appears to me that this solution gives an inaccurate result, especially in the case of 1. Suppose for example that $u$ is increasing at $r_{min}$. Then the proposed solution above would make $\frac{\partial^2 u}{\partial r^2}$ extremely positive.
Are there other reasons that the author proposed these solutions? Why not calculate the virtual points based on $u(s_{max})$ and $\frac{\partial u}{\partial s}$ at $s_{max}$ instead?
