Apple juice evaluation: Qualitative analysis and microsatellite traceability

1.   Introduction

Stochastic volatility models driven by fractional Brownian motion have been gaining attention in the area of mathematical finance, since many empirical studies find that the decay in the autocorrelation function of the volatility is better modeled by a power function than an exponential function and fractional stochastic volatility models generate better fits to the observed market implied volatility surface. Refer to, for example, Comte and Renault ^[9], Al$\mathrm{\grave{o}}$s et al. ^[2], Bayer et al. ^[4], Garnier and Solna ^[17], Forde and Zhang ^[12], Gatheral et al. ^[19], Guennoun et al. ^[20], Kim et al. ^[21], Bennedsen et al. ^[5], Fukasawa and Gatheral ^[16], Shi and Yu ^[30], and Cont and Das ^[10], among others. While Comte and Renault ^[9] stressed that fractional Brownian motion was relevant to capture long memory in stochastic volatility, Gatheral et al. ^[19] showed that log-volatility behaves essentially as a fractional Brownian motion with a Hurst exponent of order 0.1, at any reasonable time scale. A recent paper by Wang et al. ^[33] also found that the logarithmic daily realized volatility series of various financial assets have rough sample paths. However, fractional Brownian motion is not a semimartingale unless the Hurst exponent is $1/2$, as shown in Rogers ^[28], leading to a possible arbitrage opportunity (a free lunch with vanishing risk, as a general term). "Arbitrage" means profiting from a price gap between a derivative and a portfolio of assets that replicates the derivative's cash flows. So, to get around this problem while taking into account the long or short memory property, there have been studies using the mixed fractional Brownian motion of Cheridito ^[8], the approximate fractional Brownian motion of Thao ^[32], or the generalized fractional Brownian motion introduced by Pang and Taqqu ^[24].

We are interested in the approximate fractional Brownian motion introduced in Thao ^[32] in the context of modeling stochastic volatility for pricing derivatives. This process is a semimartingale even if the Hurst exponent is not $1/2$ (contrary to the original fractional Brownian motion) and so the no-arbitrage theory can be applied to obtain a partial differential equation (PDE) for the option price. There are several works which have used this process for underlying asset price models or stochastic volatility models. Refer to, for instance, Dung ^[11] for the Black-Scholes model, Sattayatham and Intarasit ^[29] for a jump-diffusion model, Pospisil and Sobotka ^[26] for the Heston stochastic volatility model, and Chang et al. ^[7] for the double Heston stochastic volatility model. In the present paper, we apply approximate fractional Brownian motions to the multiscale stochastic volatility model of Fouque et al. ^[14]. The merit of this stochastic volatility model is that the resulting option price approximation is independent of the particular details of the volatility model and leads to more flexibility in the parametrization of the implied volatility surface. We take two approximate fractional Brownian motions with two Hurst exponents (instead of two standard Brownian motions) corresponding to two characteristic (fast and slow) time scales of the multiscale volatility model, respectively. As far as we know, this approach does not exist in the literature. However, this type of volatility formulation is consistent with some of the previous works related to fractional Brownian motion. According to Xiao and Yu ^[35,36], the asymptotic distribution for the estimator of the persistence parameter is different when the Hurst exponent is less than 1/2 from that when it is larger than 1/2 in the fractional Vasicek model. Al$\mathrm{\grave{o}}$s and Leon ^[3] found, based on the Clark-Ocone-Haussman formula for the integrated variance, that the volatility can be composed of terms with a Hurst index less than $1/2$ being more relevant at short scales and terms with Hurst index greater than $1/2$ being more relevant at long scales. Also, Bennedsen et al. ^[5] discovered evidence consistent with the hypothesis that time series of realized volatility are both rough and very persistent. On the other hand, Cont and Das ^[10] observed interestingly that even when the instantaneous volatility has the same roughness as Brownian motion, the realized volatility exhibits behavior corresponding to a Hurst exponent significantly smaller than $1/2$. This observation supports our use of approximate fractional Brownian motion instead of fractional Brownian motion. The approximate fractional Brownian motion is thought of as between the standard Brownian motion and the fractional Brownian motion. It is a stochastic process equipped with a Hurst parameter, i.e., a measure of long-term memory of the time series. It is, however, a semimartingale in the form of a Brownian motion plus a time (Riemann) integral of an adapted process. So, the arbitrage opportunity can be excluded from the fundamental theorem of asset pricing and we are allowed to use the replicating portfolio method to obtain the corresponding PDEs for European vanilla and exotic options. The contribution of this work is as follows. We obtain approximate closed-form formulas for the prices of European vanilla and two exotic options. Our results cover both long- and short-memory properties of volatilities and control the skew slope by selecting appropriate Hurst exponents of the fast and slow scale volatility movements. It unifies two previously known results regarding the Hurst parameter dependence of the blow-up and slow-flattening behavior of skews and smiles of implied volatility surfaces. Consequently, the implied volatility surfaces can be calibrated over a wide range of time-to-maturities. Also, we provide a calibration method by representing the observed SPX option prices in terms of the term structure of the implied volatility formula. Based on the calibration result, we find that the implied volatility becomes higher when the fast-scale motion of the (spot) volatility becomes "rougher" and the slow-scale motion of the volatility becomes "smoother", which in turn supports the necessity of a multiscale modeling framework for stochastic volatility.

The paper is organized as follows. In Section 2, we use approximate fractional Brownian motions to establish a stochastic volatility model. In Section 3, we apply the replicating portfolio method to obtain the corresponding PDE formula for the price of a European vanilla option and derive explicitly a closed-form formula for the approximate option price using the combination of singular and regular perturbations and the Mellin transform method. Subsequently, a closed-form formula for the implied volatility corresponding to a European call option is obtained in Section 4. We check the accuracy of the pricing formula, show how to calibrate the pricing parameters, and investigate the sensitivity of the implied volatilities to the Hurst exponents in Section 5. We extend the vanilla option price formula to two exotic-option cases in Section 6. Finally, Section 7 provides some concluding remarks.

2.   Model formulation

A fractional Brownian motion $B^{H}_t$ with a Hurst exponent $H$, $0 < H < 1$, is defined by a centered Gaussian process satisfying the covariance function

The process $B^{H}_t$ is a self-similar process but is neither a semimartingale nor a Markov process except in the case where $H = 1/2$. Mandelbrot and van Ness ^[22] gave an integral representation of the general fractional Brownian motion as follows:

where $W_t$ is a standard Brownian motion. The last integral part of (2.1) is a self-similar Gaussian process, which becomes a Brownian motion for $H = 1/2$ and has non-stationary increments for $H\neq 1/2$. It is a truncated version of the general fractional Brownian motion and is usually called Riemann-Liouville fractional Brownian motion with a Hurst index $H$. This type of fractional Brownian motion has been widely used in the modeling of volatilities. See, for example, Comte and Renault ^[9], Al$\mathrm{\grave{o}}$s et al. ^[2], Bayer et al. ^[4], and Gatheral et al. ^[19], among others. It has a simple representation but it is not a semimartingale for $H\neq 1/2$. Thus Thao ^[32] used a perturbation parameter, say $\gamma$, to introduce an approximate fractional Brownian motion defined as

where $W_t$ is a standard Brownian motion, and proved that $B^{\gamma, H}_t$ is a semimartingale and converges to the last integral of (2.1) in $L^2(\Omega)$ sense and uniformly with respect to $t\in[0, T]$ for any fixed positive number $T$ when $\gamma$ goes to 0.

In this paper, we use the process $B^{\gamma, H}_t$ as a random source of the volatility of the underlying risky asset return to introduce a new model given by

under a risk-neutral probability measure $\mathbb{Q}$. The model (2.2) is the same as the multiscale stochastic volatility model in Fouque et al. ^[14] except that the standard Brownian motions driving the two stochastic volatility factors $Y_t$ and $Z_t$ are now replaced by the approximate fractional Brownian motions. Since, by the It$\hat{\mbox{o}}$ formula (see Oksendal ^[23], for example), the differential of the approximate fractional Brownian motion is

the initial model (2.2) becomes

where $\phi_{1, t}$ and $\phi_{2, t}$ are defined by

where $(t-s+\gamma)^{H-\frac{3}{2}}$ does not blow up at any $s\in [0, t]$ and $H\in(0, 1)$. In the model (2.3), we assume that $0 < \delta\ll \epsilon \ll \sqrt{\delta} < 1$, $W_t^x$, $W_t^y$, and $W_t^z$ are standard Brownian motions defined on a filtered probability space $(\Omega, \mathcal{F}, \mathcal{F}_t, \mathbb{Q})$, they have a correlation structure given by $d\langle W^x, W^y\rangle_t\! = \!\rho_{xy}dt$, $d\langle W^x, W^z\rangle_t\! = \!\rho_{xz}dt$, and $d\langle W^y, W^z\rangle_t\! = \!\rho_{yz}dt$, and $H_1\in(0, \frac{1}{2})$ and $H_2\in(\frac{1}{2}, 1)$. If $Y_t$ is a mean-reverting process, its mean itself is also mean-reverting slowly, and it is driven by $Z_t$, then this situation is called double-mean-reverting (cf. Gatheral ^[18]). $W^y$ and $W^z$ are assumed to be correlated here so that the model can somewhat capture the double-mean-reverting property of stochastic volatility even if the dependence of $Y_t$ on $Z_t$ is not explicitly specified. The functions $f$, $\alpha$, $\beta$, $g$, and $h$ are assumed to satisfy necessary smooth and boundedness conditions for the stochastic differential equation for $X_t$ to have a unique solution. The volatility factor $Y_t$ in (2.3) reflects rapid variation (for example, rapid mean reversion) while the volatility factor $Z_t$ represents slow variation because those processes correspond to the solutions of stochastic differential equations in which time $t$ is replaced by $t/\epsilon$ (sped up) and $\delta t$ (slowed down), respectively. Particularly, the process $Y_t$ is assumed to be ergodic and have an invariant distribution, denoted by $\Phi$, which allows us to use averaging principles in Fouque et al. ^[15] to approximate the option price. Of course, the model (2.3) is reduced to the model of Fouque et al. ^[14] if both $H_1$ and $H_2$ are equal to $1/2$.

While fractional Brownian motion was stressed to capture long memory in the early age of fractional stochastic volatility model development such as the study of Comte and Renault ^[9], it has been discovered empirically since then that those models are valid only for long-term behavior of volatility, while some rough volatility models are more appropriate in the short run (see, in particular, Gatheral et al. ^[19]). This has led several authors to introduce volatility models incorporating both roughness, meaning exponentially decaying autocorrelation, and long memory, meaning non-integrable autocorrelation, corresponding to two different Hurst exponents. Refer to Al$\mathrm{\grave{o}}$s and Leon ^[3] and Bennedsen et al. ^[5], for instance. This paper seeks to relate two characteristic (fast and slow) time scales of the multiscale stochastic volatility model of Fouque et al. ^[15] to the roughness and the long memory, respectively. This approach allows us to obtain an option pricing formula that can be calculated easily starting from the Black-Scholes price.

3.   Option price formula

3.1.   Singular perturbation problem

In this paper, the following lemmas are useful for asymptotic analysis. They are the solvability condition of a Poisson equation and the growth condition related to the infinitesimal generator of the ergodic process $Y_t$, $\frac{1}{\epsilon}\mathcal{A}_0$, where $\mathcal{A}_0$ is a differential operator defined as

Lemma 3.1. The Poisson equation

has a solution $p(t, x, y, z)$ if and only if the function $q$ is centered with respect to the invariant distribution $\Phi$ of the process $Y$, i.e.,

Proof. This is a version of the Fredholm alternative. Refer to Section 3.2 in Fouque et al. ^[15]. □

Lemma 3.2. Assume that equation $\mathcal{A}_0 p(t, x, y, z) = 0$ admits only solutions that do not grow as fast as

Then the solution $p$ does not depend on $y$.

Proof. Solving the equation $\mathcal{A}_0 p = 0$ directly leads to this result. □

Since the approximate fractional Brownian motion $B^{\gamma, H}_t$ is a semimartingale, the no-arbitrage theory is allowed for the model (2.3) by the fundamental theorem of asset pricing (see Pascucchi ^[25], for example). So, one can use the replicating portfolio approach to value the options. If $P^{\epsilon, \delta}(t, x, y, z; \phi_1, \phi_2)$ denotes the option price with a payoff function $H(x)$ under the model (2.3) when $X_t = x$, $Y_t = y$, $Z_t = z$, $\phi_{1, t} = \phi_1$, and $\phi_{2, t} = \phi_2$, then the no-arbitrage argument with the self-financing assumption and the It$\hat{\mbox{o}}$ formula leads to a final value problem expressed by

where the multiscale operator $\mathcal{A}^{\epsilon, \delta}$ is

where the operator symbol $\mathcal{D}_{n}$ is defined by

3.2.   Approximation

From now on, the dependence of $P^{\epsilon, \delta}$ on $\phi_1$ and $\phi_2$ is omitted for notational simplicity. Since the PDE problem (3.1) is a singular perturbation problem, we are interested in an asymptotic solution of the form

Following the multiscale asymptotic analysis of Fouque et al. ^[15] and using the operator $\mathcal{A}_{ \text{ BS}}$ defined by

one can find that $P_{00}$, $P_{01}$, and $P_{10}$ are independent of the variable $y$ and they satisfy the PDE problems

respectively, where the functions $A_{1}^{H_1}$, $A_{2}^{H_1}$, $B_{1}^{H_2}$, and $B_{2}^{H_2}$ are

respectively. Here, $\psi(y, z)$ is a function defined by the solution to

Note that $A_{1}^{H_1}(z)$ and $A_{2}^{H_1}(z)$ are related to the fast variation of volatility and a Hurst exponent less than $\frac{1}{2}$ while $B_{1}^{H_2}$ and $B_{2}^{H_2}$ are connected with the slow-scale variation of volatility and a Hurst exponent larger than $\frac{1}{2}$.

Since $\mathcal{A}_{ \text{ BS}}$ is the differential operator $\partial_t$ plus the infinitesimal generator of a geometric Brownian motion solving the stochastic differential equation

as its notation suggests, the PDE problem $\mathcal{A}_{ \text{ BS}} P_{00}(t, x, z) = 0$ with the final condition $P_{00}(T, x, z) = H(x)$ gives us that $P_{00}$ is the Black-Scholes option price (cf. Black and Scholes ^[6]) with constant volatility replaced by $z$-dependent volatility, and subsequently we use notation $P_{ \text{ BS}}(t, x, z)$ instead of $P_{00}(t, x, z)$ from now on.

By solving the PDE problems in (3.4) for $P_{ \text{ BS}}$, $P_{01}$, and $P_{10}$, we obtain the following European option price formula.

Proposition 3.1. Under the dynamics of (2.3) of the underlying asset price, the option price $P^{\epsilon, \delta}$ is approximated by $\ddot{P}^{\epsilon, \delta}\!: = \! P_{ \mathit{\text{BS}}}+\sqrt{\epsilon}P_{01}+\sqrt{\delta}P_{10}$, that is

where $\mathcal{D}_{1}^n\!: = \!(\mathcal{D}_{1})^n\! = \! \left(x\frac{\partial}{\partial x}\right)^n$ and $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$ are given by

respectively.

Proof. To solve the PDE problems in (3.4) for $P_{ \text{ BS}}$, $P_{01}$, and $P_{10}$, we use the Mellin transform and its inverse transform defined by

respectively, where $a$ is a real number, and obtain the ODE problems for $\hat{P}_{01}$ and $\hat{P}_{10}$ as follows.

where $\hat{h}(\omega)$ is the Mellin transform of $h(x)$ and the functions $\lambda(\omega, z)$, $\eta_{1}(\omega, z)$, and $\eta_{2}(\omega, z)$ are

respectively. The solutions of (3.8) are given by

Substituting (3.9) into (3.10), we obtain the following Mellin transform of $\ddot{P}^{\epsilon, \delta}$ explicitly:

The pricing formula (3.11) is given by a linear combination of terms that are in the form of a product of $\hat{P}_{ \text{ BS}}$ and a polynomial function of $\omega$. So, we can calculate the formula explicitly through the following property of the Mellin transform: $\mathcal{M}((\mathcal{D}_{1})^{n} f)(\omega) = \, (-\omega)^{n}\hat{f}(\omega)$. Using an inverse Mellin transform on (3.11), we obtain a closed-form formula given in the proposition. □

Therefore, once the Black-Scholes option price $P_{ \text{ BS}}$ is given, we can calculate the approximation $\ddot{P}^{\epsilon, \delta}$ by computing the derivatives of $P_{ \text{ BS}}$ in the formula (3.7) and plugging the estimated group parameters $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$ into (3.7). Note that the original model parameters and functions such as $\gamma$, $\epsilon$, $\delta$, $\alpha$, $\beta$, $g$, and $h$ are not required to be directly chosen for the purpose of the option price approximation. The group parameters are required to be chosen for calculating $\ddot{P}^{\epsilon, \delta}$. We use the implied volatility term structure of SPX call options to estimate those four group parameters. See Section 5.2 for details.

Remark: One can obtain the second-order terms $P_{02}$, $P_{11}$, and $P_{20}$ further and a subsequent formula for the approximation $\dddot{P}^{\epsilon, \delta}: = P_{ \text{ BS}}+\sqrt{\delta}P_{10}+\sqrt{\epsilon}P_{01}+\delta P_{02}+\sqrt{\delta\epsilon}P_{11}+\epsilon P_{20}$ as shown in Appendix A, where the approximation error is also given. We note that the same approach (approximation to the Black-Scholes price) used in this section and the following section are also used in the more simple framework of the Heston model driven by the standard Brownian motion in Al$\mathrm{\grave{o}}$s et al. ^[1].

4.   Implied volatility

All of the original model parameters given in the model (2.3) are not required to price derivatives. In fact, from Proposition 3.1, we notice that $\bar{\sigma}_{f}$, $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$ are the ones required to be estimated. In order to estimate those pricing parameters, we can utilize calibration from near-the-money European call option implied volatilities. The volatility $I^{\epsilon, \delta}$ implied by the pricing formula (3.7) is defined by the solution to the equation $P_{ \text{ BS}}(t, x; I^{\epsilon, \delta}) = \ddot{P}^{\epsilon, \delta}(t, x, z),$ where $P_{ \text{ BS}}(t, x; \sigma)$ stands for the classical Black-Scholes call option price formula with volatility $\sigma$. Then the two correction terms $I_{01}$ and $I_{10}$ of the asymptotic expansion

are given by $I_{01}\! = \! \left(\partial_{\sigma} P_{ \text{BS}} \right)^{-1}P_{01}$ and $I_{10}\! = \! \left(\partial_{\sigma} P_{ \text{BS}} \right)^{-1}P_{10}$, respectively, and the leading term $I_{00}$ is defined as

Using the vega-gamma and speed-gamma relationships, i.e.,

we can obtain an approximate implied volatility surface $\ddot{I}^{\epsilon, \delta}\!: = \! I_{00}+\sqrt{\epsilon}I_{01}+\sqrt{\delta}I_{10}$ given by

We are interested in the slope behavior of the implied volatility skew with respect to time-to-maturity and the Hurst exponent. The slope of the implied volatility skew is given by

where $k: = \log K$.

Based on the calibration result shown in Section 5.2 for SPX options, we do a numerical experiment to show how the skew slope behaves against time-to-maturity. Figure 1 presents the experimental result for the skew slope term structure. It shows that the slope tends to blow up as time-to-maturity becomes shorter.

Comte and Renault ^[9] considered a Hurst exponent with $H > \frac{1}{2}$ to explain the slow flattening of skews and smiles of the implied volatility surface when time-to-maturity increases, while Alòs et al. ^[2] gave a better description of the short time-to-maturity blow-up of the implied volatility surface with a Hurst exponent $H < \frac{1}{2}$. So, our result is consistent with these two results in that the skew slope becomes small when time-to-maturity increases while it becomes large when time-to-maturity decreases as seen in Figure 1. Our result unifies these two separate previous results. This is desirable in practice as the market volatility has both long- and short-memory properties depending on the situation. As described by (4.2), the skew slope $\frac{\partial \ddot{I}^{\epsilon, \delta}}{\partial k}$ can behave flexibly depending on the appropriate Hurt parameters $H_1$ and $H_2$ in the range of $(0, \frac{1}{2})$ and $(\frac{1}{2}, 1)$, respectively. As a consequence, the implied volatility surface can be calibrated over a wide range of time-to-maturities.

5.   Performance, calibration, sensitivity, and accuracy

In this section, we check the accuracy and performance of the price formula given by (3.7) in Proposition 3.1 for European call options via Monte Carlo simulation. We give an example of calibrating the pricing parameters for three different time-to-maturities. We also investigate the sensitivity of the implied volatility to the Hurst exponents $H_1$ and $H_2$.

Using the well-known Greeks in the Black-Scholes model, one can verify easily that the derivative price $\ddot{P}(t, x, z)$ given by (3.7) satisfies the following identity for a European call option with strike $K$:

where $\mathcal{N}$, $d_1$, $d_2$, and $\varphi$ are given by

respectively.

5.1.   Performance of the formula

In this section, we calculate option prices under three different models and a Monte Carlo simulation result (with 1 million simulations) and compare them with real market data, that is, the SPX option data observed on December 7, 2022. Notations $P_{\text{market}}$, $P_{\text{MC}}$, $P_{\text{fMSV}}$, $P_{\text{MSV}}$, and $P_{\text{BS}}$ stand for the market option price, the Monte Carlo (MC) simulation result, the option price computed by the formula (5.1) corresponding to the fractional multiscale stochastic volatility (fMSV) model (2.3), the option price computed under a multiscale stochastic volatility (MSV) model corresponding to the case of $H_{1} = H_{2} = \frac{1}{2}$ in the formula (5.1), and the Black-Scholes option price, respectively.

For Monte Carlo simulation, we need random numbers generated by the stochastic processes $\phi_{1, t}$ and $\phi_{2, t}$ in (2.4). The time interval $[0, T]$ is discretized into $t_{0}( = 0)$, $t_1$, $t_2$, $\cdots$, $t_{n}( = T)$ satisfying $t_0 < t_1 < \cdots < t_n$ with $\Delta t = t_{i}-t_{i-1}$, $i = 1, 2, \cdots, n$. Omitting the subscript number and superscript letter of $\phi_{1, t}$, $\phi_{2, t}$, $H_1$, $H_2$, $W^{y}_s$, and $W^{z}_s$ (and so using $\phi_{t}$, $H$, and $W_s$), the random source $\phi_{t}$ satisfies the following recursive equation by the It$\hat{\rm{o}}$ calculus:

where the notation $\stackrel{d}{ = }$ denotes distributional equality and the $Z_i$'s are independent and identically (standard normal) distributed. Thus, we generate the random sources $\phi_{t_i}$, $i = 1, 2, \cdots, n$, recursively, by the following algorithm:

Applying this algorithm to (2.3), we can obtain a Monte Carlo simulation result $P_{\text{MC}}$ for the European options. Table 1 represents the setting of the related parameters and functions for finding the price $P_{\text{MC}}$.

We calculate $P_{\text{market}}$, $P_{\text{MC}}$, $P_{\text{fMSV}}$, $P_{\text{MSV}}$, and $P_{\text{BS}}$ for a European call option with two time-to-maturities, where $P_{\text{MC}}$ is obtained based on the setting in Table 1. Using these results, we compute the square norms $\parallel P_{\text{market}}-P_{\text{MC}} \parallel$, $\parallel P_{\text{market}}-P_{\text{fMSV}} \parallel$, $\parallel P_{\text{market}}-P_{\text{MSV}} \parallel$, and $\parallel P_{\text{market}}-P_{\text{BS}} \parallel$ for the purpose of comparing the Monte Carlo simulation and the three different volatility models. Table 2 presents the result. It shows that the fMSV model outperforms the other models including the Monte Carlo simulation. This tends to be the case more conspicuously when time-to-maturity becomes shorter. Note that shorter time-to-maturity options tend to have higher trading volume in general. Therefore, using the approximate fractional Brownian motion instead of the standard Brownian motion for volatility seems to provide a great advantage in option pricing. On the other hand, Table 3 provides the elapsed time of option pricing based on the Monte Carlo simulation (reputation number = $1,000,000$ and $dt$ = $1/25200$), and the three different volatility models. The computer used for the computation is specified as Intel(R) Core(TM)-i9-10900 CPU, Windows 10 Pro O/S and 64GB RAM. Moreover, the program used for the computation is MATLAB R2022b. The Monte Carlo simulation method takes much more time than the analytic methods based on the three different volatility models while the three different models are relatively similar to each other in terms of the elapsed pricing time.

Among those parameters required to be estimated, i.e., $\bar{\sigma}_{f}$, $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$, the parameter $\bar{\sigma}_{f}$ is first estimated from historical SPX data over a period of time in the near past, where the slow-scale variable ($z$) dependence of $\bar{\sigma}_{f}$ accounts for updating the long-run average from time to time. To estimate the group parameters $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$, we rewrite the implied volatility surface (4.1) as

where the parameters $a^{\epsilon}$, $b^{\epsilon}$, $c^{\delta}$, and $d^{\delta}$ are related to the pricing parameters $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$ through the relationship

So, once $a^{\epsilon}$, $b^{\epsilon}$, $c^{\delta}$, and $d^{\delta}$ are estimated from calibration to the implied volatility term structure of SPX call options, one can use the relationship (5.3) to estimate the pricing group parameters $A^{\epsilon, H_1}_1$, $A^{\epsilon, H_1}_2$, $B^{\delta, H_2}_1$, and $B^{\delta, H_2}_2$ and calculate the derivative price $\ddot{P}^{\epsilon, \delta}(t, x, z)$ obtained in Proposition 3.1.

More concretely, the averaged volatility $\bar{\sigma}_{f}$ is first estimated using the 10-day historical volatilities calculated from the SPX data obtained from the site http://www.investing.com and then the parameters $a^{\epsilon}$, $b^{\epsilon}$, $c^{\delta}$, and $d^{\delta}$ are estimated using the SPX call option data obtained from the site http://www.barchart.com. Figure 2 shows the implied volatilities of the SPX option in the real market and the curve $\ddot{I}^{\epsilon, \delta}$ fitted to the market data. From this fit, $a^{\epsilon}$, $b^{\epsilon}$, $c^{\delta}$, and $d^{\delta}$ are estimated and then the pricing group parameters $A_{1}^{\epsilon, H_1}$, $A_{2}^{\epsilon, H_1}$, $B_{1}^{\delta, H_2}$, and $B_{2}^{\delta, H_2}$ are determined. Figure 3 demonstrates the corresponding result observed at one day.

5.2.   Sensitivity to $H_1$ and $H_2$

The Hurst exponents $H_1$ and $H_2$ related to the fast- and slow-scale variations of the volatility, respectively, are important parameters in our underlying asset price model (2.3). To investigate the dependence of the implied volatility (4.1) on the Hurst exponents, we rewrite (4.1) as follows:

where $a_1^{\epsilon}(H_1)$, $a_2^{\epsilon}(H_1)$, $b_1^{\delta}(H_2)$, and $b_2^{\delta}(H_2)$ are

respectively.

In Figure 4, we demonstrate the behavior of the implied volatility $\ddot{I}$ with respect to the Hurst parameters $H_1$ and $H_2$. Figure 4(a) shows $\ddot{I}$ against $H_1$ for a variety of $H_2$ while Figure 4(b) shows $\ddot{I}$ against $H_2$ for a variety of $H_1$. The figures indicate that $\ddot{I}$ increases as $H_1$ decreases or $H_2$ increases for any fixed $H_2$ or $H_1$, respectively. So, the implied volatility becomes higher when the fast-scale motion of the volatility becomes "rougher" and the slow-scale motion of the volatility becomes "smoother". This is an interesting result in view of modeling stochastic volatility. This provides us with one of the reasons why we need a multiscale framework for stochastic volatility.

5.3.   Accuracy

In this section, we check the accuracy of the approximation $P_{\text{fMSV}}$ with respect to the parameters $\gamma$, $\epsilon$, and $\delta$. We choose a very small value, close to zero, of each of the parameters and show that the approximations converge to the Monte Carlo simulation result, denoted by $P_{\text{MC*}}$, with the small parameter(s).

Table 4 shows how $P_{\text{fMSV}}$ converges to $P_{\text{MC*}} = 41.2643$ as $\gamma$ goes to $0.0001$. On the other hand, Table 5 shows how the approximations $P_{\text{fMSV}}$ move to the Monte Carlo simulation result $P_{\text{MC*}} = 41.2576$ when $\epsilon$ and $\delta$ go to $1.0000 \times 10^{-6}$ and $1.0000 \times 10^{-9}$, respectively. We note that the correlation term $\rho_{yz}$ does not appear in the first-order approximation of our interest in this article but it would appear in higher order approximation. We assume $\rho_{yz} = 0$ in the numerical experiment in order to match the situation given by our first-order approximation.

6.   Extension to exotic options

In general, the prices of exotic options can be determined after option pricing models are calibrated to market data of plain vanilla options. In this section, we extend the pricing result for European vanilla options under the fractional multiscale stochastic volatility model (2.3) to two types of path-dependent exotic options, i.e., barrier and lookback options.

6.1.   Barrier option

Barrier options are similar to vanilla options but they only become activated or extinguished when the underlying asset hits a specific price level (the so-called "barrier"). So, the value of barrier options can jump up or down greatly. This type of option is commonly traded in the foreign exchange and equity markets.

Given the model (2.3), let $U^{\epsilon, \delta}(t, x, y, z)$ be the price of a down-and-out (D/O) barrier option, where a payoff function is given by

with a strike price $K$, a barrier level $B$, and an expiration time $T$. From the no-arbitrage theory with the self-financing condition and the It$\hat{\mbox{o}}$ formula, $U^{\epsilon, \delta}(t, x, y, z)$ satisfies the PDE problem

where the multiscale operator $\mathcal{A}^{\epsilon, \delta}$ is defined by (3.2).

We are going to derive a solution of the form

as an approximate solution of the PDE problem (6.2). Substituting the series expansion (6.3) into the PDE problem (6.2) and using the same methodology as used for Proposition 3.1, one can have the following PDE problems for the terms $U_{ij}(t, x, y, z)$, $(i, j) \in \left\{(0, 0), (0, 1), (1, 0) \right\}$:

where $A_{1}^{H_1}(z)$, $A_{2}^{H_1}(z)$, $B_{1}^{H_2}(z)$, and $B_{2}^{H_2}(z)$ are given by (3.5).

The following lemma is useful to solve the PDE problems for the terms $U_{01}$ and $U_{10}$.

Lemma 6.1. Consider the PDE problems

If $\xi$ satisfies the equation $\mathcal{A}_{ \mathit{\text{BS}}} \xi(t, x, z) = 0$, then the solution $u(t, x, z)$ can be decomposed into

where $u_{1}$ and $u_{2}$ are solutions to the PDE problems given by

respectively. Furthermore, the solutions $u_{1}$ and $u_{2}$ are given by

respectively.

Proof. This lemma is related to the Black-Scholes framework with volatility $\bar{\sigma}_{f}$. Refer to Section 6.2 in Fouque et al. ^[15] for a proof. □

Based on Lemma 6.1, we obtain the following semi-closed form formula for an approximate value of $U^{\epsilon, \delta}(t, x, y, z)$.

Proposition 6.1. Under the dynamics of (2.3), the option price $U^{\epsilon, \delta}(t, x, y, z)$ is approximated by $\ddot{U}^{\epsilon, \delta}\!: = \! U_{ \mathit{\text{BS}}}+\sqrt{\epsilon}U_{01}+\sqrt{\delta}U_{10}$, that is

where $U_{ \mathit{\text{BS}}}(t, x, z)$ is defined by

and $\mathcal{H}^{\epsilon, H_1}_{ \mathit{\text{01}}}$ and $\mathcal{H}^{\delta, H_2}_{ \mathit{\text{10}}}$ are differential operators defined by

respectively. Here, $P_{ \mathit{\text{BS}}}\left(t, x, z; \varpi\right)$ is the Black-Scholes call option price given by

where $d_{i}(s, z; \varpi)$ $(i = 1, 2)$ are defined in (5.1). The group parameters $A^{\epsilon, H_1}_1(z)$, $A^{\epsilon, H_1}_2(z)$, $B^{\delta, H_2}_1(z)$, and $B^{\delta, H_2}_2(z)$ are defined in Proposition 3.1.

Proof. First of all, since $U_{00}$ satisfies a PDE problem for a D/O barrier option under the Black-Scholes model with volatility $\bar{\sigma}_{f}$ as seen in (6.4), it becomes $U_{ \text{ BS}}$ defined in the proposition.

On the other hand, from (6.4), we have the following PDE problems for $\ddot{U}_{01}: = \sqrt{\epsilon}U_{01}$ and $\ddot{U}_{10}: = \sqrt{\delta}U_{10}$:

respectively. Applying Lemma 6.1 to (6.7) directly, one can obtain solutions for $\ddot{U}_{01}$ and $\ddot{U}_{10}$ which lead to $\ddot{U}^{\epsilon, \delta}(t, x, z)$ given in the proposition when they are added to $U_{00}$. □

6.2.   Lookback option

A lookback option is an exotic option that allows the holder to exercise an option at the most favorable (minimum or maximum) price of the underlying asset over the life of the option. The floating strike lookback option eliminates the risk associated with the market entry time. In this section, we obtain a pricing formula for the floating strike lookback call option under the model (2.3). In terms of a stochastic process defined by

(the minimum value from the contract time 0 until the current time $t$), we let $V^{\epsilon, \delta}(t, m, x, y, z)$ denote the price of the lookback call option, where the payoff function $H(x, m)$ is given by

From the no-arbitrage theory with the self-financing condition and the It$\hat{\mbox{o}}$ formula, $V^{\epsilon, \delta}(t, m, x, y, z)$ satisfies the PDE problem

We solve the PDE problem (6.9) for $V^{\epsilon, \delta}(t, x, m, y, z)$ using the asymptotic series expansion

Plugging (6.10) into the PDE in (6.9) and applying the same methodology as used for Proposition 3.1 yield that $V_{ij} (t, x, m, y, z)$, $(i, j)\in \left\{(0, 0), (0, 1), (1, 0) \right\}$, satisfy the PDE problems

respectively.

The following lemma is useful to solve the PDE problems for the terms $V_{01}$ and $V_{10}$.

Lemma 6.2. Consider a PDE problem given by

Then the solution $v(t, x, z)$ is given by

where $G(\tau, w, x, z)$ is given by

Here, the function $\varphi$ is defined in (5.2).

Proof. As the PDE in (6.12) is a non-homogeneous linear PDE with a Neumann boundary condition, by the change of variables, it can become a non-homogeneous heat equation whose solution is well-known. Refer to Polyanin and Nazaiknskii ^[27] for details. □

Based on Lemma 6.2, we obtain the following formula for an approximate value of $V^{\epsilon, \delta}(t, x, m, z)$.

Proposition 6.2. Under the dynamics of (2.3), the option price $V^{\epsilon, \delta}(t, x, m, z)$ is approximated by $\ddot{V}^{\epsilon, \delta}\!: = \! V_{ \mathit{\text{BS}}}+\sqrt{\epsilon}V_{01}+\sqrt{\delta}V_{10}$, that is

where $V_{ \mathit{\text{BS}}}(t, x, m, z)$ is defined as

$\mathcal{H}^{\epsilon, H_1}_{ \mathit{\text{01}}}$ and $\mathcal{H}^{\delta, H_2}_{ \mathit{\text{10}}}$ are the differential operators defined by (6.5), and $G$ is the function defined in Lemma 6.2. Here, $P_{ \mathit{\text{BS}}}(t, x, z; \varpi)$, $\varphi\left(\varpi, \mu, \sigma\right)$, and $d_1(x, z; \varpi)$ are defined in (6.6) and (5.2), respectively.

Proof. Above all, as the solution $V_{00}$ of the first PDE problem in (6.11) is exactly the price of the floating strike lookback call option under the Black-Scholes model whose volatility is $\bar{\sigma}_f(z)$, it is the same as $V_{ \text{ BS}}$ given by (6.14) which can be found in, for instance, Wilmott ^[34].

To derive approximate closed-form formulas for $V^{\epsilon, \delta}$, we let $\ddot{V}_{01}: = \sqrt{\epsilon}V_{01}$ and $\ddot{V}_{10}: = \sqrt{\delta}V_{10}$, and apply the reduction method of the dimension used in Shreve ^[31] to the PDE problems in (6.11). Letting $w: = \frac{x}{m}$ and $\ddot{V}_{ij}(t, x, m, z) = m\ddot{V}_{ij}\left(t, \frac{x}{m}, 1, z\right) = :mW_{ij}(t, w, z)$, $(i, j) = \{(0, 1), (1, 0)\}$, (6.11) becomes

where

Applying Lemma 6.2 to the PDE problems (6.15) for $W_{10}$ and $W_{01}$ directly, one can obtain the solutions corresponding to the formula (6.13). □

7.   Conclusions

In this paper, we have introduced a semimartingale approximation of fractional stochastic volatility in terms of two approximate fractional Brownian motions corresponding to two characteristic time scales. Based on the semimartingale property, we make use of the replicating portfolio method to obtain the parabolic PDE problems for European vanilla, barrier, and lookback options, and then solve those problems explicitly and derive approximate closed-form formulas for the option prices. The mixture of the Hurst parameters and the multiple time scales of volatility can unify effectively the previously known separate results about the time-to-maturity dependence of the blow-up or flattening behavior of the skews of the implied volatility. So, knowing that stochastic volatility models driven by fractional Brownian motions can generate better fits to implied volatility surfaces, our uniform approximation result can contribute to the situation that the volatility parameters including the Hurst exponent should be calibrated over a wide range of time-to-maturities.

Author contributions

Min-Ku Lee: Resources, software, data curation, methodology, validation, visualization, and writing-original draft; Jeong-Hoon Kim: Conceptualization, formal analysis, writing-original draft, writing-review & editing, supervision, project administration, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Data and code availability

Datasets analyzed during the current study are available through the sites http://www.investing.com and http://www.barchart.com. Code for parameter calibration is provided as part of the replication package and is available at https://sites.google.com/view/two-scale-sv-fbms.

Acknowledgments

The authors deeply thank the four anonymous reviewers for their valuable comments and suggestions to improve this paper.

The author J.-H. Kim gratefully acknowledges the financial support from the National Research Foundation of Korea grant NRF2021R1A2C1004080.

Conflict of interest

The authors certify that there is no actual or potential conflict of interest in relation to this article.

Appendix

A.   Second-order approximation

In this section, we obtain the second-order terms $P_{02}$, $P_{11}$, and $P_{20}$ in the asymptotic series expansion (3.3). The corresponding PDE problems have the terminal conditions $\langle P_{02}(T, x, \cdot, z) \rangle = 0$, $P_{11}(T, x, y, z) = 0$, and $P_{20}(T, x, y, z) = 0$, respectively. Here, the averaged terminal condition $\langle P_{02}(T, x, \cdot, z) \rangle = 0$ and the condition

are imposed based on the terminal layer analysis given by Fouque et al. ^[13].

We first obtain the following lemma for the terms $P_{02}$, $P_{12}$, and $P_{03}$ in (3.3).

Lemma A.1. The second-order terms $P_{02}$, $P_{12}$, and $P_{03}$ in (3.3) can be expressed as

for some functions $F_{02}$, $F_{12}$, and $F_{03}$ independent of the variable $y$, where $\eta$, $\xi$, and $\zeta$ are the solutions of

respectively.

Proof. Putting (3.3) into (3.1), Lemmas 3.1 and 3.2 draw forth Poisson equations given by

Using the solutions $\eta$, $\xi$, and $\zeta$ of (A.3), we can derive the solutions $P_{02}$, $P_{12}$, and $P_{30}$ in the form of (A.2) for some $y$-independent functions $F_{02}$, $F_{12}$, and $F_{03}$. □

Next, we obtain PDEs for $F_{02}$ (and so $P_{02}$), $P_{11}$, and $P_{20}$.

Proposition A.1. The second-order solutions $F_{02}$ in (A.2), $P_{11}$, and $P_{20}$ are independent on the variable $y$ and satisfy the following PDEs

with the boundary conditions $F_{02}(T, x, z) = P_{11}(T, x, z) = P_{20}(T, x, z) = 0$.

Proof. First of all, similarly with the proof of Lemma A.1, we can obtain the PDEs given by

which yields that $P_{11}$, $P_{20}$, and $P_{21}$ are independent on the variable $y$. Additionally, we can get the PDE given by

Applying Lemma 3.1 and the $y$-independence of $P_{11}$, $P_{20}$, and $P_{21}$ into (A.6), we can have the PDEs

and

Thus, we can obtain the PDEs for $P_{11}$ and $P_{02}$ in (A.5).

On the other hand, putting $P_{02}$ in (A.2) into (A.7), we can obtain

where the assumption (A.1) has been used. Thus the first equation for $F_{02}$ in (A.5) has been derived. □

Finally, the second-order approximation $\dddot{P}(t, x, y, z): = P_{ \text{ BS}}+\sqrt{\delta}P_{10}+\sqrt{\epsilon}P_{01}+\delta P_{02}+\sqrt{\delta\epsilon}P_{11}+\epsilon P_{20}$ and its accuracy are obtained as follows by solving the PDEs (A.5) in Proposition A.1.

Proposition A.2. Under the dynamics of (2.3) of the underlying asset price, the option price $P^{\epsilon, \delta}$ is approximated by $\dddot{P}(t, x, y, z): = P_{ \mathit{\text{BS}}}+\sqrt{\delta}P_{10}+\sqrt{\epsilon}P_{01}+\delta P_{02}+\sqrt{\delta\epsilon}P_{11}+\epsilon P_{20}$, where the first approximation part $\ddot{P}(t, x, y, z): = P_{ \mathit{\text{BS}}}+\sqrt{\delta}P_{10}+\sqrt{\epsilon}P_{01}$ is given by Proposition 3.1 and $P_{02}$, $P_{11}$, and $P_{20}$ are given by

where $A^{k}_{ij}$, $B^{k}_{ij}$, $C^{k}_{ij}$, and $D^{k}_{ij}$ are set aside in Appendix B for comfortable readability. Moreover, the approximation has the accuracy

for any $l < 1$.

Proof. First of all, we can rewrite the solutions $P_{01}$ and $P_{10}$ in Proposition 3.1 as

where $\mathcal{B}$ and $\mathcal{C}$ are the operators defined by

respectively. Putting the solutions (A.10) for $P_{01}$ and $P_{10}$ and (A.2) for $P_{02}$, $P_{12}$, and $P_{03}$ into (A.5), we have the following PDEs

where $\mathcal{G}$, $\mathcal{H}$, $\mathcal{M}$, and $\mathcal{N}$ are the differential operators defined by

Here, $D^{H_1, H_2}_{i}$, $D^{H_1}_{i}$, $D^{H_2}_{i}$, and $E^{H_2}_{i}$ are defined in Appendix B. By applying the following results (motivational equations) in Fouque et al. ^[13] to the PDEs (A.11), we solve the PDEs (A.11) to obtain the solutions (A.9).

The proof of the error estimate should be similar to the proof in Fouque et al. ^[13] and so we omit the proof. □

B.   The constants in the solutions $F_{02}$, $P_{11}$, and $P_{20}$