Unable to obtain correct Finite Difference Results

Question

A rather general question regarding a specific problem I am facing with my Matlab implementation of the implicit FD method for this PDE:

\begin{equation} \frac{\sigma_s^2}{2}\frac{\partial^2 V}{\partial S^2} + \rho(t) \sigma_S \sigma_\alpha\frac{\partial^2 V}{\partial S \partial \alpha} + \frac{\sigma_\alpha^2}{2}\frac{\partial^2 V}{\partial \alpha^2} + \mu_s \frac{\partial V}{\partial S} + \mu_\alpha \frac{\partial V}{\partial \alpha} + \frac{\partial V}{\partial t} - r(t)V \end{equation}

When running my code, I obtain negative and non-monotonic-increasing values. However, they should be all positive and monotonic-increasing throughout the two-dimensions of the grid ($F$ and $\alpha$). I checked my boundary conditions again and again but they seem fine to me. Same holds for FD discretization of the PDE.

My question is, what could be the cause of my problem?

For example,

• Is there a possibility that FD formula and boundary conditions are correct, but that a particular choice of step sizes is ruining the method?

• How can I know if the problem is coding-related or that, for example, my FD matrix (A) is misspecified?

• Am I correct that FD matrix A is not dependent on time, hence can be calculated outside the FD loop?

A more detailed version of my question is given here

Edit June 7: Extended the code to handle larger grids. Applied variable transformation to both F and vol (alpha) similar to here formula (2.16). (Obviously, the upper bound on F/x changed accordingly to account for transformation from F to x). I now obtain 'nicer' results, in the sense that values are positive and for large part of the grid monotonic-increasing. However, when $\beta=1$ and $\nu=0$, I do not obtain similar results as BS. Also, for the lower part of the $V$ matrix (upper part F/x domain), values are no longer monotonic increasing (along $\alpha$/column-dimension), see left down corner V matrix (near $F_{max}$ and $\alpha_{min}$) and right down corner V matrix (near $F_{max}$ and $\alpha_{max}$). $A$ does seem to satisfy the stability constraint: $||A^{−1}||_{\infty} = 1 \leq 1$. I have noted that the diagonals of $A$ tend to increase rapidly suddenly (to an order of magnitude +04) and then decrease, after it increases again, etc (example). Is this how $A$ would normally behave (after variable transformation)?

Answer 1

a few pointers:

• did you use appropriate boundary conditions ?
• how did you truncate the space domain ? in particular it is common to first do a change of variable to get a better behaved $A$ matrix with all terms roughly of the same size. For instance in the "log normal" SABR case ($\beta=1$) it would be appropriate to work in $(x, y) \in ]-\infty, +\infty[ \times ]-\infty, +\infty[$ space with $x = \ln(F/F_0)$ and $y = \ln(\alpha/\alpha_0)$, and to truncate the space domain to $[-b_x, +b_x] \times [-b_y, +b_y]$ with appropriate bounds $b_x$ and $b_y$.
• how many grid points in each direction did you use ? For practical purpose a $50 \times 50$ discretization would be a minimum.

Also from a practical standpoint since the state space is 2 dimensional the matrix $A$ is not tridiagonal and thus very costly to invert. This is why ADI (alternate direction implicit for which you will find plenty of references) are preferred.