SPX Calibration of Option Approximations under Rough Heston Model

: The volatility of stock return does not follow the classical Brownian motion, but instead it follows a form that is closely related to fractional Brownian motion. Taking advantage of this information, the rough version of classical Heston model also known as rough Heston model has been derived as the macroscopic level of microscopic Hawkes process where it acts as a high-frequency price process. Unlike the pricing of options under the classical Heston model, it is signiﬁcantly harder to price options under rough Heston model due to the large computational cost needed. Previously, some studies have proposed a few approximation methods to speed up the option computation. In this study, we calibrate ﬁve different approximation methods for pricing options under rough Heston model to SPX options, namely a third-order Padé approximant, three variants of fourth-order Padé approximant, and an approximation formula made from decomposing the option price. The main purpose of this study is to ﬁll in the gap on lack of numerical study on real market options. The numerical experiment includes calibration of the mentioned methods to SPX options before and after the Lehman Brothers collapse.


Introduction
Stock price has been modelled by many stochastic processes/models, e.g., geometric Brownian motion (Black-Scholes model) [1], Merton jump diffusion model [2], Hull-White model [3], classical Heston model [4], and many more models, but they do not reflect the high-frequency movement of the stocks. In particular, most of the stochastic volatility models failed in modelling short-term implied volatility behaviour or more specifically, the term structure at-the-money (ATM) implied volatility skew. ATM implied volatility skew refers to the absolute rate of change of the implied volatility with respect to the log-moneyness k = log(K/S) with k = 0. The empirical implied volatility generates a skew of somewhere near O(T −1/2 ) when T → 0, whereas the above-mentioned models all have ATM implied volatility skew of O(T) when T → 0.
A study done by Fukasawa [5] has made a large impact in the finance field of volatility modelling. In particular, the study indicates that stochastic volatility model where the volatility process is driven by fractional Brownian motion (fBm) with Hurst parameter H can generate ATM implied volatility skew of O(T H−1/2 ) when T → 0. This is significant and different from the local-stochastic volatility model [6] as it only possesses an extra Hurst parameter H in addition to the normal stochastic volatility model, whereas the local-stochastic volatility model would normally contain many parameters, which is not very ideal for a model. In other words, the stochastic volatility models driven by fBm is a parsimonious model.
The fact that the volatility is empirically tested to be rough is largely due to Gatheral and his co-authors [7]. Their main result of the study indicates that the log-volatility follows a generalization of Brownian motion also known as fractional Brownian motion [8] with Hurst parameter 0 < H < 1/2. A model that contains fractional Brownian motion's features (explosive behaviour) with Hurst parameter 0 < H < 1/2 is now known as "rough volatility model". The name is given due to its anti-persistency trajectory.
Currently, there are couple of rough volatility models, namely the rough Bergomi model [9], rough Heston model [10,11], and rough SABR model [12]. In this study, we will focus on rough Heston model, which is also a rough version of Heston model [4]. It is noted that the rough Heston model is actually a microscopic Hawkes price model on a macroscopic scale [10]. In particular, the model aims to replicate several high-frequency behaviour of the trading activities, e.g., the absence of statistical arbitrage opportunities in trading, no real economic motivation for most of the trades (market is highly endogenous), the presence of metaorders (trading algorithms are splitting large transactions of trading orders to smaller-sized orders), and liquidity asymmetry in bid-ask orders (a sell order might decrease the stock price more than an increase for a purchase order).
The characteristic function of rough Heston model is derived by El Euch and Rosenbaum [11]. The main component of the characteristic function involves a fractional Riccati equation under the Riemann-Liouville fractional calculus. As of now, there is no closedform solution for the fractional Riccati equation, therefore, we must rely on numerical experiment to solve for its solution. Fractional Adams-Bashforth-Moulton is first studied by Diethelm et al. [13], it involves using the quadrature theory. The fractional Adams method is similar to the construction of product rectangle and trapezoidal method, but it has an iterative procedure to solve the fractional Riccati equation's solution (presence of Volterra integral equation). In addition, fractional Adams method has computational complexity of O(N 2 ) with N as the number of time steps.
To price an option under the rough Heston model, we would need to use Fourier inversion method through Carr and Madan's method [14] or Lewis' method [6,15]. Therefore, the computational complexity to price an option would be O(n L N 2 ) (if Lewis' method is used), where n L is number of space steps in Lewis' integral. Before we discuss the main problem and objective of this study, we would like to mention some works related to ours. There are numerous studies that approximate fractional Riccati equation's solution through infinite series expansions, e.g., Adomian decomposition method [16], variational iteration method [17], homotopy perturbation method [18], and homotopy analysis method [19].
The main aim of this study is to further test the option price computed through thirdorder Padé approximant (GPadé) [20,21] and fourth-order Padé approximant (FPadé) [22] as well as the approximation formula for option price under rough Heston model [23] on real option data as the methods are lacking real option data to verify their effectiveness in practice. The Padé approximation method is used as a rational approximation of fractional Riccati equation's solution, then we compute the option price through Lewis' method, whereas the approximation formula [22] is a truncated decomposition formula of rough Heston's option price. The numerical experiment of S&P 500 (SPX) options will be carried out on two separate scenarios, the first scenario is on 18 June 2021 and the other scenario is before and after the Lehman collapse crisis (12 and 15 September 2008). The paper is organized as follows. Section 2 provides an introduction of rough Heston model and its characteristic function. In Section 3, we study the Padé methods [20,22] and the approximation formula [23]. In Section 4, the numerical experiment of the mentioned methods on SPX options will be discussed.

Rough Heston Model and Its Characteristic Function
Rough Heston model is a form of stochastic Volterra equation, in other words, the rough Heston model is not a Markovian model as the variance process relies on its past self. This also indicates that if the Monte Carlo method is available, one simulation of variance process requires O(N 2 M ) time complexity, where N M is the discretization steps. In this section, we will mainly discuss the rough Heston model as well as its characteristic function [11].
We define the rough Heston model as follows: Definition 1 (Rough Heston model). Suppose that S t is the stock price at time t, X t is log-stock price X t = log(S t ), and V t is the variance process at time t, then the rough Heston model follows the stochastic process: and where ρ is the correlation of stock return and variance process, θ is variance's mean reversion level, λ is the magnitude of pull on the variance's mean reversion, ν is the magnitude of variance's stochastic movement, H is the Hurst parameter (used to control roughness of variance process), and W 1 , W 2 are independent Brownian motions in filtered probability space (Ω, F , Q) with (F t ) t≥0 satisfying the usual conditions (complete and right-continuous). Please note that F := F W 1 ∨ F W 2 where F W 1 and F W 2 are the filtrations generated by W 1 and W 2 respectively, and Q is the risk neutral probability measure. Definition 1 is the macroscopic level of the microscopic Hawkes process [10]. We can actually shorten the variance process of rough Heston model to a stochastic equation that consist of forward variance curve by letting θ = θ 0 (u) and satisfies some constraints [24,25]. Equation (2) in terms of forward variance curve is defined as follows: where ξ 0 (t) is the expected variance at time t conditional on information on time zero. We will use Equation (3) instead of Equation (2) as its characteristic function is in simpler terms and used in many other studies [7,26,27]. We move on to the discussion of the characteristic function of rough Heston model. It is derived by El Euch and Rosenbaum [11] and further developed (forward variance form) in the study [24]. In this study, we denote E[·] as the expectation with respect to probability measure Q. The characteristic function of rough Heston model in terms of forward variance curve is as follows: (3), then the characteristic function of rough Heston model is

Theorem 1 (Characteristic function of rough Heston model). Consider the rough Heston model from Definition 1 and Equation
where i is imaginary number, ξ 0 (u) is the forward variance curve, and D α h(a, T − s) is fractional Riccati equation in the form of: Note that h(a, t) is the unique continuous solution of Equation (5).
Now that we have the characteristic function of rough Heston model, all we must do for computing the option price is to plug in the characteristic function to Lewis' method [6,15]. Recently, a study by Baschetti et al. [28] has noted that SINC approach [29] is more suitable to price option under rough Heston model. The efficiency of the SINC method triumphs Carr-Madan's method [14] and Lewis' method [15].

Option Pricing Methods under Rough Heston Model
In this section, we discuss about the Padé method [20,22] and option approximation under rough Heston model [23]. As compared to fractional Adams method, these methods can efficiently approximate option price under rough Heston model. In particular, instead of time complexity O(n L N 2 ) needed to price option using fractional Adams method, Padé methods [20,22] only require time complexity O(n L ) and the approximation formula [23] require time complexity O(n v ) (n v is number of size steps for simple integration function of forward variance curve; if flat variance curve is used, it only require time complexity O(1)). Please note that the time complexity O(1) refers to the constant time needed to execute or compute the algorithm/solution; as for our study, methods that have time complexity of O(1) are computed almost instantaneously (less than 10 milliseconds).

Padé Approximants
We introduce the generic definition of Padé approximant as follows: where p z and q z are variables to be solved.

Remark 1.
In the study [20], the authors have noted that only diagonal Padé approximant is admissible for fractional Riccati equation's solution, in other words, the term i is equivalent to term j in Definition 2. For simplicity purposes, we call it k-th order Padé approximant, where k = i = j.
Let us give the overall process of how the method works. Our main aim is to approximate the solution of fractional Riccati equation as follows: where h(a, t) is defined in Equation (5). The fractional Riccati equation actually possesses an infinite series expansion solution [27], and the solution will be used as short-time series expansion. Furthermore, in the study [20], the authors have derived the long-time series expansion. Now, the p z and q z are solved in the following linear equations: and where h s (a, t) and h l (a, t) are the short-time and long-time series FRE solution respectively, and the values (m, n) are the truncated terms of order (m, n). Let us call the third-order Padé approximant introduced by Gatheral [20] as GPadé and fourth-order Padé approximant [22] as FPadé. Please note that the GPadé and FPadé are constructed differently (see [22]). In particular, we can construct three different FPadé, i.e., FPadé (m, n) = (5, 4), (6, 3), (7,2). The full steps to obtain the p z and q z for z = 0, 1, ..., k are outlined in the study [22]. Now, we can compute for fractional Riccati equation's solution in O(1) time complexity using the Padé methods instead of O(N 2 ) through the fractional Adams method. The use of this Padé approximation will greatly reflect in the time needed to compute for option price under rough Heston model.

Approximation Formula of Option Price under Rough Heston Model
Unlike the Padé approximation method, the approximation formula derived in this section does not require the Lewis' method [15] to compute for the option price, it is a standalone equation. The approximation formula is obtained through the decomposition of options price under the rough Heston model. We introduce some definitions and notations before touching upon the decomposition theorem.

Definition 3.
Suppose that V t is the variance in rough Heston model (Definition 1), we define the future expected variance as follows: We define a similar equation called the "martingale": From Definition 3, we can establish the following relationship between w t (T) and M t as: and We will assume zero interest rate in this section as it is easier to deal with the equations. For the inclusion of interest rate, readers may see the study [30,31].
The decomposition formula will involve the Black-Scholes call option pricing equation, therefore we introduce it as follows: where w t := w t (T), and d 1 (a), d 2 (a) are defined as and Using the Black-Scholes pricing equation V(X t , w t ), we can define the decomposition of option price under rough Heston model as follows: Theorem 2 (Decomposition formula). Let V t := V(X t , w t ) from Equation (14). Then, the decomposition formula of an option price for t ∈ [0, T] can be written as: where X is the log-stock price, M is the martingale (Definition 3), the partial derivative ∂ ww = ∂ 2 ∂w 2 , and another ∂ xw = ∂ 2 ∂xw .
Proof. The proof can be found in [23,27].
Theorem 2 describes the option price under a certain process is governed by conditional expectation of two integrals that involve partial derivatives of future expected variance and log-stock price. Although now we have a decomposition formula for option price under a stochastic volatility process (including rough volatility model), the decomposition formula is still computationally expensive. Therefore, we will derive an approximation of the decomposition option pricing formula based on rough Heston model under an assumption. We assume the following inequality to hold when deriving the approximation: where E u (·) := E(·|F u ) is the conditional expectation on filtration F u . Then, by considering the assumption above, we can derive an approximation option pricing formula under rough Heston model: where Equation (19) is an approximation of the decomposition Formula (17). Notice that the first term of the right side of Equation (19) is the Black-Scholes call pricing equation, the second and third terms can be computed in O(1) or O(N a ) time complexity depending on whether flat or market forward variance curve is used. Either way, Equation (19) can be efficiently computed. The error of the approximation Formula (19) is discussed in the study [23]. Similar to the option price under rough Heston model, its approximation formula has the capability of generating term structure of at-the-money skew of order O(T H−1/2 ). The approximation formula shows a great resemblance of the option price under rough Heston model after a calibration [23].

Numerical Experiment
In this section, we will deal with the numerical experiment of Padé methods and the approximation formula for option price. The numerical experiment will be carried out in two vastly different scenarios, i.e., the first scenario is on 18 June 2021 (a normal day) and the second scenario is during the Lehman collapse event (12 and 15 September 2008). We use the parameterized forward variance curve of the rough Heston model (Definition 1): where E α,β (·) is the generalized Mittag-Leffler function [22]. In addition, let OGPadé and OFPadé be the option price under rough Heston model where the fractional Ric-cati equations in their characteristic functions are approximated by GPadé and FPadé respectively.

Case of 18 June 2021
The first calibration will be carried out on an ordinary day-18 June 2021. We use the option for the top 500 companies in US market, i.e., the S&P 500 option or SPX option. Four different maturities ranging from T = 0.5 to T = 1.99 will be considered; the calibration of the parameters will be based on the minimization of root-mean squared error of the implied volatility. The methods GPadé, FPadé (5,4), (6,3), (7,2), and the approximation Formula (19) will be considered for the numerical experiment.
Let us interpret Table 1. The option price computed by Padé methods are calibrated to similar parameters. It can be seen that Hurst parameters are low; this indicates that the variance of the stock return has strong dependency on its past variance movement. The low Hurst parameter H observed on the OGPadé and OFPadé would mean that they are very consistent to the rough Heston's option price computed by fractional Adams method [20,22], in other words, it is no coincident that the rest of the parameters (besides Hurst parameter H) are similar too. If the calibrated Hurst parameter H is actually large, for instance near 0.5, the other calibrated parameters of OFPadé method are very likely to differ from OGPadé's parameters. In addition, the parameters of the approximation Formula (19) are quite different from Padé methods, this is mostly due to the approximation error of option price under rough Heston model.
The Figure 1 shows the performance of OGPadé, OFPadé (5,4), (6,3), (7,2), and approximation Formula (19). A close inspection indicates that OGPadé and OFPadé (5,4), (6,3), (7,2) are extremely identical to one another. On top of that, the Padé methods perform extremely well, i.e., matching the SPX implied volatility surface across four different maturities almost perfectly. On the other hand, the approximation formula (gray asterisk marker) generates matching implied volatility curves of the SPX implied volatility, but it clearly has poorer performance as compared to the option price computed through Padé methods.

Case of Lehman Collapse Event
By definition, subprime is referring to a loan for an individual that has below-average credit rating. The rising of U.S. home prices in the early 2000s has led to an increase of individuals or groups taking out mortgages to buy houses, this has led to a housing bubble. In September 2006, U.S. home prices dropped the first time in 16 years, this event subsequently led to a large ripple effect soaring through the risk takers which would eventually affect banks that issued the mortgages. The large exposure and position in the mortgage market by the Lehman Brothers caused them to declare the largest ever bankruptcy in U.S. history on 15 September 2008.
In this section, we calibrate five different methods, namely the OGPadé, OFPadé (5,4), (6,3), (7,2), and the approximation formula (19) to the SPX option before and after the Lehman Brothers collapse (12 and 15 September 2008). The calibration approach will be the same as the case of SPX option on 18 June 2021. One thing to take note specifically is that we are using option with the same maturity date, but with different starting date, therefore it has different time to maturity, 1 day lesser to be precise.    From Tables 2 and 3, we first discuss about the methods OGPadé and OFPadé (5,4), (6,3), (7,2). It is noticeable that all parameters of OGPadé and OFPadé are extremely identical on its calibration date. Comparing the parameters on 12 and 15 September 2008, we can first infer from the increase in Hurst parameter H and the increase in "volatility of volatility" ν such that the reliant on past volatility of stock return has decreased, but the anti-persistency effect still remains. The correlations of variance and the stock return ρ have decreased and the initial variances V 0 have variance have also increased; this is expectable as a black swan event has just happened. As for the mean reversion variance level θ, it has decreased, but it should remain unchanged or differ by a small amount. The speed of variance heading to mean reversion level λ has just increased by a small amount. The changes of parameters of the approximation formula (19) before and after the Lehman Brothers collapse actually behave in the same order as the OGPadé and OFPadé methods, but one parameter in particular (λ) has changed largely. We will further discuss this change in the next paragraph.
Let us observe the Figure 2. Starting with before the Lehman Brothers collapse (12 September 2008, dashed red line), we can see that the methods OGPadé, OFPadé (5,4), (6,3), (7,2), and the approximation formula (19) are quite consistent to the SPX implied volatility. Right after the Lehman Brothers collapse (15 September 2008, dashed blue line), the implied volatilities of short-expiry SPX option (T = 0.10, 0.51) of SPX option have shifted upwards significantly when compared to before the collapse, whereas the implied volatilities of long-expiry SPX option (T = 1.26, 1.75) have only slight alteration (near unchanged). As for the performance of OGPadé, OFPadé (5,4), (6,3), (7,2), and approximation formula, there is a decrease in compatibility for matching the SPX implied volatility. In particular, while the overall shape of calibrated implied volatilities of the methods are there, the calibrated implied volatilities of the methods appear to have decreased too quickly as time increases. Notice that this phenomenon is likely the reason the parameter λ of the approximation formula (19) has increased significantly. Overall speaking, OGPadé and OFPadé performed better than approximation formula in terms of performance for non-extreme and extreme market condition, but they do not differ by much.

Conclusions
This study is mainly filling in the absence of numerical experiment of some methods [22,23] on real option data. It mainly provides a review and numerical performance of five different methods on SPX options in different scenarios. In Section 1, we give some literature review of the current development of stochastic volatility models, and in particular, the rough Heston model. In addition, we have discussed some available methods to compute for the characteristic function of rough Heston model. Furthermore, we have discussed the definition of rough Heston model and its characteristic function in Section 2. Then, in Section 3, we move on to the option pricing methods of rough Heston model, where we have briefly introduced and explained five approximation methods, i.e., OGPadé, OFPadé (5,4), (6,3), (7,2), and approximation formula (19).
We perform numerical experiment in Section 4. The numerical experiment is tested on two different scenarios, a normal day and a major catastrophic event, the Lehman Brothers collapse. It is found that under non-extreme market condition, the OGPadé, OFPadé (5,4) (6,3), (7,2) performed similarly to one another and it has better performance than the approximation formula (19). In extreme market condition, the methods OGPadé, OFPadé, and approximation formula have reduced accuracy in matching the SPX implied volatility, but nevertheless, they are still capable of mimicking the SPX implied volatility surface.

Acknowledgments:
The article was partially funded by the Geran Putra Berimpak (GPB) with project number GPB/2017/9543000.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The