Approximation Formula for Option Prices under Rough Heston Model and Short-Time Implied Volatility Behavior

Abstract

Rough Heston model possesses some stylized facts that can be used to describe the stock market, i.e., markets are highly endogenous, no statistical arbitrage mechanism, liquidity asymmetry for buy and sell order, and the presence of metaorders. This paper presents an efficient alternative to compute option prices under the rough Heston model. Through the decomposition formula of the option price under the rough Heston model, we manage to obtain an approximation formula for option prices that is simpler to compute and requires less computational effort than the Fourier inversion method. In addition, we establish finite error bounds of approximation formula of option prices under the rough Heston model for $0.1\le H<0.5$ under a simple assumption. Then, the second part of the work focuses on the short-time implied volatility behavior where we use a second-order approximation on the implied volatility to match the terms of Taylor expansion of call option prices. One of the key results that we manage to obtain is that the second-order approximation for implied volatility (derived by matching coefficients of the Taylor expansion) possesses explosive behavior for the short-time term structure of at-the-money implied volatility skew, which is also present in the short-time option prices under rough Heston dynamics. Numerical experiments were conducted to verify the effectiveness of the approximation formula of option prices and the formulas for the short-time term structure of at-the-money implied volatility skew.

1. Introduction

Stochastic volatility models and jump diffusion models such as classical Heston model [1] and Merton jump diffusion model [2] have played an important role in the option pricing theory. The models have aimed at replicating the stochastic volatility effect (along with the mean-reversion) and the jump-effect as displayed in the financial market. For decades, the models have been proven to be useful in pricing options, but they are deemed as inadequate for modeling the short-time behavior of the implied volatility of the options using the real data.

It was discovered by Alos et al. [3] that the introduction of fractional components with Hurst parameter $0<H<0.5$ is capable of generating the term structure of at-the-money skew of order $T^{\delta }$ for every $\delta >-0.5$ and small time to maturity T. In addition, Fukasawa [4] managed to show that (through an example of martingale expansion’s application) the stochastic volatility model with volatility term driven by fractional Brownian motion with Hurst parameter H can generate at-the-money volatility skew of order $T^{H-0.5}$ for the small time to maturity T. This motivated the authors of [5] to further explore through real empirical data to verify that the use of fractional Brownian motion in a stochastic volatility model is useful and adequate in replicating the roughness behavior of volatility in the financial market. Fractional Brownian motion with Hurst parameter $0<H<0.5$ in volatility model is subsequently named as rough volatility. On a side note, Hurst parameter H is used to control the roughness of the fractional Brownian motion, i.e., when $0<H<0.5$, the rough volatility displays anti-persistency effect, and, for $0.5<H<1$, the volatility would display long-memory effect. Prior to the use of rough volatility, fractional Brownian motion with Hurst parameter $0.5<H<1$ was proposed by Comte and Renault [6] to model the long-memory effect of volatility. A linear sum of classical Brownian motion and fractional Brownian motion called mixed fractional Brownian motion $0.5<H<1$ was also proposed and developed by Cheridito et al. [7] to price model stock movement and price options.

Rough Bergomi model by Bayer, C., & Friz, P., Gatheral, J. [8] was among the first models proposed to incorporate the rough volatility component. The initial model proposed is slow in simulation such that the authors ruled out the possibility of optimization in practice. Nevertheless, the authors proposed an efficient expansion method for the implied volatility to the second order in the volatility of volatility component. Furthermore, further advancements in reducing the computational costs were made by the authors of [9,10], i.e., they introduced a faster scheme that utilizes hybrid methods for the Brownian semistationary process in the simulation and variance reduction techniques for a significant runtime reduction in the simulation process.

Rough Heston model is our main discussion of this paper. It possesses the rough quality of the classical Heston model and was introduced by El Euch and Rosenbaum [11] as a stochastic Volterra equation. The model is used to describe the high degree of endogeneity (orders that have no real economic motivation) of the liquid market and the no statistical arbitrage mechanism which should present on a highly liquid electronic markets. In addition, the liquid asymmetry (average impact of a sell order is greater than the impact of a buy order) and presence of metaorders (large orders whose execution is split in time by trading algorithm) are also features of rough Heston model. Accordingly, the authors of [11] obtained the characteristic function of the rough Heston model to compute the option price using the Fourier inversion method. Coincidentally, the computation for option price under the rough Heston model is also computationally expensive. One of the brilliant achievements of El Euch and Rosenbaum [11] is that, under the numerical experiments, the authors were able to numerically verify that pricing option price under rough Heston model consists of explosive behavior for term-structure at-the-money implied volatility skew (the same goes to the rough Bergomi model) as $T\to 0$. Explosive behavior for the short-time term structure of at-the-money implied volatility skew is one of the crucial components in pricing options for short maturity time (see [4,8,12]). In addition, Gatheral and Radoicic [13] used rational approximations on fractional Riccati equation which in return substantially reduces the computational cost (see [14] for a review and [15] for an improvement on the method). Recent advancements by the authors of [16,17] made effort in reducing the computational cost of the rough Heston model.

We study an approximation formula for option prices under the rough Heston model. Section 1 serves as prerequisites to our study. In particular, the approximation method is obtained using the decomposition formula proposed in [18,19,20]. Accordingly, the error bounds of the approximation formula are studied in Section 2 under a simple assumption. Details of the proof in Section 2 are listed in Appendix A. Subsequently, in Section 3, we propose a second-order approximation to the implied volatility similar to the work of Alos, Elisa and De Santiago, Rafael and Vives, Josep [19]. We then match the coefficients and thus obtain some of the short-time implied volatility behaviors for our approximation formula for implied volatility. Ultimately, one important result that we obtain is that the approximation formula for implied volatility does display explosive behavior for the term structure at-the-money volatility skew. Through the numerical experiments, as shown in Section 4, we manage to verify that the option prices under rough Heston model computed using the approximation formula and the Fourier inversion method with fractional Adams scheme [11,21] are extremely consistent with one another for short maturity time. In addition, we find that, by comparing the implied volatility computed by Fourier inversion method with fractional Adams scheme [11,21] under rough Heston model, approximation formula for option prices (transformed back to implied volatility) and short-time term structure of the at-the-money skew equation have roughly the same order of explosive behavior for the term structure of at-the-money implied volatility skew. Lastly, we give a conclusion to our work in Section 5.

Prerequisites

Consider the price model process with rough Heston dynamics:

$$\frac{d{S}_u}{S_u}={\sigma }_u\left[\rho d{W}_u^{*}+\sqrt{1-{\rho }^2}d{W}_u^{\perp }\right]$$

(1)

$${\sigma }_u^2={\sigma }_t^2+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}\lambda \left[\theta -{\sigma }_s^2\right]ds+\frac{\nu }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\sigma }_{s}d{W}_s^{*}$$

(2)

where S is the stock price, $\rho$ is the correlation between stock return and the volatility movement, $0.5<\alpha =H+0.5<1$, $\theta >0$ is the mean reversion level, $\lambda >0$ is the rate at which the process ${\sigma }_s^2$ reverts to the mean reversion level $\theta$, $\nu >0$ is the magnitude of random movement, $W^{*}$ and $W^{\perp }$ are independent Brownian motion defined in a filtered probability space ($\mathsf{\Omega },\mathcal{F},\mathbb{Q}$) with $\mathcal{F}:={\mathcal{F}}^{W^{*}}\vee {\mathcal{F}}^{W^{\perp }}$, and $\sigma$ is the volatility of the return that is positive and square-integrable function adapted to the filtration.

Let $X_u=log\left(S_u\right)$ and we describe Equation (2) in the forward variance form:

$$d{X}_u={\sigma }_u\left[\rho d{W}_u^{*}+\sqrt{1-{\rho }^2}d{W}_u^{\perp }\right]-\frac{1}{2}{\sigma }_u^{2}du$$

(3)

$${\sigma }_u^2={\xi }_t\left(u\right)+\frac{\nu }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\sigma }_{s}d{W}_s^{*}$$

(4)

where ${\xi }_t\left(u\right):=\mathbb{E}\left[{\sigma }_u^2\right|{\mathcal{F}}_t]$ is the forward variance curve. Following the work of Alos et al. [20], we denote the future expected variance $w_t\left(T\right)$ as

$$w_t\left(T\right):=\mathbb{E}\left[{\int }_t^T{\sigma }_s^{2}ds|{\mathcal{F}}_t\right]$$

(5)

and the martingale as

$$M_t:=\mathbb{E}\left[{\int }_0^T{\sigma }_s^{2}ds|{\mathcal{F}}_t\right]$$

(6)

The connection between $w_t\left(T\right)$ and $M_t$ can be easily established as

$$w_t\left(T\right)=M_t-{\int }_0^T{\sigma }_s^{2}ds$$

(7)

then it follows that

$$d{w}_t\left(T\right)=d{M}_t-{\sigma }_t^{2}dt$$

(8)

The Delta–Gamma–Vega relationship can be established as

$$\frac{1}{2}({\partial }_{xx}-{\partial }_x)B(X_t,w_t)-{\partial }_{w}B(X_t,w_t)=0$$

(9)

where $B_t:=B(X_t,w_t\left(T\right))$ is some function that solves the Black–Scholes equation with X being the log-spot price with the dynamics specified in Equation (3) and w depends on $\sigma$. Otherwise, $B(X_t,w_t)$ can be formulated as

$$B(X_t,w_t\left(T\right))=e^{X_t}·N\left(d_{+}\left(w_t\left(T\right)\right)\right)-K·N\left(d_{-}\left(w_t\left(T\right)\right)\right)$$

(10)

where $N(·)$ is the cumulative distribution function of standard normal distribution and

$$d_{\pm }\left(w_t\left(T\right)\right)=\frac{X_t-lnK}{\sqrt{w_t\left(T\right)}}\pm \frac{\sqrt{w_t\left(T\right)}}{2}$$

(11)

with the assumption of risk-free interest rate $r=0$ throughout the paper.

Theorem 1

(Decomposition formula). Assume $B_t:=B(X,w)$ as some function that solves the Black–Scholes equation where x is the log-spot price with the dynamics specified in Equation (3) and w depends on the ${\sigma }^2$ from Equation (4). Then, for all $t\in [0,T],$

$$\begin{array}{cc}\hfill \mathbb{E}\left[B_T\right|{\mathcal{F}}_t]& =B_t+\mathbb{E}\left[{\int }_t^T{\partial }_{xw}{B}_{s}d\langle X,M\rangle |{\mathcal{F}}_t\right]\hfill \\ \hfill & \phantom{=}+\frac{1}{2}\mathbb{E}\left[{\int }_t^T{\partial }_{ww}{B}_{s}d{\langle M,M\rangle }_s|{\mathcal{F}}_t\right]\hfill \end{array}$$

(12)

Proof. 

Applying Itô’s formula to $B_t$, we have

$$\begin{array}{cc}\hfill d{B}_t=& {\partial }_x{B}_{t}d{X}_t+{\partial }_w{B}_{t}d{w}_t\left(T\right)+\frac{1}{2}{\partial }_{xx}{B}_{t}d{\langle X,X\rangle }_t\hfill \\ \hfill & +{\partial }_{xw}{B}_{t}d{\langle X,w\left(T\right)\rangle }_t+\frac{1}{2}{\partial }_{ww}{B}_{t}d{\langle w\left(T\right),w\left(T\right)\rangle }_t\hfill \end{array}$$

(13)

We use Equations (3), (8) and (9) to simplify

$$\begin{array}{cc}\hfill d{B}_t& ={\partial }_x{B}_t\left({\sigma }_t\left[\rho d{W}_t^{*}+\sqrt{1-{\rho }^2}d{W}_t^{\perp }\right]-\frac{1}{2}{\sigma }_t^{2}dt\right)+{\partial }_w{B}_t\left(d{M}_t-{\sigma }_t^{2}dt\right)\hfill \\ & \phantom{=}+\frac{1}{2}{\partial }_{xx}{B}_t{\sigma }_t^{2}dt+{\partial }_{xw}{B}_{t}d{\langle X,M\rangle }_t+\frac{1}{2}{\partial }_{ww}{B}_{t}d{\langle M,M\rangle }_t\hfill \\ & ={\partial }_x{B}_t{\sigma }_t\left[\rho d{W}_t^{*}+\sqrt{1-{\rho }^2}d{W}_t^{\perp }\right]+{\sigma }_t^2\left(\frac{1}{2}({\partial }_{xx}-{\partial }_x)B_t-{\partial }_w{B}_t\right)dt\hfill \\ & \phantom{=}+{\partial }_w{B}_{t}d{M}_t+{\partial }_{xw}{B}_{t}d{\langle X,M\rangle }_t+\frac{1}{2}{\partial }_{ww}{B}_{t}d{\langle M,M\rangle }_t\hfill \\ & ={\partial }_x{B}_t{\sigma }_t\left[\rho d{W}_t^{*}+\sqrt{1-{\rho }^2}d{W}_t^{\perp }\right]+{\partial }_w{B}_{t}d{M}_t\hfill \\ & \phantom{=}+{\partial }_{xw}{B}_{t}d{\langle X,M\rangle }_t+\frac{1}{2}{\partial }_{ww}{B}_{t}d{\langle M,M\rangle }_t\hfill \end{array}$$

(14)

Integrating up to time T and taking conditional expectation of Equation (14) on ${\mathcal{F}}_t$, the result follows.□

Theorem 1 is the decomposition formula from [20] with slight changes of the notation as corresponds to the original work in [18]. Equation (12) is exact rather than an approximation, i.e., it describes the interaction between the movement of log-stock movement and the future expected variance from time t to T. It is known that Equation (12) under conditional expectation can be difficult to compute, therefore approximation of the equation is needed. In this study, we provide an approximation formula for Equation (12) using the rough Heston model dynamics under a simple assumption. It is found that the errors are bounded and decreasing as the time to maturity gets smaller. Let ${\mathbb{E}}_t(·):=\mathbb{E}[·|{\mathcal{F}}_t]$ as a simplification.

We consider the following assumption throughout this study:

Assumption 1

(A1). ${\mathbb{E}}_t\left({\sigma }_u^2\right){\int }_t^u{(u-s)}^{\alpha -1}ds\ge {\int }_t^u{(u-s)}^{\alpha -1}{\mathbb{E}}_t\left({\sigma }_s^2\right)ds$

The relevance of Assumption 1 for the practical application is discussed at the end of Section 4. We prove some bounds under Assumption 1.

Lemma 1.

Assume the model in (2) and (3) and Assumption 1; then, for all $T-s\in [0,1]$,

$${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du\ge {\sigma }_s^2\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha +1)+\lambda }(T-s)$$

(15)

$${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du\ge \frac{\lambda \theta {(T-s)}^{\alpha +1}}{\mathsf{\Gamma }(\alpha +2)+\lambda (\alpha +1)}$$

(16)

and for $T-s\in (1,\infty )$

$${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du\ge {\sigma }_s^2\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha +1)+\lambda }$$

(17)

$${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du\ge \frac{\lambda \theta }{\mathsf{\Gamma }(\alpha +2)+\lambda (\alpha +1)}$$

(18)

Proof. 

Since ${\mathbb{E}}_t\left[{\sigma }_u^2\right]$ is ${\mathcal{F}}_t$-martingale, it is not hard to see that from Equations (2) and (3)

$$\begin{array}{cc}\hfill {\mathbb{E}}_t\left({\sigma }_u^2\right)& ={\sigma }_t^2+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}\lambda [\theta -{\mathbb{E}}_t\left({\sigma }_s^2\right)]ds\hfill \\ & ={\sigma }_t^2+\frac{\lambda \theta }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}ds-\frac{\lambda }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\mathbb{E}}_t\left({\sigma }_s^2\right)ds\hfill \\ & \ge {\sigma }_t^2-\frac{\lambda }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\mathbb{E}}_t\left({\sigma }_s^2\right)ds\hfill \end{array}$$

(19)

Based on Assumption 1,

$$\begin{array}{cc}\hfill {\mathbb{E}}_t\left({\sigma }_u^2\right)& \ge {\sigma }_t^2-\frac{\lambda }{\mathsf{\Gamma }\left(\alpha \right)}{\mathbb{E}}_t\left({\sigma }_u^2\right){\int }_t^u{(u-s)}^{\alpha -1}ds\hfill \\ \hfill & ={\sigma }_t^2-\frac{\lambda {(u-t)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)}{\mathbb{E}}_t\left({\sigma }_u^2\right)\hfill \end{array}$$

(20)

Remembering that $0.5<\alpha <1$, rearranging we obtain

$$\begin{array}{cc}\hfill {\mathbb{E}}_t\left({\sigma }_u^2\right)& \ge {\sigma }_t^2\left(\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha +1)+\lambda {(u-t)}^{\alpha }}\right)\hfill \\ \hfill & \ge {\sigma }_t^2\mathsf{\Gamma }(\alpha +1)\left[\left(\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda }\right){𝟙}_{u-t\le 1}+\left(\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda (u-t)}\right){𝟙}_{u-t>1}\right]\hfill \end{array}$$

(21)

For $T-s\le 1$, we compute ${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du$ as

$$\begin{array}{cc}\hfill {\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du& \ge {\int }_s^T{\sigma }_s^2\mathsf{\Gamma }(\alpha +1)\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda }du\hfill \\ \hfill & ={\sigma }_s^2\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha +1)+\lambda }(T-s)\hfill \end{array}$$

(22)

For $T-s>1$,

$$\begin{array}{cc}\hfill {\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du& \ge {\sigma }_s^2\mathsf{\Gamma }(\alpha +1)\left[{\int }_s^{s+1}\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda }du+{\int }_{s+1}^T\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda (u-s)}du\right]\hfill \\ & ={\sigma }_s^2\mathsf{\Gamma }(\alpha +1)\left[\frac{1}{\mathsf{\Gamma }(\alpha +1)+\lambda }+\frac{1}{\lambda }ln\left(\frac{\mathsf{\Gamma }(\alpha +1)+\lambda (T-s)}{\mathsf{\Gamma }(\alpha +1)+\lambda }\right)\right]\hfill \\ & \ge {\sigma }_s^2\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha +1)+\lambda }\hfill \end{array}$$

(23)

Furthermore, we can obtain two more bounds by considering the following

$$\begin{array}{cc}\hfill {\mathbb{E}}_t\left({\sigma }_u^2\right)& \ge \frac{\lambda \theta }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}ds-\frac{\lambda }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\mathbb{E}}_t\left({\sigma }_s^2\right)ds\hfill \\ \hfill & =\frac{\lambda \theta }{\mathsf{\Gamma }(\alpha +1)}{(u-t)}^{\alpha }-\frac{\lambda }{\mathsf{\Gamma }\left(\alpha \right)}{\int }_t^u{(u-s)}^{\alpha -1}{\mathbb{E}}_t\left({\sigma }_s^2\right)ds\hfill \end{array}$$

(24)

By Assumption 1,

$$\begin{array}{cc}\hfill {\mathbb{E}}_t\left({\sigma }_u^2\right)& \ge \frac{\lambda \theta }{\mathsf{\Gamma }(\alpha +1)}{(u-t)}^{\alpha }-\frac{\lambda }{\mathsf{\Gamma }(\alpha +1)}{(u-t)}^{\alpha }{\mathbb{E}}_t\left({\sigma }_u^2\right)\hfill \\ & =\frac{\lambda \theta {(u-t)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)+\lambda {(u-t)}^{\alpha }}\hfill \\ & \ge \left(\frac{\lambda \theta {(u-t)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)+\lambda }\right){𝟙}_{u-t\le 1}+\left(\frac{\theta \lambda }{\mathsf{\Gamma }(\alpha +1)+\lambda }\right){𝟙}_{u-t>1}\hfill \end{array}$$

(25)

Again, for $T-s\le 1$, we compute ${\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du$ as

$$\begin{array}{cc}\hfill {\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du& \ge {\int }_s^T\frac{\lambda \theta {(u-s)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)+\lambda }du\hfill \\ \hfill & =\frac{\lambda \theta {(T-s)}^{\alpha +1}}{\mathsf{\Gamma }(\alpha +2)+\lambda (\alpha +1)}\hfill \end{array}$$

(26)

For $T-s>1$,

$$\begin{array}{cc}\hfill {\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du& \ge {\int }_s^{s+1}\frac{\lambda \theta {(u-s)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)+\lambda }du+{\int }_{s+1}^T\frac{\theta \lambda }{\mathsf{\Gamma }(\alpha +1)+\lambda }du\hfill \\ & =\frac{\lambda \theta (1+(\alpha +1\left)\right(T-s-1\left)\right)}{\mathsf{\Gamma }(\alpha +2)+\lambda (\alpha +1)}\hfill \\ & \ge \frac{\lambda \theta }{\mathsf{\Gamma }(\alpha +2)+\lambda (\alpha +1)}\hfill \end{array}$$

(27)

and the result follows.□

Lemma 2.

Let $0\le t\le s\le T$. Then, for every $n\ge 0$, there exists $C=C(n,\rho )$ such that

$$|{\partial }_x^{n}G(X_s,w_s)|\le C{\left({\int }_s^T{\mathbb{E}}_s\left({\sigma }_u^2\right)du\right)}^{-\frac{1}{2}(n+1)}$$

(28)

where $G(X_s,w_s)=({\partial }_{xx}^2-{\partial }_x)B(X_s,w_s)$

Proof. 

See Lemma 2.1 in [18] or similar result in Lemma 4.1 of [3].□

Proposition 1.

Assuming the volatility model (2), the conditional expectation is

$${\mathbb{E}}_0\left({\sigma }_t^2\right)={\sigma }_0^2-({\sigma }_0^2-\theta ){\int }_0^t{g}^{\alpha ,\lambda }\left(s\right)ds$$

(29)

where

$$g^{\alpha ,\lambda }\left(t\right):=\lambda t^{\alpha -1}{E}_{\alpha ,\alpha }(-\lambda t^{\alpha })$$

(30)

with $E_{\alpha ,\beta }\left(x\right)={\sum }_{n=0}^{\infty }\frac{x^n}{\mathsf{\Gamma }(\alpha n+\beta )}$ as the generalized Mittag–Leffler function.

Proof. 

See pages 7–8 in [22], Proposition 3.1 in [23] or Proposition 2.1 in [24]. □

2. Approximation for Option Pricing Formula

In this section, we prove that the errors for approximation formula of rough Heston model in the range of $0.1\le H<0.5$ is bounded.

Theorem 2

(Approximation formula). Assume model (3) and (4) with Assumption 1 satisfied. Then, for all $t\in [0,1]$ and $\alpha \in [0.6,1)$, we have

$$\begin{array}{cc}\hfill & \left|{\mathbb{E}}_t[B(X_T,w_T\left(T\right)]-B(X_t,w_t\left(T\right))-H(X_t,w_t\left(T\right))U_t-K(X_t,w_t\left(T\right))R_t\right|\hfill \\ \hfill & \phantom{V_t}\le C(\nu ,\rho ,\alpha ,\lambda ,\theta ){(T-t)}^{\frac{5}{2}\alpha -\frac{3}{2}}\hfill \end{array}$$

(31)

and for $t\in (1,\infty )$ and $\alpha \in (0.5,1)$, we have

$$\begin{array}{cc}\hfill & \left|{\mathbb{E}}_t[B(X_T,w_T\left(T\right)]-B(X_t,w_t\left(T\right))-H(X_t,w_t\left(T\right))U_t-K(X_t,w_t\left(T\right))R_t\right|\hfill \\ \hfill & \phantom{V_t}\le C_1(\nu ,\rho ,\alpha ,\lambda ,\theta )+C_2(\nu ,\rho ,\alpha ,\lambda ,\theta ){(T-t)}^{4\alpha +1}\hfill \end{array}$$

(32)

where $C(\nu ,\rho ,\alpha ,\lambda ,\theta )$, $C_1(\nu ,\rho ,\alpha ,\lambda ,\theta )$ and $C_2(\nu ,\rho ,\alpha ,\lambda ,\theta )$ are functions depending on the parameters $\nu ,\rho ,\alpha ,\lambda ,\theta$, but not T and t. The terms $H(X_t,w_t\left(T\right))$, $K(X_t,w_t\left(T\right))$, $U_t$ and $R_t$ are defined as

$$\begin{array}{cc}\hfill H(X_t,w_t\left(T\right))& =\frac{{\partial }^2}{\partial xw}B(X_t,w_t\left(T\right))\hfill \end{array}$$

$$\begin{array}{cc}& =\frac{1}{2}\left(\frac{{\partial }^3}{\partial x^3}-\frac{{\partial }^2}{\partial x^2}\right)B(X_t,w_t\left(T\right))\hfill \\ \hfill K(X_t,w_t\left(T\right))& =\frac{{\partial }^2}{\partial w^2}B(X_t,w_t\left(T\right))\hfill \end{array}$$

(33)

$$\begin{array}{cc}& =\frac{1}{4}\left(\frac{{\partial }^4}{\partial x^3}-2\frac{{\partial }^3}{\partial x^3}+\frac{{\partial }^2}{\partial x^2}\right)B(X_t,w_t\left(T\right))\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill U_t& ={\mathbb{E}}_t\left[{\int }_t^{T}d{\langle X,M\rangle }_s\right]\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill R_t& =\frac{1}{2}{\mathbb{E}}_t\left[{\int }_t^{T}d{\langle M,M\rangle }_s\right]\hfill \end{array}$$

(36)

Proof. 

See Appendix A.□

Remark 1.

Theorem 2 shows that the error bound does not diverge for any $\alpha \in [0.6,1]$ (provided that $\lambda >0$ and $\theta >0$). Coincidentally, for $t\in [0,1]$, we can obtain the same order convergence when $\alpha =1$ in the case of classical Heston model as shown in [18,19]. Furthermore, if Assumption 1 is satisfied, then the errors of the approximation formula can be effectively reduced when ν and $T-t$ are small.

3. Second-Order Approximation and Small-Time Behavior of Implied Volatility

The outline of this section is roughly the same as that of Alos, Elisa and De Santiago, Rafael and Vives, Josep [19], but we focus on the rough Heston model rather than classical Heston model. We expand the function f with respect to the two scales in order to obtain the second-order approximation for the implied volatility. Consider the asymptotic sequence ${\left\{{\delta }^i\right\}}_{i=0}^{\infty }$ and the function f can be expanded as $f=f_0^{ϵ}+\delta f_1^{ϵ}+{\delta }^2{f}_2^{ϵ}+\dots$. Then, consider the scale function ${\left\{{ϵ}^i\right\}}_{i=0}^{\infty }$ and expand each $f_i^{ϵ}$ as $f_i^{ϵ}=f_{0,i}+epsilon{f}_{1,i}+{ϵ}^2{f}_{2,i}+\dots$. Combining the expansions $f_i^{ϵ}$ and f, we obtain the following

$$f=f_{0,0}+ϵf_{1,0}+\delta f_{0,1}+O({(ϵ+\delta )}^2)$$

(37)

We thus let $ϵ=\rho \nu /\alpha$ and $\delta ={\nu }^2/\alpha$ with $I(T,K)$ be the implied volatility function expanded to these two scales:

$$I(T,K)={\nu }_0+\frac{\rho \nu }{\alpha }{I}_1(T,K)+\frac{{\nu }^2}{\alpha }{I}_2(T,K)+O({(\frac{\rho \nu +{\nu }^2}{\alpha })}^2)$$

(38)

Let

$$\widehat{I}(T,K)={\nu }_0+\frac{\rho \nu }{\alpha }{I}_1(T,K)+\frac{{\nu }^2}{\alpha }{I}_2(T,K)$$

(39)

as the second-order approximation to the implied volatility. Then, let $v_t=\sqrt{\frac{1}{T}{w}_t\left(T\right)}$ as the future average volatility and $V(X_T,v_t)=B(X_T,w_t\left(T\right))$ such that $\widehat{V}(X_T,v_t\left(T\right))$ is the approximation formula to the option price obtained in Theorem 2. Ultimately, we have

$$\begin{array}{cc}\hfill \widehat{V}(X_T,v_T)& :=V(X_0,v_0)+H(X_0,v_0^{2}T)U_0+K(X_0,v_0^{2}T)R_0\hfill \\ & =V(X_0,v_0)+\frac{e^x{e}^{-\frac{d_{+}^2\left(v_0^{2}T\right)}{2}}}{2v_0\sqrt{2\pi T}}\left(1-\frac{d_{+}\left(v_0^{2}T\right)}{v_0\sqrt{T}}\right)U_0\hfill \\ & \phantom{=}+\frac{e^x{e}^{-\frac{d_{+}^2\left(v_0^{2}T\right)}{2}}}{4v_0\sqrt{2\pi T}}\left(\frac{d_{+}^2\left(v_0^{2}T\right)-v_0\sqrt{T}{d}_{+}\left(v_0^{2}T\right)-1}{v_0^{2}T}\right)R_0\hfill \\ & =V(X_0,v_0)+\frac{1}{2v_{0}T}\frac{\partial V(X_0,v_0)}{\partial \sigma }\left(1-\frac{d_{+}\left(v_0^{2}T\right)}{v_0\sqrt{T}}\right)U_0\hfill \\ & \phantom{=}+\frac{1}{4v_{0}T}\frac{\partial V(X_0,v_0)}{\partial \sigma }\left(\frac{d_{+}^2\left(v_0^{2}T\right)-v_0\sqrt{T}{d}_{+}\left(v_0^{2}T\right)-1}{v_0^{2}T}\right)R_0\hfill \end{array}$$

(40)

where $\frac{\partial V(X_0,v_0)}{\partial \sigma }=\frac{e^x{e}^{-\frac{d_{+}^2\left(v_0^{2}T\right)}{2}}\sqrt{T}}{\sqrt{2\pi }}$ is the Vega of the Black–Scholes formula $V(X_0,v_0)$. Similarly, we can expand the $V(X_0,v_0)$ around $v_0$ using the Taylor Expansion and Equation (37) as such

$$\begin{array}{cc}\hfill V(X_0,I(T,K))& =V(X_0,v_0)+\frac{\partial V(X_0,v_0)}{\partial \sigma }\left(\frac{\rho \nu }{\alpha }{I}_1(T,K)+\frac{{\nu }^2}{\alpha }{I}_2(T,K)+\dots \right)\hfill \\ & \phantom{=}+\frac{1}{2}\frac{{\partial }^{2}V(X_0,v_0)}{\partial {\sigma }^2}{\left(\frac{\rho \nu }{\alpha }{I}_1(T,K)+\frac{{\nu }^2}{\alpha }{I}_2(T,K)+\dots \right)}^2+\dots \hfill \\ & =V(X_0,v_0)+\frac{\rho \nu }{\alpha }\frac{\partial V(X_0,v_0)}{\partial \sigma }{I}_1(T,K)+\frac{{\nu }^2}{\alpha }\frac{\partial V(X_0,v_0)}{\partial \sigma }{I}_2(T,K)+.\dots \hfill \end{array}$$

(41)

Let ${\widehat{I}}_1(T,K):=\frac{\rho \nu }{\alpha }{I}_1(T,K)$ and ${\widehat{I}}_2(T,K):=\frac{{\nu }^2}{\alpha }{I}_2(T,K)$; then, comparing the first-order approximation of Equation (41) to Equation (40), we can obtain the following

$$\begin{array}{cc}\hfill {\widehat{I}}_1(T,K)& =\frac{1}{2v_{0}T}\left(1-\frac{d_{+}\left(v_0^{2}T\right)}{v_0\sqrt{T}}\right)U_0\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\widehat{I}}_2(T,K)& =\frac{1}{4v_{0}T}\left(\frac{d_{+}^2\left(v_0^{2}T\right)-v_0\sqrt{T}{d}_{+}\left(v_0^{2}T\right)-1}{v_0^{2}T}\right)R_0\hfill \end{array}$$

(43)

Similarly, as noted by Alos, Elisa and De Santiago, Rafael and Vives, Josep [19], the smile effect of the implied volatility can be shown by the quadratic term $d_{+}^2\left(v_0^{2}T\right)$ which has impact of the ${(X_t-lnK)}^2$ in its term. Suppose we let $\epsilon =V(X_0,v_0)-\widehat{V}(X_0,v_0)$, and we let $I^{*}={\widehat{I}}_1(T,K)+{\widehat{I}}_2(T,K)$. The error between the option pricing formula using the approximation of the implied volatility in Equation (38) and the option price under rough Heston model is as follows:

$$\begin{array}{cc}\hfill \left|V(X_0,\widehat{I}(T,K))-V(X_0,v_0)\right|& =V(X_0,v_0)+\frac{\partial V(X_0,v_0)}{\partial \sigma }{I}^{*}+\frac{1}{2}\frac{{\partial }^{2}V(X_0,v_0)}{\partial {\sigma }^2}{\left(I^{*}\right)}^2+\cdots \hfill \\ & \phantom{=}-\left(V(X_0,v_0)+H(X_0,v_0^{2}T)U_0+K(X_0,v_0^{2}T)R_0\right)\hfill \\ & \le \left|\sum _{i=2}^{\infty }\frac{1}{i!}\frac{{\partial }^{i}V(X_0,v_0)}{\partial {\sigma }^i}{\left(I^{*}\right)}^i\right|+\left|\epsilon \right|\hfill \end{array}$$

(44)

The first term of Equation (44) has an error of $\mathcal{O}({\left(I^{*}\right)}^2)$ or $\mathcal{O}\left({\left(\frac{\rho \nu +{\nu }^2}{\alpha }\right)}^2\right)$. Under Assumption 1, and by Theorem 2, we know that the second term of (44) has a small time (when $0\le T\le 1$) error of $\mathcal{O}\left(C\right(\nu ,\rho ,\alpha ,\lambda ,\theta \left)\right)$ which is a decreasing term proportional to $\rho$ and $\nu$. The total errors of option prices on the approximation to the implied volatility are of the order of $\mathcal{O}\left({\left(\frac{\rho \nu +{\nu }^2}{\alpha }\right)}^2+C(\nu ,\rho ,\alpha ,\lambda ,\theta )\right)$. Noticeably, the errors converges to zero as $\rho \to 0$ and $\nu \to 0$. We prove some results that would assist in obtaining the short-time behavior on the approximation of implied volatility.

Proposition 2.

Assume volatility dynamics of Equation (4); then,

$$\underset{T\to 0}{lim}{v}_0={\sigma }_0$$

(45)

Proof. 

From the definition of $v_0$ and Proposition 1, we can obtain the following

$$\begin{array}{cc}\hfill v_0& :=\sqrt{\frac{1}{T}{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds}\hfill \\ & =\sqrt{\frac{1}{T}{\int }_0^T{\sigma }_0^2-({\sigma }_0^2-\theta ){\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)duds}\hfill \\ & =\sqrt{{\sigma }_0^2-\frac{{\sigma }_0^2-\theta }{T}{\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)duds}\hfill \end{array}$$

(46)

Let ${lim}_{T\to 0}{v}_0$; we can subsequently notice that it is sufficient to prove that

$$\underset{T\to 0}{lim}\frac{1}{T}{\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)duds=0$$

(47)

Note that

$$\begin{array}{cc}\hfill \sum _{n=0}^{\infty }{\int }_0^s\left|\frac{{(-\lambda )}^n{u}^{n\alpha +\alpha -1}}{\mathsf{\Gamma }(\alpha n+\alpha )}\right|du& =\sum _{n=0}^{\infty }\frac{{\lambda }^n{u}^{n\alpha +\alpha }}{\mathsf{\Gamma }(\alpha n+\alpha +1)}\hfill \\ & =\lambda \sum _{n=0}^{\infty }\frac{{\lambda }^{n+1}{u}^{\alpha (n+1)}}{\mathsf{\Gamma }\left(\alpha \right(n+1)+1)}\hfill \\ & =\lambda (-1+E_{\alpha ,1}\left(\lambda u^{\alpha }\right))\hfill \\ & <\infty \hfill \end{array}$$

(48)

Then, by the above result and Fubini’s theorem, we can obtain the following

$$\begin{array}{cc}\hfill \underset{T\to 0}{lim}\frac{1}{T}{\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)duds& =\lambda \underset{T\to 0}{lim}\frac{1}{T}{\int }_0^T{\int }_0^s\sum _{n=0}^{\infty }\frac{{(-\lambda )}^n{u}^{n\alpha +\alpha -1}}{\mathsf{\Gamma }(\alpha n+\alpha )}duds\hfill \\ & =\lambda \underset{T\to 0}{lim}\frac{1}{T}\sum _{n=0}^{\infty }{\int }_0^T{\int }_0^s\frac{{(-\lambda )}^n{u}^{n\alpha +\alpha -1}}{\mathsf{\Gamma }(\alpha n+\alpha )}duds\hfill \\ & =\lambda \underset{T\to 0}{lim}\frac{1}{T}\sum _{n=0}^{\infty }{\int }_0^T\frac{{(-\lambda )}^n{s}^{n\alpha +\alpha }}{\mathsf{\Gamma }(\alpha n+\alpha +1)}ds\hfill \\ & =\lambda \underset{T\to 0}{lim}\sum _{n=0}^{\infty }\frac{{(-\lambda )}^n{T}^{n\alpha +\alpha }}{\mathsf{\Gamma }(\alpha n+\alpha +2)}\hfill \\ & =0\hfill \end{array}$$

(49)

which is the desired result.□

Proposition 3.

Assume volatility dynamics of Equation (4); for $k>-\alpha -1$ and $k\ne -1$,

$${\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right){(T-s)}^{k}ds\sim {\sigma }_0^2{\int }_0^T{(T-s)}^{k}ds$$

(50)

as $T\to 0$.

Proof. 

From Proposition 1, we can obtain the following

$$\begin{array}{cc}\hfill {\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right){(T-s)}^{k}ds& ={\int }_0^T\left[{\sigma }_0^2-({\sigma }_0^2-\theta ){\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)du\right]{(T-s)}^{k}ds\hfill \\ \hfill & ={\sigma }_0^2{\int }_0^T{(T-s)}^{k}ds-({\sigma }_0^2-\theta ){\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)du{(T-s)}^{k}ds\hfill \end{array}$$

(51)

From Proposition 2 and by the Fubini’s theorem again, we can obtain the following:

$$\begin{array}{cc}\hfill {\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)du{(T-s)}^{k}ds& =\lambda {\int }_0^T\sum _{n=0}^{\infty }\frac{{(-\lambda )}^n{s}^{n\alpha +\alpha }}{\mathsf{\Gamma }(\alpha n+\alpha +1)}{(T-s)}^{k}ds\hfill \\ & =\lambda \sum _{n=0}^{\infty }{\int }_0^T\frac{{(-\lambda )}^n{s}^{n\alpha +\alpha }}{\mathsf{\Gamma }(\alpha n+\alpha +1)}{(T-s)}^{k}ds\hfill \\ & =\lambda \sum _{n=0}^{\infty }\frac{{(-\lambda )}^n}{\mathsf{\Gamma }(\alpha n+\alpha +1)}{\int }_0^T{s}^{n\alpha +\alpha }{(T-s)}^{k}ds\hfill \end{array}$$

(52)

Substitute $y=s/T$ and notice that it can be solved using the Beta function, such that

$${\int }_0^T{s}^{n\alpha +\alpha }{(T-s)}^{k}ds=T^{n\alpha +\alpha +k+1}\frac{\mathsf{\Gamma }(n\alpha +\alpha +1)\mathsf{\Gamma }(k+1)}{\mathsf{\Gamma }(n\alpha +\alpha +k+2)}$$

(53)

We then have

$${\int }_0^T{\int }_0^s{g}^{\alpha ,\lambda }\left(u\right)du{(T-s)}^{k}ds=\lambda \sum _{n=0}^{\infty }\frac{{(-\lambda )}^n\mathsf{\Gamma }(k+1)T^{n\alpha +\alpha +k+1}}{\mathsf{\Gamma }(n\alpha +\alpha +k+2)}$$

(54)

Then, for $k>-\alpha -1$ and $k\ne -1$, we can deduce the following asymptotic expansion

$${\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right){(T-s)}^{k}ds\sim {\sigma }_0^2{\int }_0^T{(T-s)}^{k}ds$$

(55)

as $T\to 0$.□

Remark 2.

Proposition 2 is equivalent to the case of the classical Heston model. Furthermore, by utilizing Proposition 3, we ultimately yield the same result from [19] in the case of $\alpha =1$ for the rest of the result in this paper.

3.1. Limiting Behavior of $\widehat{I}(T,K)$ When $T\to 0$

We utilize some of the derivations from [19]. The limiting behavior of $\widehat{I}(T,K)$ as $T\to 0$ can be computed in term of the asymptotic formula as follows

Lemma 3.

Assume model (3) and (4). Then, we can obtain the asymptotic formula of ${\widehat{I}}_1(T,K)$ as

$${\widehat{I}}_1(T,K)\sim -\frac{\rho \nu (X_0-lnK)}{2\mathsf{\Gamma }(\alpha +2){\sigma }_0}{T}^{\alpha -1}$$

(56)

as $T\to 0$.

Proof. 

From Equation (A7) and Proposition 3, we can obtain the limiting behavior of $U_0$ when $T\to 0$ as

$$\begin{array}{cc}\hfill U_0& \sim \frac{\nu \rho {\sigma }_0^2}{\mathsf{\Gamma }(\alpha +1)}{\int }_0^T{(T-s)}^{\alpha }ds\hfill \\ \hfill & \sim \frac{\nu \rho {\sigma }_0^2}{\mathsf{\Gamma }(\alpha +2)}{T}^{\alpha +1}\hfill \end{array}$$

(57)

Then, the limiting behavior of ${\widehat{I}}_1(T,K)$ can be obtained as following

$$\begin{array}{cc}\hfill \underset{T\to 0}{lim}{\widehat{I}}_1(T,K)& =\underset{T\to 0}{lim}\frac{1}{2v_{0}T}\left(1-\frac{d_{+}\left(v_0^{2}T\right)}{v_0\sqrt{T}}\right)U_0\hfill \\ \hfill & =\underset{T\to 0}{lim}\frac{1}{4v_{0}T}{U}_0-\underset{T\to 0}{lim}\frac{(X_0-lnK)}{2v_0^3{T}^2}{U}_0\hfill \end{array}$$

(58)

The result follows from the limiting behavior of Equation (57) and Proposition 2.

We focus on the asymptotic formula of ${\widehat{I}}_2(T,K)$ in the next lemma.□

Lemma 4.

Assume model (3) and (4). Then, we can obtain the asymptotic formula of ${\widehat{I}}_2(T,K)$ as

$${\widehat{I}}_2(T,K)\sim \frac{{\nu }^2{(X_0-lnK)}^2}{8{\sigma }_0^3\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha -2}$$

(59)

as $T\to 0$.

Proof. 

From Equation (A7) and Proposition 3, the limiting behavior of $R_0$ when $T\to 0$ can be obtained as

$$\begin{array}{cc}\hfill R_0& \sim \frac{{\nu }^2{\sigma }_0^2}{2\mathsf{\Gamma }{(\alpha +1)}^2}{\int }_0^T{(T-s)}^{2\alpha }ds\hfill \\ \hfill & \sim \frac{{\nu }^2{\sigma }_0^2}{2\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha +1}\hfill \end{array}$$

(60)

Then, the limiting behavior of ${\widehat{I}}_2(T,K)$ can be obtained as follows

$$\begin{array}{cc}\hfill \underset{T\to 0}{lim}{\widehat{I}}_2(T,K)& =\underset{T\to 0}{lim}\frac{1}{4v_{0}T}\left[\frac{d_{+}^2\left(v_0^{2}T\right)-v_0\sqrt{T}{d}_{+}\left(v_0^{2}T\right)-1}{v_0^{2}T}\right]R_0\hfill \\ \hfill & =\underset{T\to 0}{lim}\left[\frac{{(X_0-lnK)}^2}{4v_0^5{T}^3}-\frac{1}{16v_{0}T}-\frac{1}{4v_0^3{T}^2}\right]R_0\hfill \end{array}$$

(61)

Similarly, the result follows from the limiting behavior of Equation (60) and Proposition 2.□

Remark 3.

From Lemmas 3 and 4, when the option is close to maturity $(T\to 0)$, the second-order approximation to the implied volatility can be written

$$\widehat{I}(T,K)\sim {\sigma }_0-\frac{\rho \nu (X_0-lnK)}{2\mathsf{\Gamma }(\alpha +2){\sigma }_0}{T}^{\alpha -1}+\frac{{\nu }^2{(X_0-lnK)}^2}{8{\sigma }_0^3\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha -2}$$

(62)

as the asymptotic formula. Furthermore, the expression (62) becomes exact and equivalent to the result found by Alos, Elisa and De Santiago, Rafael and Vives, Josep [19] when $\alpha =1$.

3.2. Term Structure of At-The-Money Implied Volatility Skew

This subsection is crucial in showing that the approximation formula for the implied volatility under rough Heston model is capable of replicating the explosive behavior of the implied volatility under the rough volatility model as stated in [4,5]. When the option is at-the-money, we have $lnK=X_0$ and this leads to $d_{+}\left(v_0^{2}T\right)=\frac{v_0}{2}\sqrt{T}$. Consequently, the term structure of the at-the-money implied volatility skew can be computed as

Lemma 5.

When the option is at-the-money, the term structure of implied volatility skew has the asymptotic formula of

$$\frac{\partial }{\partial lnK}\widehat{I}(T,K){|}_{lnK=X_0}\sim \frac{\rho \nu }{2\mathsf{\Gamma }(\alpha +2){\sigma }_0}{T}^{\alpha -1}$$

(63)

as $T\to 0$.

Proof. 

From the result of Lemmas 3 and 4, we have the following as $T\to 0$:

$$\begin{array}{cc}\hfill \widehat{I}(T,K)& \sim {\sigma }_0-\frac{\rho \nu (X_0-lnK)}{2\mathsf{\Gamma }(\alpha +2){\sigma }_0}{T}^{\alpha -1}+\frac{{\nu }^2{(X_0-lnK)}^2}{8{\sigma }_0^3\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha -2}\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial lnK}\widehat{I}(T,K)& \sim \frac{\rho \nu }{2\mathsf{\Gamma }(\alpha +2){\sigma }_0}{T}^{\alpha -1}-\frac{{\nu }^2(X_0-lnK)}{4{\sigma }_0^3\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha -2}\hfill \end{array}$$

(65)

Let $lnK=X_0$ and the result follows.□

Remark 4.

From Lemma 5, we can also observe the following: for the uncorrelated case ($\rho =0$) and correlated case $(\rho \ne 0)$ in model (3) and (4), the slope of implied volatility (with respect to $lnK$ and $T\to 0$) is positive when $\ln K>X_0$ and negative when $lnK<X_0$. The result is slightly different than the implied volatility behavior when $\alpha =1$, where $\rho \ne 0$ yields different positive and negative region (slope is negative when $lnK<X_0-3\rho {\sigma }_0^2/\nu$ and positive when $lnK>X_0-3\rho {\sigma }_0^2/\nu$). The empirical result in Figure 8 (page 11) of [8] shows that Lemma 5 is valid as $T\to 0$. On a side note, Lemma 5 confirms that the term structure of at-the-money approximation of implied volatility under rough Heston model has explosive behavior of order $T^{\alpha -1}$, i.e., it possesses an important feature which the classical Heston model does not possess. The presence of explosive behavior in the term structure of at-the-money implied volatility skew fits the empirical observation from [4,5].

3.3. Limiting Behavior of Rate of Change of $\widehat{I}(T,K)$ with Respect to T When the Option Is At-The-Money

In this section, we discuss the behavior of the rate of change of implied volatility $\widehat{I}(T,K)$ with respect to T when $T\to 0$ and the option is at-the-money.

Lemma 6.

When the option is at-the-money, the skew of the approximation of implied volatility at $T\to 0$ has the asymptotic formula of

$$\frac{\widehat{I}(T,K)-{\sigma }_0}{T}\sim \frac{\lambda (\theta -{\sigma }_0^2)\mathsf{\Gamma }(\alpha +2)+\rho \nu {\sigma }_0^2\mathsf{\Gamma }(\alpha +1)}{4{\sigma }_0\mathsf{\Gamma }(\alpha +1)\mathsf{\Gamma }(\alpha +2)}{T}^{\alpha -1}-\frac{{\nu }^2}{8{\sigma }_0\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{\left(T^{\alpha -1}\right)}^2$$

(66)

Proof. 

Approximation to implied volatility when the option is at-the-money can be formulated as

$$\begin{array}{cc}\hfill \widehat{I}(T,K)& =v_0+\frac{1}{2v_{0}T}\left(1-\frac{d_{+}\left(v_{0}T\right)}{v_0\sqrt{T}}\right)U_0+\frac{1}{4v_{0}T}\left(\frac{d_{+}^2\left(v_0^{2}T\right)-v_0\sqrt{T}{d}_{+}\left(v_0^{2}T\right)-1}{v_0^{2}T}\right)R_0\hfill \\ \hfill & =v_0+\frac{1}{4v_{0}T}{U}_0-\frac{1}{4v_{0}T}\left(\frac{1}{4}+\frac{1}{v_0^{2}T}\right)R_0\hfill \end{array}$$

(67)

From Lemmas 3 and 4 and Remark 3, it is easy to see that, when the option is at-the-money ($lnK=X_0$), $\widehat{I}(0,K)={\sigma }_0$. To find the rate of change of $\widehat{I}(T,K)$ with respect to T, we use the definition of derivative as follows

$$\begin{array}{cc}\hfill \underset{T\to 0}{lim}\frac{\widehat{I}(T,K)-{\sigma }_0}{T}& =\underset{T\to 0}{lim}\left[\frac{v_0-{\sigma }_0}{T}+\frac{1}{4v_0{T}^2}{U}_0-\frac{1}{4v_0{T}^2}\left(\frac{1}{4}+\frac{1}{v_0^{2}T}\right)R_0\right]\hfill \end{array}$$

(68)

We solve the above expression separately. The first term of the above expression can be computed using L’Hopital’s rule repetitively as follows

$$\begin{array}{cc}\hfill \underset{T\to 0}{lim}\frac{v_0-{\sigma }_0}{T}& =\underset{T\to 0}{lim}\frac{\sqrt{\frac{1}{T}{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds}-{\sigma }_0}{T}\hfill \\ & =\underset{T\to 0}{lim}\frac{-\frac{1}{T^2}{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds+\frac{1}{T}{\mathbb{E}}_0\left({\sigma }_T^2\right)}{2\sqrt{\frac{1}{T}{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds}}\hfill \\ & =\underset{T\to 0}{lim}\frac{-\frac{1}{T^2}\left[{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds-T{\mathbb{E}}_0\left({\sigma }_T^2\right)\right]}{2\sqrt{\frac{1}{T}{\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds}}\hfill \\ & =\underset{T\to 0}{lim}-\frac{1}{2{\sigma }_0{T}^2}\left({\int }_0^T{\mathbb{E}}_0\left({\sigma }_s^2\right)ds-T{\mathbb{E}}_0\left({\sigma }_T^2\right)\right)\hfill \\ & =\underset{T\to 0}{lim}-\frac{1}{4{\sigma }_{0}T}\left({\mathbb{E}}_0\left({\sigma }_T^2\right)-\left[T\frac{\partial }{\partial T}{\mathbb{E}}_0\left({\sigma }_T^2\right)+{\mathbb{E}}_0\left({\sigma }_T^2\right)\right]\right)\hfill \\ & =\underset{T\to 0}{lim}\frac{1}{4{\sigma }_0}\frac{\partial }{\partial T}{\mathbb{E}}_0\left({\sigma }_T^2\right)\hfill \end{array}$$

(69)

We have ${\mathbb{E}}_0\left({\sigma }_T^2\right)={\sigma }_0^2+(1/\mathsf{\Gamma }\left(\alpha \right)){\int }_0^T{(T-s)}^{\alpha -1}\lambda (\theta -{\mathbb{E}}_0\left({\sigma }_s^2\right))ds$. Then, it follows that

$$\underset{T\to 0}{lim}\frac{{\mathbb{E}}_0\left({\sigma }_T^2\right)-{\sigma }_0^2}{T}=\underset{T\to 0}{lim}\frac{\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^T{(T-s)}^{\alpha -1}\lambda (\theta -{\mathbb{E}}_0\left({\sigma }_s^2\right))ds}{T}$$

(70)

By Proposition 3, we have the following

$$\begin{array}{cc}\hfill \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^T{(T-s)}^{\alpha -1}\lambda (\theta -{\mathbb{E}}_0\left({\sigma }_s^2\right))ds& \sim \frac{\lambda (\theta -{\sigma }_0^2)}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^T{(T-s)}^{\alpha -1}ds\hfill \end{array}$$

(71)

$$\begin{array}{cc}& \sim \frac{\lambda (\theta -{\sigma }_0^2)}{\mathsf{\Gamma }(\alpha +1)}{T}^{\alpha }\hfill \end{array}$$

(72)

as $T\to 0$. Hence,

$$\frac{{\mathbb{E}}_0\left({\sigma }_T^2\right)-{\sigma }_0^2}{T}\sim \frac{\lambda (\theta -{\sigma }_0^2)}{\mathsf{\Gamma }(\alpha +1)}{T}^{\alpha -1}$$

(73)

as $T\to 0$. Finally,

$$\frac{v_0-{\sigma }_0}{T}\sim \frac{\lambda (\theta -{\sigma }_0^2)}{4{\sigma }_0\mathsf{\Gamma }(\alpha +1)}{T}^{\alpha -1}$$

(74)

as $T\to 0$. From the previous result in Lemma 3, we know that $U_0\sim \frac{\nu \rho {\sigma }_0^2}{\mathsf{\Gamma }(\alpha +2)}{T}^{\alpha +1}$ as $T\to 0$; then, along with Proposition 2, the second term of Equation (63) can be computed in terms of asymptotic formula

$$\frac{1}{4v_0{T}^2}{U}_0\sim \frac{\rho \nu {\sigma }_0}{4\mathsf{\Gamma }(\alpha +2)}{T}^{\alpha -1}$$

(75)

as $T\to 0$. Similarly, by the previous result in Lemma 3 such that $R_0\sim \frac{{\nu }^2{\sigma }_0^2}{2\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha +1}$ as $T\to 0$ and Proposition 2, the third term of Equation (63) can be computed as

$$\left(-\frac{1}{16v_0{T}^2}-\frac{1}{4v_0^3{T}^3}\right)R_0\sim -\frac{{\nu }^2{\sigma }_0}{32\mathsf{\Gamma }{(\alpha +1)}^2(2\alpha +1)}{T}^{2\alpha -1}$$

(76)

Then, the result follows with the sum of three computed terms (74)–(76). □

Remark 5.

Lemma 6 satisfies the result of $\alpha =1$ from the work of Alos, Elisa and De Santiago, Rafael and Vives, Josep [19]. Furthermore, when the option is at-the-money, we can obtain a rough approximation of implied volatility in terms of T as contrary to Remark 3 which provides constant term when the option is at-the-money. An integration of the result in Equation (66) with respect to T near $T=0$ gives us the following approximation:

$$\widehat{I}(T,K)\approx {\sigma }_0+\frac{\lambda (\theta -{\sigma }_0^2)\mathsf{\Gamma }(\alpha +2)+\rho \nu {\sigma }_0^2\mathsf{\Gamma }(\alpha +1)}{4{\sigma }_0\alpha \mathsf{\Gamma }(\alpha +1)\mathsf{\Gamma }(\alpha +2)}{T}^{\alpha }-\frac{{\nu }^2}{8{\sigma }_0\mathsf{\Gamma }{(\alpha +1)}^2(4{\alpha }^2-1)}{T}^{2\alpha -1}$$

(77)

4. Numerical Experiments

Through numerical experiments, we can show that the approximation formula for the rough Heston model possesses the explosive behavior of at-the-money implied volatility skew, which is desirable in a stochastic volatility model to match the empirical implied volatility data in the market. Although Theorem 2 shows that, under Assumption 1, the approximation formula for the rough Heston model is bounded from the rough Heston model using the same set of parameters, we decided to conduct the numerical experiments using calibration approach, i.e., calibrate a new set of parameters used on approximation formula based on the artificial data given by rough Heston model. We wish to note that the components or terms of $v_0,U_0$ and $R_0$ can be obtained using Proposition 1. The following parameters for the rough Heston model are set:

$$\begin{array}{cc}\hfill {\sigma }_0^2& =0.04,\theta =0.04,\alpha =0.6,\hfill \\ \hfill \lambda & =1,\phantom{xxl}\nu =0.1,\phantom{l}\rho =-0.65\hfill \end{array}$$

(78)

Then, we calibrate the approximation Formula (40) to the rough Heston model through minimizing the squared errors of implied volatility and subsequently we obtain the following parameters:

$$\begin{array}{cc}\hfill {\sigma }_0^{2*}& =0.039671,{\theta }^{*}=0.039230,{\alpha }^{*}=0.568425,\hfill \\ \hfill {\lambda }^{*}& =1.942492,{\nu }^{*}=0.058033,{\rho }^{*}=-0.911687\hfill \end{array}$$

(79)

Figure 1 shows a calibrated option prices from $T=0.01$ to $T=0.1$. After the calibration of the parameters, the approximation formula for option prices under the rough Heston model are capable of displaying a great fit to the option prices computed using Fourier inversion method with the fractional Adams scheme [21], as displayed in [11].

The implied volatility plots are shown in Figure 2. In the plot, we can observe that the implied volatilities computed using the approximation formula (40) match the ones produced by the Fourier inversion method with the fractional Adams scheme [11,21] under the rough Heston model very well, except for at $T=0.01$ where the implied volatility is understated when the option is deep out-of-money and in-the-money.

One result that sets apart from many other stochastic volatility and jump diffusion models is the short-time behavior of term-structure at-the-money implied volatility skew, i.e., derivative of implied volatility with respect to log-strike for the at-the-money call option. In Figure 3, we display the term structure of the at-the-money implied volatility skew using three different methods—Fourier inversion method with fractional Adams scheme [11,21], approximation formula for option prices (40) (transformed back to implied volatility) and Equation (63). The result is pretty much consistent with the ones demonstrated by the approximation formula for rough Heston model (40) and Equation (63). Explosive behavior on the short-time term-structure of at-the-money implied volatility skew is a remarkable feature on a stochastic volatility model and the behavior is commonly observed in the empirical data in the financial market.

It is important to realize that our work is relying on Assumption 1. In terms of practical application, Assumption 1 will be fulfilled whenever the current variance process ${\sigma }_t^2$ is equivalent or below the mean reversion level $\theta$. In other words, the usual flat forward variance curve (${\mathbb{E}}_0\left({\sigma }_t^2\right)={\mathbb{E}}_0\left({\sigma }_0^2\right)$ for all t) frequently used by many practitioners will satisfy Assumption 1 too. Other than that, we would say that Assumption 1 will be satisfied roughly $50\%$ of the time under the normal market condition.

5. Conclusions

The main contribution of this paper lies within the bounded gaps of option prices on approximation formula for the option prices of the rough Heston model under the simple Assumption 1. Due to the nature of fractional terms in the volatility component of the rough Heston model, the errors bound incurred by the decomposition formula are significantly higher than the classical version of the approximation formula for the decomposition formula of classical Heston model [18,19]. In other words, this disadvantage decreases the usability of the approximation formula in practice due to the errors incurred. Nevertheless, in the numerical experiments, we showed that it is possible to match the option prices and implied volatility quite well by re-calibrating the parameters of the approximation formula for option prices.

In addition, we propose a second-order approximation to the implied volatility of the rough Heston model. Through the approximation of implied volatility of rough Heston model, we manage to prove some important results which includes the explosive behavior of short-time term structure of at-the-money implied volatility skew using the approximation formula of option prices for the rough Heston model, which is desirable in the empirical financial market. Numerical experiments were conducted and we subsequently discovered that the approximation formula for option prices using the rough Heston model does have explosive behavior in its term-structure of at-the-money implied volatility skew and roughly at the same rate as the rough Heston model after recalibrating the parameters.