Short Maturity Implied Volatility in Heston Model

I am trying to compute Implied Volatilities from Heston OTM Call and OTM Put prices for short maturities ($$\tau=$$ 1 week). I am using the Carr & Madan representation, the little Heston trap formulation of the characteristic function, the Lord & Kahl algorithm for the optimal damping factor and the Kahl & Jäckel transformation in order to avoid any truncation error when integrating over the domain $$[0,\infty)$$.

The software of my choice is MATLAB_R2020a and for the numerical integration I am using the built-in function integral. I obtain positive prices, however, for deep OTM options the algorithm for implying the volatilities is failing to produce non-zero implied volatilities. I think my implementation is correct, as I am lining up with Lewis reference prices up to 12 digits, and I don't think that implementing my own Gauss-Lobatto scheme can improve the price accuracy. I ran out out ideas on how to solve this, any help is much appreciated. Thanks!

My reference Heston parameters are $$\kappa = 1.0$$, $$\theta = 0.1$$, $$\sigma = 1.0$$, $$\rho = -0.5$$ and $$v_0 = 0.1$$, $$r=0\%$$ and $$F=1$$ (same as in Lord & Kahl). I was naively taking for granted that MATLAB's implementation of Jäckel's method for implying volatilities is precise for each strike and maturity. This is not the case! I also overlooked the fact that, considering the relations $$\begin{array}{l}P_{B S}(\sigma=0, T, K, S)=\left(K e^{-r(T-t)}-S\right)^{+} \\ P_{B S}(\sigma=\infty, T, K, S)=K e^{-r(T-t)}\\ C_{B S}(\sigma=0, T, K, S)=\left(S-K e^{-r(T-t)}\right)^{+} \\ C_{B S}(\sigma=\infty, T, K, S)=S\end{array}$$ as well as that the European Put and Call option price is a monotone increasing and continuous function of volatility, obtaining positive prices for OTM options means that there exist non-zero implied volatilities. Indeed, using my implementation of the bisection method for computing implied volatilities, I obtain a nice IV smile. • Someone answering their own question gets a +1 from me! Well done Jun 21 '20 at 19:57

As a coincidence I have also started to learn about and implement stochastic vol models. So I first want to say that I am no expert in the field. However, I just struggled through implementations of calibration and simulation of the Heston model. For the calibration I used the article "Full and fast calibration of the Heston stochastic volatility model" by Y.cui, et al. They have a formulation of the characteristic function they claim lack the numerical instabilities other earlier articles suffers from (due to some ways of rewriting the expressions involving the complex logarithm to avoid discontinuities due to branch cuts).

However, I also ran into the same problems as you. There are formulas that express the option prices in terms of integrals involving the characteristic function. I guess those are the ones you also use for calculations. While the characteristic function itself might be fine, the integrals needs to be computed numerically. The authors of the article "full and fast.." recommend standard Gauss quadrature on $$[-1,1]$$ rescaled to the half axis.

I noticed that in some cases, even though I was using quite many nodes, the numerical integrals results here were very bad. So I investigated the issue a bit for some limit extreme cases. It turned out the integrals we needed to compute were essentially behaving like $$\int_0^{\infty} \frac{\sin{x}}{x} \, dx$$ This is a tough integral to compute numerically and it is not suitable to use ordinary Gauss quadrature. The convergence is very slow, similar to the series $$\sum_{k=1}^{\infty} \frac{(-1)^k}{k}$$ , I would guess.

The problem here is both the slow convergence and that it is a oscillatory integral which is bad for standard Gauss quadrature if the oscilliations are too frequent compared to the number of quadrature nodes.

But on the other hand you say that your prices match the reference prices to 12 decimals, so one has to assume the numerical integrals have converged nicely in your case.

What procedure are you using to calculate the implied vol from the prices? I also often get problems in extreme cases (which are not that extreme). The problem for me trying compute it numerically is that while evaluating the Black formula, the $$d_1$$ and $$d_2$$ becomes very large or very negative so that the normal cdf becomes equal to $$1$$ or $$0$$ and cannot separate different volatilites anymore due to lack of precision. Then the numerical root finder might try to search for negative volatilities.

I also got some strange cases where the integrals seemed to converge, but option prices did not correspond to any positive Black vol. Those cases were of the type where (assuming zero interest rate) the underlying had value $$S_0$$ say, the strike was $$K$$, but the calculated call option price became slighly below $$S_0-K$$. This seems impossible to match with a positive volatility since higher vol just increases the option price and zero vol would give the option price $$S_0-K$$.

As a first step I would therefore calculate the Black-prices using zero Black volatility (which will just be based on the forward value then since there is no uncertainity). If, in the call option case, the obtained values corresponding to vol=0 are larger than your calculated values from the Heston model, then the problem of calculating the implied vol is unsolvable and there is some issue in the option pricing in the Heston model. However, if all the black call prices using vol=0 are smaller than your Heston option prices, then the issue lies in your calculation of implied Black vol.

• Thanks Jesper. The formulation of the Heston CF by Albrecher et. al is also not affected by discontinuities, as also pointed out in the paper you mentioned. To handle the oscillatory integrants I suggest you to take a look at Lord & Kahl. I am using MATLAB's implementation of Jäckel's method for computing implied volatilities which, to the best of my knowledge, is the standard market practice (...) Jun 21 '20 at 7:16
• (...) However, if I fix $\sigma_{BS}=0.7$, compute the BS for $K=0.1$ and $\tau=1/52$ , I realised that the MATLAB implementation of Jäckel's method is failing to back out the $0.7$ IV. With my bisection implementation I am able to imply the volatilities from the Heston OTM Call and OTM put prices and I obtain a nice IV smile. Thanks for pointing me to the problem, I was taking for granted that MATLAB's implementation of Jäckel's method was precise for each strike and maturity. Jun 21 '20 at 7:27
• Very nice that you solved your problem. Thanks for the reference regarding oscillatory integrals. I will take a look at that paper. Jun 21 '20 at 16:17