Put-Call parity and the log-transformed Black-Scholes PDE
We will assume zero interest rates and no dividends on the asset \(S\) for clarity. The results can be easily generalized to the case with non-zero interest rates and dividends. Under those assumptions, the Black-Scholes PDE is: $$ \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} = 0.$$
An implicit Euler discretisation on a uniform grid in \(S\) of width \(h\) with linear boundary conditions (zero Gamma) leads to:
$$ V^{k+1}_i - V^{k}_i = \frac{1}{2}\sigma^2 \Delta t S_i^2 \frac{V^{k}_{i+1}-2V^k_{i}+V^{k}_{i-1}}{h^2}.$$ for \(i=1,…,m-1\) with boundaries $$ V^{k+1}_i - V^{k}_i = 0.$$ for \(i=0,m\).
This forms a linear system \( M \cdot V^k = V^{k+1} \) with \(M\) is a tridiagonal matrix where each of its rows sums to 1 exactly. Furthermore, the payoff corresponding to the forward price \(V_i = S_i\) is exactly preserved as well by such a system as the discretized second derivative will be exactly zero. The scheme can be seen as preserving the zero-th and first moments.
As a consequence, by linearity, the put-call parity relationship will hold exactly (note that in between nodes, any interpolation used should also be consistent with the put-call parity for the result to be more general).
This result stays true for a non-uniform discretisation, and with other finite difference schemes as shown in this paper.
It is common to consider the log-transformed problem in \(X = \ln(S)\) as the diffusion is constant then, and a uniform grid much more adapted to the process. $$ \frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 \frac{\partial^2 V}{\partial X^2}-\frac{1}{2}\sigma^2 \frac{\partial V}{\partial X} = 0.$$
An implicit Euler discretisation on a uniform grid in \(X\) of width \(h\) with linear boundary conditions (zero Gamma in \(S\) ) leads to:
$$ V^{k+1}_i - V^{k}_i = \frac{1}{2}\sigma^2 \Delta t \frac{V^{k}_{i+1}-2V^k_{i}+V^{k}_{i-1}}{h^2}-\frac{1}{2}\sigma^2 \Delta t \frac{V^{k}_{i+1}-V^{k}_{i-1}}{2h}.$$ for \(i=1,…,m-1\) with boundaries $$ V^{k+1}_i - V^{k}_i = 0.$$ for \(i=0,m\).
Such a scheme will not preserve the forward price anymore. This is because now, the forward price is \(V_i = e^{X_i}\). In particular, it is not linear in \(X\).
It is possible to preserve the forward by changing slightly the diffusion coefficient, very much as in the exponential fitting idea. The difference is that, here, we are not interested in handling a large drift (when compared to the diffusion) without oscillations, but merely to preserve the forward exactly. We want the adjusted volatility \(\bar{\sigma}\) to solve $$\frac{1}{2}\bar{\sigma}^2 \frac{e^{h}-2+e^{-h}}{h^2}-\frac{1}{2}\sigma^2 \frac{e^{h}-e^{-h}}{2h}=0.$$ Note that the discretised drift does not change, only the discretised diffusion term. The solution is: $$\bar{\sigma}^2 = \frac{\sigma^2 h}{2} \coth\left(\frac{h}{2} \right) .$$ This needs to be applied only for \(i=1,…,m-1\).
This is actually the same adjustment as the exponential fitting technique with a drift of zero. For a non-zero drift, the two adjustments would differ, as the exact forward adjustment will stay the same, along with an adjusted discrete drift.