Implicit finite difference method always guarantees positive and stable price of derivative?

Question

For the following black scholes pde $$ f_t + rSf_S+\frac{1}{2}\sigma^2S^2f_{SS} = rf $$

By denoting $f_{i}^{n} = $ Price of derivative at price node $i$ and time node $n$ and assume uniform grid, the corresponding implicit scheme would be $$ a_if_{i-1}^n + b_if_{i}^n + c_if_{i+1}^n = f_i^{n+1} $$ where $$ a_i = -\frac{1}{2}\Delta t \left( \frac{\sigma^2S_i^2}{\Delta S^2} - \frac{rS_i}{\Delta S} \right) = -\frac{1}{2}\Delta t(\sigma^2i^2 - ri)\\ b_i = 1+\Delta t \left( \frac{\sigma^2S_i^2}{\Delta S^2}+r \right) = 1+\Delta t(\sigma^2i^2 + r) \\ c_i = -\frac{1}{2}\Delta t \left( \frac{\sigma^2S_i^2}{\Delta S^2} + \frac{rS_i}{\Delta S} \right) = -\frac{1}{2}\Delta t(\sigma^2i^2 + ri) $$

In matrix form, $$ CF_n + K_n = F_{n+1} \\ F_n = C^{-1}\left( F_{n+1}-K_n \right) $$ where $$ F_n= \begin{pmatrix} f_1^n \\ f_2^n \\ \vdots \\ f_{M-1}^{n} \end{pmatrix}\\ C = \begin{pmatrix} b_1 & c_1 & 0 & \cdots & 0 & 0 \\ a_2 & b_2 & c_2 & \cdots & 0 & 0 \\ 0 & a_3 & b_3 & \cdots & 0 & 0 \\ 0 & \vdots &\vdots &\ddots &\vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{M-1} & b_{M-1} \end{pmatrix} $$ $$ K_n = \begin{pmatrix} a_1f_0^n \\ 0 \\ \vdots \\ 0 \\ c_{n-1}f_M^n \end{pmatrix} $$ where $f_0$ and $f_M$ are two ends of the price grid with some boundary conditions.

There are two questions to ask

1. All the cofficients should be greater than or equal to zero to guarentee that the pricing of derivative is alwasy positive, because referecnece I have read so far mentions that for explicit scheme the cofficients must be greater than equal to zero but not for implicit scheme. I guess it is not necessary, since $a_i \geq 0$ when $$ \frac{\Delta S}{S_i} \geq \frac{\sigma^2}{r} $$ and this would hold for small $S_i$.
2. For stability, I think that $\left\|C\right\|_{\infty} \geq 1$ as we take inverse of $C$. When $a_i < 0$ and $c_i \geq 0$, $$ \begin{align} |a_i|+|b_i|+|c_i| &= \frac{1}{2}\Delta t(\sigma^2i^2 - ri) + 1 + \Delta t(\sigma^2i^2 + r) - \frac{1}{2}\Delta t(\sigma^2i^2 + ri) \\ &= -\Delta t \frac{rS_i}{\Delta S} + 1 + \Delta t\frac{\sigma^2S_{i}^{2}}{\Delta S^2} + \Delta r \end{align} $$ and it should be greater than or equal to 1. $$ \begin{align} & -\Delta t \frac{rS_i}{\Delta S} + 1 + \Delta t\frac{\sigma^2S_{i}^{2}}{\Delta S^2} + \Delta r \geq 1 \\ \implies & -\frac{rS_i}{\Delta S} + \frac{\sigma^2S_{i}^{2}}{\Delta S^2}+r \geq 0 \\ \implies & -rS_i\Delta S + \sigma^2S_{i}^{2} + \Delta S^2r \geq 0 \end{align} $$ By letting $g(S_i, \Delta S) = -rS_i\Delta S + \sigma^2S_{i}^{2} + \Delta S^2r$, it requires minimum of $g$ greater than or equal to 0. $$ g_{S_i} = -r\Delta S + 2\sigma^2S_i = 0 \implies S_{i}^{*} = \frac{r\Delta S}{2\sigma^2} $$ and $$ \begin{align} g(S_{i}^{*},\Delta S) &= -\frac{r^2\Delta S^2}{2\sigma^2} + \frac{r^2\Delta S^2}{4\sigma^2} + r\Delta S^2 \\ &= -\frac{2r}{4\sigma^2} + \frac{r}{4\sigma^2} + 1 \\ &= -\frac{r}{4\sigma^2} + 1 \geq 0 \\ \implies & \frac{\sigma^2}{r} \geq \frac{1}{4} \end{align} $$ Therefore, I think that the iteration is not stable for $\frac{\sigma^2}{r} < \frac{1}{4}$.

I have tried to find references, but most of them used change of variables to transform black scholes pde into normal heat equation and used von-neumann stability analysis, so I could not find an answer. Thank you in advance.

Edit: $c_i \geq 0$ is impossible since $$ c_i \geq 0 \implies \sigma^2i^2+ri \leq 0 \implies i \leq -\frac{r}{\sigma^2} $$ Hence, $|a_i|+|b_i|+|c_i| > 1$ for any $a_i$. Please ignore the second question.

Answer 1

For the original PDE, the positivity can be deduced from the maximum principle for a parabolic operator. There is also a discrete version of the maximum principle for the finite difference parabolic operator as for example stated in Hung-Ju Kuo and N. S. Trudinger, On the discrete maximum principle for parabolic difference operators which can be applied to prove the positivity of the implicit finite difference scheme of the PDE.

It is easier to show the discrete maximal principle for the canonical heat equation which is obtained by transforming the original PDE through eliminating both $r$ and the $S$ in the coefficients.

Answer 2

For positivty, I think it will also depend on the boundary condition. The idea is briefly stated as following.

Let say we have a boundary condition $f(S_T,T) = min\{S_T - K,0\}$, $f(b,0) = 0$ and $f(a,t) = 0$. By Feynman-Kac Theorem, the solution of the original PDE can be written as:

$f = e^{-rT}\mathbb{E}[min\{S_T-k,0\}]$ , where $dS_t = rS_tdt + \sigma S_tdW_t$

We know $min\{S_T-k,0\} \leq 0, \forall S_T \in [0,\inf)$. Therefore, we have $f \leq 0$.

If your numerical method converges, your solution, in this case, will be less than or equal to 0 too.

For stability, please see PDE methods for pricing barrier options.