Summary
For Heston model parameters that render the variance process constant, the solution should revert to plain Black-Scholes. Closed from solutions to the Heston model don't seem to do this, even though the dynamics clearly do. Can anyone see how to resolve this? I care because I'm getting incorrect results for four independent implementations spanning three completely different methods of approaching the problem, but not getting incorrect results for Monte Carlo pricing based on time evolution of the SDEs.
Details
I am implementing some option pricers based on the Heston model and using different methods:
- Monte Carlo simulation
- The original closed-form solution presented in Heston 1993
- The Green's function as presented in this very nice MSc thesis
- Fourier methods
Unfortunately what I am finding is that while the final three (semi-)analytical solutions agree with one another, they do not agree with the Monte Carlo simulation and my intuition about the problem.
Given the Heston dynamics: $$ dS_t = \mu S_t,dt + \sqrt{\nu_t} S_tdW^S_t,\\ d\nu_t = \kappa(\theta - \nu_t)dt + \sigma \sqrt{\nu_t}dW^{\nu}_t,\\ dW^S_t dW^{\nu}_t = \rho dt\\ $$ it is clear that if we simply choose our parameters for the variance process $\nu_t$ such that it remains constant ($d\nu_t = 0, \nu_t = \nu_0$), we should return to the plain Black-Scholes dynamics $$ dS_t = \mu S_tdt + \sqrt{\nu_0} S_tdW^S_t . $$ From this we can derive the analytical solution in the usual way and up with the following $$ C(S,t) = \mathcal{N}(d_1)~S-\mathcal{N}(d_2)~K e^{-r \tau} \\ d_1 = \frac{\ln(\frac{S}{K})+(r+\frac{\nu_0}{2})(\tau)}{\sqrt{\nu_0}\sqrt{\tau}} \\ d_2 = \frac{\ln(\frac{S}{K})+(r-\frac{\nu_0}{2})(\tau)}{\sqrt{\nu_0}\sqrt{\tau}} = d_{1}-\sqrt{\nu_0}\sqrt{\tau}. $$ This is all well and good and is the reason that I feel comfortable in claiming that my three different implementations of the semi-analytical Heston solution must be wrong, despite their accord with one another. Further, the Monte Carlo simulation of the model does give the same price as the Black-Scholes equation in this case.
Analogous to the Black-Scholes model, the solution of Heston model as given in the original paper is of the form $$ C(S,v,t) = S P_1 - K e^{-r(\tau)} P_2, $$ where $$ P_j = \frac{1}{2} + \frac{1}{\pi} \int^\infty_0 Re [ \frac{e^{-i \phi ln(K)} f_j(x,\nu,T,\phi)}{i \phi} ] d \phi, \\ f_j = e^{C + D \nu + i \phi x}, \\ C = r \phi i \tau + \frac{\kappa \theta}{\sigma^2} {(b_j - \rho \sigma \phi i + d)\tau - 2 \log [\frac{1-ge^{d \tau}}{1-g} ] }, \\ D = \frac{b_j - \rho \sigma \phi i + d}{\sigma^2} [ \frac{1-e^{d \tau}}{1-ge^{d \tau}} ], \\ g = \frac{b_j - \rho \sigma \phi i + d}{b_j - \rho \sigma \phi i - d}, \\ d = \sqrt{(\rho \sigma \phi i - b_j)^2 - \sigma^2(2 u_j \phi i - \phi^2)}, \\ b_1 = \kappa - \rho \sigma \; \; \; b_2 = \kappa, \\ u_1 = \frac{1}{2} \; \; \; u_2 = - u_1. $$ so it seems rather apparent that for cases where the Heston model is equivalent to Black-Scholes, we must have $P_1 = \mathcal{N}(d1)$ and $P_2 = \mathcal{N}(d2)$.
To demonstrate, let's take $ \theta = \nu_0, \sigma = 0 $, implying $\forall \nu_0, d\nu_t = 0$.
When just considering the Heston dynamics, these parameters work nicely and everything is well behaved. The semi-analytical solutions, however, are plagued with division-by-zero (ex. $g$) and can't be handled directly. So instead we simply want to take them in the limit as $\sigma \to 0$. This is where it gets tricky now, as we're left with some very complicated formulas and no apparent link to $\mathcal{N}(d_j)$. I tried Mathematica's rather powerful symbolic calculus capabilities, but it was unable to show equality.
Testing some other limiting cases, I've made the following observations:
- For $\nu_0 \to 0$, all methods converge on the correct price (Black Scholes).
- For low strikes, $K \to 0$, all methods converge on the correct price ($S_0$).
To me, this suggests that the problem is in $P_2$, but that's all I have so far.
What direction should I take next?