A Solution to the Time-Scale Fractional Puzzle in the Implied Volatility

In the option pricing literature, it is well known that (i) the decrease in the smile amplitude is much slower than the standard stochastic volatility models and (ii) the term structure of the at-the-money volatility skew is approximated by a power-law function with the exponent close to zero. These stylized facts cannot be captured by standard models, and while (i) has been explained by using a fractional volatility model with Hurst index H > 1 / 2 , (ii) is proven to be satisfied by a rough volatility model with H < 1 / 2 under a risk-neutral measure. This paper provides a solution to this fractional puzzle in the implied volatility. Namely, we construct a two-factor fractional volatility model and develop an approximation formula for European option prices. It is shown through numerical examples that our model can resolve the fractional puzzle, when the correlations between the underlying asset process and the factors of rough volatility and persistence belong to a certain range. More specifically, depending on the three correlation values, the implied volatility surface is classified into four types: (1) the roughness exists, but the persistence does not; (2) the persistence exists, but the roughness does not; (3) both the roughness and the persistence exist; and (4) neither the roughness nor the persistence exist.


Introduction
In the finance literature, there has been a general consensus that volatility is highly persistent. There are numerous pieces of evidence that the price dynamics of financial products are consistent with fractional Brownian motion (fBM) volatility models with Hurst index H > 1/2, which implies that the volatility has a long memory. See, e.g., [1] for the existence of the long memory features in stock market volatilities. However, inconsistent with this stylized fact, [2] recently find that the log-volatility behaves essentially as an fBM with H close to zero at any reasonable time scale, by estimating the volatility from high frequency data. This puzzle (the word "puzzle" is used in the context of using a fractional volatility model in option pricing, but not used in the context of finance in general) has not been resolved, although one possible explanation may be the smoothing effect by sampling intervals of data. Recall that if H = 1/2, the fBM is the standard Brownian motion (BM), which is a Markov process with short memory. Hence, the fact H < 1/2 implies that volatility is rough.
On the other hand, in the context of option pricing, there are also seemingly two inconsistent stylized facts. Namely, (i) the decrease in the smile amplitude is much slower than the standard stochastic volatility models; and (ii) the term structure of the at-the-money volatility skew is well approximated by a power-law function with the exponent close to zero. It is well recognized that the two phenomena cannot be captured by standard models. Nevertheless, in the option pricing literature, these problems have been studied separately.

The Setup
Inspired by the model proposed in [4], we consider the following two-factor fractional volatility model. Namely, the underlying asset price S t and its volatility σ t = σ(X 1 t , X 2 t ) are modeled by the stochastic differential equations (SDEs):

Integral Representation
The fBM w H i t can be represented in terms of the stochastic integral with respect to another standard BM. In this study, we find it useful for the development of our approximation to employ the Mandelbrot-Van Ness representation of fBMs. Namely, we have: where Γ(a) = ∞ 0 t a−1 e −t dt denotes the gamma function and where w i t is a standard BM with constant correlations dw 1 t dw 2 t = ρ 1,2 dt and dw t dw i t = ρ i dt, i = 1, 2. By defining: and: we then have the following representation: Note that, given w i 0 , the future behavior of w i t , t > 0, is independent of the past, because the BM w i t is a Markov process. Hence, supposing that the past w i t , t ≤ 0 has been observed, the quantity g H i t is considered to be a deterministic function of time t ≥ 0 with g H i 0 = 0. The fractional mean-reverting process X i t given in (1) can be solved as: where: By applying the integration-by-parts formula and changing the order of integration, the volatility can be rewritten as: where: See [4] for the detailed derivation. Now, we introduce the market prices of risks η t andη i t to define the standard BMs W t and W i t under a martingale measure Q. Namely, define: respectively. It follows from (1) and the standard argument that the asset price S t follows the SDE: under Q, where r denotes the risk-free spot rate, which we assume to be constant for simplicity.
By applying Ito's formula, we obtain: where F(0, t) = S 0 e (r−q)t is the forward price of the underlying asset with delivery date t. On the other hand, the fractional mean-reverting process X i t is given by: where: Summarizing, our fractional volatility model is formulated as: for t ≥ 0, where dW 1 t dW 2 t = ρ 1,2 dt and dW t dW i t = ρ i dt, i = 1, 2.

Some Special Cases
Recall that we suppose the past w i t , t ≤ 0 has been observed and the quantity g H i t is a deterministic function of time t ≥ 0 with g H i 0 = 0. However, it seems difficult to observe the whole past of w i t in the actual market, and so, we assume that g H i t = 0 as in [3] in the rest of this paper. Furthermore, following most of the previous research, we assume that the volatility risks are fully diversified, and so, the market prices of risksη i t are zero, i.e.,η i t = 0, for the sake of simplicity. Under these specifications, the fractional mean-reverting process X i t is given by: where λ H i 2 (t) is defined in (4). Given these processes X i t , the dynamics of the asset price S t is determined by the SDE (5).
We observe from (8) that the process X i t is a sum of a deterministic function h i (t) and a stochastic convolution I i (t), where: respectively. For the purpose of Monte Carlo simulation, we consider the discrete-time process x i n = X i (n∆t), n = 1, 2, . . . , N, where T = N∆t is the option maturity for sufficiently small ∆t > 0. The stochastic convolution I i (t) can be approximated by a discrete convolution such as ∑ n The discrete convolution can be represented by a matrix product, and so, the mean-reverting process x i n , n = 1, 2, . . . , N, in discrete time can be expressed compactly in matrix form as: . .
In the following, thanks to the above specification, we adopt matrix Equation (9) for Monte Carlo simulation (see, e.g., [9] for Monte Carlo simulation methods for general fBMs).

Approximation Formula
The call option price written on S t with strike K and maturity T is given by: where (x) + = max{x, 0}, and the expectation is taken over S T under the martingale measure Q. In this case, all we need to know is the distribution of asset price S T at the maturity T. Let us define X t = S t F(0,t) − 1. From (10), the value of the European call option is given by: With the density function f T (x) of X T at hand, it follows that: In this section, we derive an approximation formula for option prices by applying the technique used in [4] for the fractional volatility model (7). To this end, we expand the underlying asset into a sum of iterative stochastic integrals with deterministic integrands. In the following, we denote The proof of the next result is given in Appendix A.
Lemma 1. For X t := S t /F(0, t) − 1, we approximate it by: where: Following the idea given in [10], the density function f T (x) can be approximated byf T (x) given in (A16) in Appendix B. By computing (11), the next result then follows. The proof is straightforward, although messy, and omitted. Proposition 1. The value of a European call option with maturity T and strike K under the fractional volatility model (7) is approximated as: where n(x; µ, σ 2 ) is the density function of the normal distribution with mean µ and variance σ 2 and Φ(x) denotes the cumulative probability function of the standard normal distribution.
The definitions of h n (x), q k (T) and Σ T are provided in Appendix B.

Numerical Examples
This section is devoted to numerical studies of our fractional volatility model by using the approximation formula given in Proposition 1. (The accuracy of our approximation formula is checked by the Monte Carlo simulation explained in Section 2.2. We note that our approximation gets gradually worse for the high volatility, long maturity and deep in-the and out-of-the money cases. For example, when the percentage volatility defined by η = (γ 1 + γ 2 )/(X 1 0 + X 2 0 ) is bigger than 1.5, our approximation seems not enough for practical uses. For such cases, higher order approximation is required.) Throughout the numerical examples, we use the parameter values listed in Table 1 as the base case. In particular, for the Hurst indexes H 1 and H 2 , Bollerslev and Mikkelsen [1] observed long memory features in stock market volatilities, and so, we set H 1 = 0.9 as the persistent volatility. On the other hand, Bayer et al. [11] claimed (they report a good fit of their rough volatility model with H 2 = 0.07 and the percentage volatility η = 1.9) that H 2 is of order 0.1, and so, we set H 2 = 0.1 as the rough volatility. Finally, we assume that the volatility function is given by:  In Figure 1, we depict the ATM (at the money) implied volatility (Figure 1a) and the ATM skew (Figure 1b) with respect to the option maturity T. It is observed that the term structure of the ATM volatility skew is a power-law function of time to maturity T. This is a typical feature of rough volatility, which is observed in the S&P index options market reported by Bayer et al. [11]. The model ATM skew is approximated by the power-law function with the order of −0.449.  Table 1.
Next, based on our approximation formula, we investigate the effects of the model parameters on option prices.

Effect of H 2
In Figure 2, we plot the skew of the European options with respect to H 2 . It is explicitly observed that, as H 2 increases, the power-law index increases to be −0.449, −0.21, 0.044 and 0.353 for H 2 = 0.1, 0.3, 0.5 and 0.8, respectively. According to [5], an empirical study shows that the index is typically given by about −0.5, and so, our model is capable of capturing this stylized fact by setting the rough volatility H 2 close to zero.  Table 1.

Effect of H 1
We first check whether the Hurst index H 1 of persistent volatility has the ability to capture the power-law of the ATM skew or not. Figure 3 shows the ATM skew of the European options with respect to H 1 . In contrast to the rough volatility index H 2 given in Figure 2, the effect of H 1 on the ATM skew is very limited or almost has no effect. On the other hand, Figure 4 shows the volatility smile with respect to the strike and maturity.   Figure 4, thanks to the long memory feature of the index H 1 , the case of H 1 = 0.9 exhibits a slower decrease of the volatility smile amplitude with respect to time to maturity T than the short memory case H 1 = 0.5 (the rough volatility case H 1 = 0.4 shows a much faster decrease). This is a preferable feature, because the observed smile amplitude decreases much more slowly with respect to maturity than that explained by standard stochastic volatility models.  Table 1.
Summarizing, an important observation at this point is that, by setting H 1 > 0.5 and H 2 < 0.5, we can realize both the long memory feature and the roughness of volatility simultaneously by using our model, which is the major contribution of this paper.

Effect of ρ 1,2 and H 1
We examine the impact of ρ 1,2 on option prices. Recall that ρ 1,2 is the correlation between the persistent volatility X 1 t and the rough volatility X 2 t . Note that the ATM skew amplitude comes from the roughness of σ(X 1 t , X 2 t ). From Figure 5, it is observed that, as the correlation decreases, the decrease of the smile amplitude is decelerated in the case of H 1 ≤ 0.5. This means that the effect of the volatility persistence survives only when the correlation ρ 1,2 is negative.  Table 1.
As we can see form the figures, depending on the correlation values, the volatility surface is classified into four types: (1) the roughness exists, but the persistence does not; (2) the persistence exists, but the roughness does not; (3) both the roughness and the persistence exist; and (4) neither the roughness nor the persistence exist.
Specifically, as ρ 1 increases, it is observed from Figures 6a and 7a that the volatility smile is preserved and the skew becomes greater, i.e., the power-law index of the ATM skew decreases, respectively. Form Figures 6b and 7b, when the absolute value of ρ 2 is small, say |ρ 2 | < 0.2, it is observed that the volatility smile is preserved, but the ATM skew disappears, respectively. On the other hand, when the absolute value of ρ 2 is high, say |ρ 2 | > 0.7, we can see that the volatility smile gradually disappears, but the ATM skew becomes greater. Furthermore, from Figures 6c and 7c, as ρ 1,2 increases, the ATM skew becomes greater, but the volatility smile disappears, respectively.  Table 1.   Table 1.

Conclusions
In this study, we extend the fractional volatility model proposed in [4] by introducing another factor of rough volatility. Through numerical experiments, we demonstrate that, when one of the Hurst indexes in fractional volatility is larger than 1/2 (volatility persistence) and the other is smaller than 1/2 (rough volatility), our model can explain both the slower decay of the smile amplitude decline and the term structure of the at-the-money volatility skew observed in the options market, simultaneously. However, the coexistence of the two stylized facts seems non-trivial. Namely, depending on the three correlation values between the underlying asset and the factors of rough volatility and volatility persistence, the implied volatility surface is classified into four types: (1) the roughness exists, but the persistence does not; (2) the persistence exists, but the roughness does not; (3) both the roughness and the persistence exist; and (4) neither the roughness nor the persistence exist.
As future work, we plan to apply our model to actual markets. Namely, we develop a fast algorithm to calibrate our model to the options market, because our model involves many parameters itself. Furthermore, it is of interest to develop a model under the physical measure P to explain the volatility persistence and the volatility roughness simultaneously. An asymmetric model between the persistent volatility and the rough volatility may be of great importance. at n = 3. If the volatilities σ n (t) are deterministic functions and σ t = max n σ n t ∈ L 2 ([0, t]) is sufficiently small, then Proposition 2.2 of [12] assures that the sum of the iterated integrals converges very quickly.
Before proceeding, since σ 0 (s) is the deterministic function of s, we can apply the Wiener-Ito chaos expansion to S (1) t to obtain Hence, we define: as an approximation for S (1) t . Summarizing, we approximate the quantity S t by: where S (1) t is given by (A4) and I n (t) by (A3). In the following, we approximate each I n (t) by an iterated integral with deterministic volatilities.
In this paper, we employ Taylor's expansion around S (1) t for this purpose. Recall that J t (σ) = t 0 f (X 1 u , X 2 u )dW u . It follows that: and where we denote: for the sake of notational simplicity. We next approximate each I n (t) by using the approximations (A6) and (A7). Since I 1 (t) = J t ( f ) − J t (σ 0 ), from (A6) and (7), we get: Moreover, by applying Ito's formula, the second term on the left-hand side of (A8) is written as: Hence, we have: From the definition, we have: However, from (A7) and (7), we obtain: By substituting the conditional expectations into (A15), the approximate density function, denoted byf X t (x), can be expressed as: where h n (x) denotes the Hermite polynomial of order n: h n (x) = (−1) n e x 2 /2 d n dx n e −x 2 /2 , n = 1, 2, . . . , with h 0 (x) = 1.