Most accurate Fourier transform method for extreme OTM options

Question

I need to calculate vanilla options prices for extreme moneyness range of e.g. (0%,1000%) under the Heston model for various parameter values that satisfy Feller.

Which Fourier method (or other method) would be the most accurate for this? For instance the COS method is nice but I think it has some issues for extreme moneyness.

I am not too concerned with speed of computation.

Thanks.

Answer 1

I'll give it a start and stick with Fourier methods. The approaches from Carr and Madan (1999) and Fang and Oosterlee (2009) are indeed known to be inaccurate for highly OTM options. I'd suggest to try out one the following three alternatives. The first one seems to be the most relevant one.

Saddlepoint Method

I begin to cite Hirsa (2013):

The saddlepoint method offers considerably better accuracy in pricing out-of-the-money options than either the FFT or the COS method. However, the algorithm’s accuracy for at the money and in the money options is somewhat lacking compared to these two other methods and like the FFT method its solution must be rederived for each different payoff.

Like almost everything in Fourier based option pricing, this idea stems from Carr and Madan (2009). It is based on Madan's et al. (2008) insight that option prices are tail probabilities.

Under the stock measure, $$\frac{C(S_0,K,T)}{S_0}=\mathbb{E}_0^\mathbb{S}\left[\frac{\max\{S_T-K,0\}}{S_T}\right]=\mathbb{E}_0^\mathbb{S}\left[\left(1-\frac{K}{S_T}\right)^+\right].$$ If $y=\ln\left(\frac{S_T}{K}\right)$ and $\frac{K}{S_T}=e^{-y}$ and $f$ is the density of $\ln\left(\frac{S_T}{K}\right)$ under $\mathbb{S}$, then, using partial integration, \begin{align*} \frac{C(S_0,K,T)}{S_0}&=\int_0^\infty (1-e^{-y})f(y)\mathrm{d}y \\ &= 1-\int_0^\infty e^{-y}F(y)\mathrm{d}y \\ &= \int_0^\infty (1-F(y))e^{-y}\mathrm{d}y. \end{align*}

Given $y$, $1-F(y)$ is just the probability (under $\mathbb{S}$) of the event $\left\{\frac{S_T}{K}>y\right\}$. Furthermore, $e^{-y}$ is just the pdf of an exponential random variable with mean 1. Thus, \begin{align*} \frac{C(S_0,K,T)}{S_0}&= P\left[\left\{\ln\left(\frac{S_T}{K}\right)>Y\right\}\right]=P\left[\left\{X-Y>\ln(K)\right\}\right], \end{align*} where $X=\ln(S_T)$ under $\mathbb{S}$ and $Y$ an independent exponential random variable.

That's it. The Saddlepoint method is based on the following idea. Let $M$ be the moment generating function of $X-Y$. Then the cumulant generating function is defined via $K(t)=\ln(M(t))$. The Lugannani Rice formula uses CGFs to approximate tail probabilities. That's precisely what we need. Carr and Madan (2009) write down all the equations. Hirsa (2013) also provides a step by step implementation.

I jump straight to the results. This is Table 2.10 from Hirsa. .

We got Monte Carlo, the COS method, fractional FFT, plain FFT and the saddlepoint method. Parameters are $S_0=100$, $r=3\%$, $\kappa=2$, $\xi=0.5$, $\bar{v}=0.04$, $v_0=0.04$, $\rho=-0.7$ and $T=\frac{1}{2}$. I can state all the equations here if you don't have Hirsa's book. I guess this method is as reasonably accurate as it gets with Fourier methods.

---

Carr Madan Time Value

Carr and Madan (1999) acknowledge that their FFT method is not flawless and propose an alternative approach that works with time values only. That could improve the pricing of OTM options. The OTM prices are simply \begin{align} C^\text{OTM}(S_0,K,T) &= \mathrm{1}_{\{K>S_0\}} \cdot C(S_0,K,T) \\ P^\text{OTM}(S_0,K,T) &= \mathrm{1}_{\{K<S_0\}}\cdot P(S_0,K,T). \end{align} Assuming that Fourier transforms exist, you can write \begin{align*} C^\text{OTM}(S_0,K,T) &= \frac{1}{\pi}\int_0^\infty \Re\left(e^{-iku} \hat{C}^\text{OTM}(u)\right)\mathrm{d}u, \\ P^\text{OTM}(S_0,K,T) &= \frac{1}{\pi}\int_0^\infty \Re\left( e^{-iku} \hat{P}^\text{OTM}(u)\right)\mathrm{d}u, \end{align*}

You can combine both equations to $O(S_0,K,T)=C^\text{OTM}(S_0,K,T)+P^\text{OTM}(S_0,K,T)$. For $K>S_0$, $O$ equals the time value of the call option and for $K<S_0$ it is the time value of a put option. The function $O(K)$ is not defined for $K=S_0$. By linearity, $\hat{O}(u)=\hat{C}^\text{OTM}(u) +\hat{P}^\text{OTM}(u)$. Carr and Madan compute that \begin{align} \hat{O}(u) &= e^{-rT}\left( \frac{S_0^{i(u-i)}}{i(u-i)}-\frac{S_0^{i(u-i)}e^{(r-q)T}}{iu}-\frac{\varphi_T(u-i)}{u(u-i)} \right), \label{Eqn: Fourier Transform of Time Value} \end{align} This is Equation 14 in Carr and Madan (1999), where I add an $S_0\neq1$ and dividends $q\neq0$. The calculation is straight forward, you only need to change the order of integration.

However, the function $O$ is obviously peaks at the money, see below (I use a GBM with $S_0=2$).

This effect gets worse as $T\to0$ because the time value rapidly converges to zero. The function $O$ looks more and more like a Dirac delta function -- and the Fourier transform spreads out further and further (uncertainty principle).

To address this issue Carr and Madan use $\sinh(\alpha(k-\ln(S_0))$ for damping, rather than the standard $e^{\alpha k}$. The Fourier transform of the damped OTM option price is \begin{align*} \psi_T(u) &= \int_\mathbb{R} e^{iku} \sinh\left(\alpha \big(k-\ln(S_0)\big)\right)O(k) \mathrm{d}k\\ &= \frac{\hat{O}(u-i\alpha)e^{-\alpha \ln(S_0)}-\hat{O}(u+i\alpha)e^{\alpha \ln(S_0)}}{2}, \end{align*} which results in \begin{align*} O(k) &= \frac{1}{\pi\sinh\left(\alpha \big(k-\ln(S_0)\big)\right)}\int_0^\infty \Re\left(e^{-iku} \psi_T(u)\right) \mathrm{d}u, \end{align*}

---

Control Variates

Somewhat related is the idea to subtract the current intrinsic value from the option to obtain the OTM option time value. Then, \begin{align*} C(S_0,K,T) = \left(S_0e^{-qT}-Ke^{-rT}\right)^+- \frac{S_0e^{-qT}}{\pi}\int_0^\infty \Re\left( e^{iu\tilde{k}} \frac{\varphi_{X_T}(u-i)-1}{u(u-i)} \right)\mathrm{d}u, \end{align*} where $\tilde{k}=\ln\left(\frac{S_0e^{-qT}}{Ke^{-rT}}\right)$.

Cont and Tankov (2004, Chapter 11.1.3) point out that's better to use a smooth function such as the BS call price. Mark Joshi, too, liked using the BS model as a control variate, see also Joshi and Yang (2011). The call option price looks somewhat like this \begin{align*} C(S_0,K,T) = C_\text{BS}(S_0,K,T)- \frac{S_0e^{-qT}}{\pi}\int_0^\infty \Re\left( e^{i\tilde{k}u} \frac{\varphi_{X_T}(u-i)-\varphi_{X_T}^\text{BS}(u-i)}{u(u-i)} \right)\mathrm{d}u, \end{align*} where $\tilde{k}=\ln\left(\frac{S_0e^{-qT}}{Ke^{-rT}}\right)$. Furthermore, $C_\text{BS}$ denotes the standard Black Scholes (1973) option price and $\varphi_{X_T}^\text{BS}(u)=e^{-\frac{1}{2}\sigma^2 T u(u+i)}$. There's some literature on how to find the optimal volatility for the BS call, see below for the Heston model:

Answer 2

The appropriate technique is the optimal $\alpha$ damping of Lord and Kahl (2007).

Since then, this technique has been refined by Andersen and Lake (2019).

The saddlepoint method is not precise. It is an approximation to the price. It only gives a rough idea of what the price should be. This is also detailed in the Lord and Kahl (2007) paper.