# Unable to obtain correct Finite Difference Results

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

$$\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$$

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)?

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.

• Thanks a lot for your reply. Previously, my code stored the FD matrix as a whole ($space_1 \times space_2 \times time$), which allowed for small grids up to 12x12 in space direction due to Memory constraints. Changed it to only save previous time matrix. Unfortunately, problem described above still exits when using 60x60 grid. I will implement a change of variable and update my post when results are in. With regards to the boundary conditions, I think they are ok, but will add the ones I used to my original post for anybody to look at.
– Pim
Jun 6, 2018 at 10:19
• I have a question regarding your proposed variable transformation though. $F_0$ can be observed as of today, but as I understood, $\alpha_0$ is a parameter that is being calibrated. Hence, how would I know its value in advance to use it for transformation to $y$?
– Pim
Jun 6, 2018 at 12:12
• The $/\alpha_0$ is here only to make $y$ symmetrical around zero so that a discretization on $[-b_y, +b_y]$ is appropriate. But you can perfectly use $y = \ln(x)$ and a discretization on $[-b^{\text{min}}_y, +b^{\text{max}}_y]$ with appropriate $b^{\text{min}}_y$ and $b^{\text{max}}_y$. For instance since after calibration $\alpha_0$ will be of the same order of magnitude as ATM implied vol you could use $b^{\text{min}}_y = \text{ATM vol} - 10 \times \nu \sqrt{T}$ and $b^{\text{max}}_y = \text{ATM vol} + 10 \times \nu \sqrt{T}$. Jun 6, 2018 at 13:06
• @Antoine_Conze I implemented the variable transformation proposed in essay.utwente.nl/59228/1/scriptie_I_Khomasuridze.pdf. For a 60x60 grid, I now obtain 'nice' positive (monotonic-increasing) values. However, in the lower part of the V matrix (high part of F / transformed x domain) I observe strange behavior (non-monotonic increasing) For two examples, see updated post. I ran the program for beta=1 and nu=0, which should give equal results to Black-Scholes. For most part of the grid values are lower than BS though. Do you have any clue about the cause?
– Pim
Jun 7, 2018 at 12:00
• Follow-up for other readers: using a larger grid and applying a variable transformation solved my problem of getting negative, non-monotonic increasing values. Although results are not fully correct yet, it solved the above question. Therefore I accepted this answer.
– Pim
Jun 18, 2018 at 14:03