Fractional Riccati Equation and Its Applications to Rough Heston Model Using Numerical Methods

: Rough volatility models are recently popularized by the need of a consistent model for the observed empirical volatility in the ﬁnancial market. In this case, it has been shown that the empirical volatility in the ﬁnancial market is extremely consistent with the rough volatility. Currently, fractional Riccati equation as a part of computation for the characteristic function of rough Heston model is not known in explicit form and therefore, we must rely on numerical methods to obtain a solution. In this paper, we will be giving a short introduction to option pricing theory (Black–Scholes model, classical Heston model and its characteristic function), an overview of the current advancements on the rough Heston model and numerical methods (fractional Adams–Bashforth–Moulton method and multipoint Padé approximation method) for solving the fractional Riccati equation. In addition, we will investigate on the performance of multipoint Padé approximation method for the small u values in D α h ( u − i/2, x ) as it plays a huge role in the computation for the option prices. We further conﬁrm that the solution generated by multipoint Padé (3,3) method for the fractional Riccati equation is incredibly consistent with the solution generated by fractional Adams–Bashforth–Moulton method.


Introduction
The classical stochastic volatility diffusion models such as [1,2] have played an important role in volatility modelling specifically to match for the smile or skew effect (variation of implied volatility with respect to strike price) as displayed in the empirical implied volatility [3]. However, the shortcoming of using these stochastic volatility models is their inability o generate reasonable shape for the term-structure (variation of implied volatility with respect to maturity time) in the implied volatility. Generally speaking, the classical stochastic volatility diffusion models have a skew slope of O(1), whereas the implied volatility based on empirical observations has a skew slope of O(τ −1/2 ) [4][5][6] with τ as the time to maturity, therefore they are not compatible as market models. In response to the shortcoming, the stochastic volatility jump diffusion models such as Bates model [7] were also introduced to fit the empirical short-time skew. Moderate but insufficient improvement was made to the problem by using the stochastic volatility jump diffusion model.
Recently, the academic contribution by [8] enables us to price assets accurately by considering rough volatility in our framework. The term "rough volatility" indicates that the volatility term is driven by fractional Brownian motion (fBm) [9] with Hurst exponent of 0 < H < 1/2, where H controls the roughness of the fBm (the lower the H, the rougher the trajectories of fBm). Intuitively, the term "rough volatility" came from its roughness behavior or anti-persistency effect. Accordingly, The usual algorithm dealing with fractional order differential equation is obtained by using the fractional Adams-Bashforth-Moulton method [28]. The error analysis is also provided in the paper. We refer to the fractional Adams-Bashforth-Moulton method as fractional Adams method from now onwards. It is noted that [14] has also demonstrated the use of fractional Adams method in its numerical application section for the rough Heston model. Although it is capable of providing extremely well solution with large time steps, the algorithm complexity to evaluate a call option is at O(N a n 2 ) where the N a is the number of space steps used for Fourier-type method and n is the number of time steps for the fractional Adams method. The sole complexity of the algorithm makes the computation of the call option to a certain accuracy not feasible to most practitioners. We will take advantage of this method by setting a high number of time steps to compare our numerical solution later in this paper.
In this paper, we will also focus heavily on the work of [29] where it uses multipoint Padé approximants on the asymptotical solutions (t → 0 and t → ∞) of the fractional Riccati equation and applying it to the rough Heston model. This is needed as the typical method (fractional Adams method) requires great computational effort, i.e., not feasible to most practitioners or perhaps even researchers. Before that, we wish to discuss some of the general aspects and history of Padé approximation method. Specifically, it was developed by Henri Padé in the 1890s with the constructing idea credited to Georg Frobenius where he investigated the usage of rational approximations of the power series. Coincidentally, Padé approximation turns out to be a great approximation to many application in physics, e.g., nuclear physics [30], kinetic electron model [31], quantum resonances [32] and many more [33]. We follow closely on the literature from [29], we note that the theory of convergence of Padé approximants can be found in [34]. Also, as noted in [33], one of the theorems (Montessus de Ballore's theorem) states that if f is analytical in a ball centered at zero except for multiples of n, the sequence of { f (m,n) } ∞ m=1 for a fixed n, converges uniformly to f in compact subsets omitting the poles. Furthermore, since we are dealing with diagonal Padé approximants for our fractional Riccati equation later on, it is only fair that we mention its convergence properties too. As noted in [35], Baker-Gammel-Wills conjecture speculates that under some certain condition, the diagonal Padé approximant { f (n,n) } ∞ n=1 converges uniformly to f . Taken from [34], a counterexample for the Baker-Gammel-Wills conjecture is provided. Consequently, Nutall-Pommerenke's theorem is provided where it states that there exists a subsequence of { f (n,n) } ∞ n=1 where it converges almost everywhere. Coupled with Monstessus de Ballor's theorem, we should be extremely careful when dealing with diagonal Padé approximant in our work. One of the functions we will mention later is the Mittag-Leffler function E α (z). The uniform convergence of Mittag-Leffler function has been proven in [36] on the compact set {|z| ≤ 1}. However, it is not necessarily compatible with the asymptotic behavior of the Mittag-Leffler function for large t, although it performs better than their corresponding truncated Taylor series.
We now move to the discussion on the multipoint Padé approximation method. The method was originally derived in [37] as a global rational approximant. As mentioned previously, it aims to match asymptotical points of f using (m, n) Padé approximant. In addition, it has been successfully applied to both the Mittag-Leffler function [38] and generalized Mittag-Leffler function [39]. The authors of [29] demonstrated the use of global rational approximations for the fractional Riccati equation and found the excellent performance of the multipoint Padé approximation (3, 3) especially when H is close to 0. However, the imaginary part of the solution deteriorates as H → 0.5 or α → 1. This poses a challenge in terms of accuracy of the model and it would be erroneous to employ the same method to approximate the solution for the fractional Riccati equation when H is somewhat near 0.5. In addition, we take a step further to investigate the accuracy of D α h(u − i/2, x) for small u values as it plays an important role for the computation of option price. Consequently, we further confirm that the multipoint Padé (3,3) approximation method performs incredibly well at small u values as compared to fractional Adams method, thereby enabling practitioners to price their options using a low computational cost method with similar accuracy.

Organization of the Paper
We have provided a short introduction to the Riemann-Liouville fractional calculus in the Appendix A and Mittag-Leffler function in Appendix B as preliminaries for interested readers. Furthermore, we intend to organize our paper as following: Section 2: Fractional Adams-Bashforth-Moulton method and its error analysis. Section 3: Black Scholes equation, implied volatility, classical Heston, rough Heston model, their characteristic functions and their connection to call option pricing. Section 4: Small and long time expansion of solution for the fractional Riccati equation. Section 5: Padé, multipoint Padé approximation method and its application to the asymptotical solutions of the fractional Riccati equation.
Section 6: Numerical experiments and performances.

Fractional Adams-Bashforth-Moulton Method and Its Error Analysis
We discuss the general method of obtaining a solution for the fractional ordinary differential equations in which we will utilize this method to obtain the solution for the fractional Riccati equation later on. In particular, we refer to [40] for the fractional Adams-Bashforth-Moulton method. The original work on the error analysis belongs to [28]. Furthermore, the full algorithmic description with pseudocode can be found in [41]. It is also noted that this method is also known as one of broader method called predictor corrector method or PECE (predict-evaluate-correct evaluate). This method relies on product trapezoidal method and product rectangle method in the quadrature theory. The fractional Adams method with fractional order of α is shown as following with a j,k and b j,k as In particular, we can let N = k as the total of time steps, pick T and ∆ such that N = T/∆, where T is the total time and ∆ is the step size (usually extremely small). Subsequently, the D α h(a, x) in Equation (1) will be our interest of study-fractional Riccati equation. It is noted that the complexity of this algorithm is at O(N 2 ). In the later section, we will discuss why employing this algorithm for the rough Heston model in the option pricing framework has higher computational complexity than O(N 2 ).
We follow similarly from [14] and note that [42] provides convergence proofs of fractional Adams method. We denote h j as the solution obtained using fractional Adams method and h(x j ) as the real solution. In particular, for x > 0 and a ∈ R and max for any ε > 0.

Option Pricing Models, Implied Volatility, Characteristic Functions of the Option Pricing Models
We discuss some of the financial related models and their aspects in this section. We intend to organize them as following 1. Introduce the famous Black-Scholes pricing model and discuss about implied volatility as well as its importance in the financial market. 2. Focus on one of the most famous financial model amongst the practitioners-classical Heston model and the newly developed rough Heston model. 3. Display its characteristic functions and its connection to call option pricing formula using the inversion of characteristic function.

Black-Scholes Model and Implied Volatility
The Nobel Prizes of economics were awarded to the Myron S. Scholes and Robert C. Merton in 1997 for their extremely marvellous derivation of the now known as Black-Scholes equations. Sadly, one of the authors Fisher Black passed away before receiving the prize. In this subsection, we briefly discuss about the option pricing equation [43] from the Black-Scholes equation derived under some conditions and the implied volatility. The importance of the discussion for Black-Scholes model is at its implied volatility computation. Although we will not display the usage of the implied volatility in this study, we believe that some readers would find its mechanism interesting. In particular, Black and Scholes derived the Black-Scholes equation using two different methods, i.e., constructing a hedged portfolio and the capital asset pricing. Both of which led to the same equation. We now provide the closed-form solution for the European call option: where S is the stock price, K is the strike price, σ is the volatility of the stock return, r is the interest rate, T is the time to maturity and N(·) is the cumulative normal density function. The short and concise Equation (9) plays an important concept in the financial industry; literally everyone who trades stock and options would have known or heard about it. However, some of the assumptions that led to Equation (9) are unrealistic. One of them is the assumption of constant volatility in which we intend to replace it with rough volatility in the rough Heston model. Nevertheless, we will discuss why Equation (9) is still popular among the practitioners. Although no practitioners would use Equation (9) to price their options in this era, many of them took interest on its implied volatility σ imp . The implied volatility has become the standard tool to evaluate the volatility surface of a particular stock. Volatility surface is a combination of volatility skew/smile and volatility term structure, normally seen as a three-dimensional graph for comparison purposes. Let the F m as the market price of the model or price of other option model, then with BS −1 (F m ) denoting the inverse of market price on the Black-Scholes function in Equation (9). While there is no explicit closed form solution for obtaining Equation (11), the monotonic relationship of the Black-Scholes pricing model and the volatility σ makes the implied volatility extremely easy to obtain. The monotonicity of σ can be seen from one of the greeks of Black-Scholes pricing model-vega where C is the call option price in equation (9). In particular, bisection method is one of the easy and fast way to obtain the implied volatility σ imp . What practitioners normally do with the implied volatility σ imp is that they compare the figures across different strike price K and maturity T in between the volatility implied by the financial market prices and the model prices. As a result, it is easy to observe whether the model is capable of fitting the empirical data and reconstruct the surfaces of the market prices.

Classical Heston Model and Rough Heston Model
We now turn to the discussion of Heston models, both classical and rough version. Suppose from [2], we let some security or stock price S and its instantaneous variance V = σ 2 such that S and V follow the stochastic process: where λ ≥ 0 is the magnitude or speed of mean reversion level, θ ≥ 0 is the mean reversion level for the volatility, ν is the volatility of the volatility of the model and finally ρ ∈ [−1, 1] is the correlation between B t and B ⊥ t as the independent Brownian motions. From [14], we discuss why the classical Heston model is popular among the practitioners. 1. The model reproduce several stylized facts of low frequency stock data, e.g., the leverage effect, time-varying volatility and fat tails. 2. It generates similar shapes and dynamics for the implied volatility surface. 3. Efficient computation for the classical Heston model using the explicit formula for the characteristic function of the asset log-price (we will discuss it later).
Subsequently, we discuss about the the rough volatility model shown in [14], the one-dimensional asset price S is in the form of with where λ ≥ 0, ρ ∈ [−1, 1] is the correlation between spot and volatility movement, Γ(·) is the Gamma function, ν is the volatility of the volatility, H is the so called Hurst parameter with the connection of α = H + 1/2 from the Appendix A and θ t (s) is the F t -measurable mean reversion level that makes the model time consistent [44]. We note that the third term in Equation (15) is the fractional Brownian motion with Hurst parameter H, when H is between (0, 1/2), the rough Heston model follows.
In particular, the Hurst parameter controls the roughness of Brownian motion. When 0 < H < 1/2, the fractional Brownian motion displays a rougher state or anti-persistency of fluctuating pattern as compared to the classical Brownian motion. Furthermore, fractional Brownian motion possesses memory property, i.e., future observation depends on the past observations, a feature seen in the financial stock market. As a result, we note that fractional Brownian motion with Hurst parameter 0 < H < 1/2 matches empirical volatility in the financial market better than the classical Brownian motion, therefore it is advantageous to use the rough Heston model rather than classical Heston model for modelling the dynamics of assets. Visit ( [8], Figure 15) for a figure illustration on the comparison of fractional Brownian motion and the market volatility. Not surprisingly, when H = 1/2, the rough Heston model will revert back to the classical Heston model. In comparison to Equation (15), under the perfect hedging scheme as shown in [44], Equation (15) can be rewritten in the compact form of where (ξ t (u)) u≥t is the forward variance curve observed at time t. This is because λθ t (s) can be inferred from the forward variance curve, therefore λ is set to 0 and rewritten as the forward variance form of Equation (16). Furthermore, from [20,45], it is possible to obtain the forward variance curve by differentiating the variance swap curve. In particular, by assuming continuous sample paths, the well-known fair value variance swaps can be obtained from an infinite log-strip of out-of-money options. The importance of Equation (16) cannot be emphasized enough, reduction of parameter estimation for λ and θ t (s) as well as an integration term help immensely in the computation or calibration aspect of the rough Heston model. Note that the absence of λ in the second term of Equation (16) is due to a different way of expressing the notation between different authors, and we will be sure to specify the appropriate form of the characteristic function in the next section accordingly.

Characteristic Functions and Their Connection to Call Option Pricing
In this subsection, we discuss both of the characteristic functions for classical Heston and rough Heston model. Subsequently, we show their connection to call option pricing. Although the characteristic function was derived in [2], we refer to [14] for an overall review for the classical Heston model. We let the characteristic function of the log-price X t = log(S t /S 0 ) as φ t (a, T), then it is given by and g(a, t) = θλ t 0 h(a, s)ds (18) where h is the unique solution to the following Riccati equation: where the parameters λ, ρ, ν and θ are taken from Equations (12) and (13). The parameter a ∈ C will be used in inversion of the characteristic function later and i is the imaginary number. As stated before, Equation (19) has an explicit closed-form solution which can be easily obtained. From page 18 of [45], the solution is given as with Eventually, the closed-form solutions from Equations (20)-(22) enable us to price the call option price with little effort. The exact pricing formula using the inverse Fourier method will be discussed after the introduction on characteristic function for rough Heston model.
Surprisingly, the rough case of Heston model (we are referring to Equation (15) first) that corresponds to the version of Equation (19) differs with little modification to its Riccati form. The only difference is that the classical Riccati equation is replaced with a fractional Riccati equation. In a precise form, the characteristic function of the rough Heston model is as following and where α = H + 1 2 , θ and λ are taken from Equation (15), whereas h is the solution for the fractional Riccati equation with the following form: with the fractional derivative D α and integral I 1−α discussed in Appendix A. All the parameters λ, ρ and ν are taken from Equation (15). Note that when α = 1, Equation (25) reverts back to Equation (19) and the rough Heston model coincides with the classical Heston model. Unfortunately, there is no known explicit closed-form solution for the fractional Riccati equation, i.e., no closed form solution when α < 1. Therefore, the only available choice is to rely on numerical methods. In addition, we wish to discuss another version of rough Heston model based on the work of [44] which is based upon Equation (16). We let x α = νt α as a dimensionless quantity and from [46], the characteristic function of Equation (16) can be denoted as following where h(a, ·) as the unique continuous solution of the fractional Riccati equation with the form of where all the parameters have the same definition as from Equation (16). Note that we are using X t = log(S t ) in Equations (26) and (27) as compared to X t = log(S t /S 0 ) in Equations (17) and (23). From now onwards, we will focus on this second version of characteristic function for rough Heston model, i.e., Equations (26) and (27).
We have now discussed the characteristic function of both the classical Heston model and rough Heston model. Now, we turn to the inversion of the characteristic function which leads to the option price, the key aspect of this work. From [47][48][49], the call option price can be denoted in the form of where S is the current stock price, q is the dividend yield, T is the time to maturity, K is the strike price of the option, r is the interest rate or borrowing rate and k = log( K S ). Now, using Equation (28), it is possible to evaluate the call option price under the classical Heston model immediately since the explicit closed-form solution for the classical Riccati equation is available. As for the case of fractional Riccati equation in Equation (27), the typical method to solve the equation by using the fractional Adams method as shown in Section 2. The fractional Adams method has been widely employed in [14,20,29].

Small and Long Time Expansion of Solution for the Fractional Riccati Equation
The short and long time expansion of the solution h(a, x) with x α = νt α are the key components for our efficient approximation method (Padé approximation method) for the solution of fractional Riccati equation which we will discuss in Section 5. Therefore, we take this opportunity to discuss on the expansions in this section.

Small Time Expansion on Solution of Fractional Riccati Equation
We first discuss on the work of [46], the short-time x → 0 expansion of h(a, x). In particular, h(a, x) can be represented as where Γ(·) is the Gamma function. Although Equation (31) appears as a complicated set of summation, we note that only first few terms will be used in the approximation method which we will discuss in the next section. Earlier terms of Equation (31) are fairly quite manageable and easy to obtain. We wish to emphasize that the combination of small time expansion and the Padé approximation method are extremely effective and consistent in obtaining the small time solution for the fractional Riccati equation later on.

Large Time Expansion on Solution of Fractional Riccati Equation
Next, we obtain the following work from [29]. Under the expectation for lim x→∞ h(a, x) = r − from classical Riccati solution, we expect the solution of fractional Riccati equation behaves the same way. Consequently, we approximately let h(a, x) to be r − in one of the term in the fractional Riccati equation as the large time expansion. Subsequently, the authors managed to derive out the large time expansion for the solution of the fractional Riccati equation. We begin by noting that Equation (27) can be rewritten as with In linear form, the fractional Riccati equation is approximated as following The solution to Equation (34) is of the form where E α (·) is the Mittag-Leffler function denoted in Appendix A. Equation (35) is apparent from the derivation of Equation (A6). We take a proposition from [29] as: Proof. See proposition 3.1 in [29] with the help of corollary A1.
From Equation (A8), it motivates the authors of [29] to derive the following result: as x → ∞ along with γ 0 , γ 1 , γ 2 in the forms of (38) and the general recursion formula as where Equations (37)-(39) are obtained by matching coefficients using Equation (32). Equation (36) is our large time expansion for the fractional Riccati equation. We now have the necessary tools to employ the approximation method on fractional Riccati equation.

Multipoint Padé Approximation Method for Fractional Riccati Equation
As we have mentioned before, the main issue for the rough Heston model is at its high computational cost. Coupled with fractional Adams method, in order to compute for the option price, a computation complexity of O(N a N 2 ) is needed, where N a is the number of space steps for Equation (28) and N arises from the use of fractional Adams method noted in Section 2. One of the important advancements for rough Heston model belongs to the use of multipoint Padé method for obtaining the fractional Riccati equation's solution at O(1) complexity. This is a huge leap in computational aspect as [29] has shown virtually that N ≈ 20,000 would result to a satisfactory solution indifferent from the N ≈ 50,000 case. By employing the multipoint Padé method, the computational complexity of the evaluation for option price using Equation (28) reduces to only O(N a ), the same as the case of the classical Heston model. We first introduce the work of [38] which concerns with the global rational approximation on h(a, x): h (m,n) (a, x) = ∑ m i=1 p m y m ∑ n j=0 q n y n (40) where m, n ∈ N + and we make the substitution y = x α on both Equations (29) and (36). Based on the argument for [29], we note that only diagonal Padé approximation is admissible. This is due to the fact when x → ∞, Padé approximation should approach as following where 0 < |r − | < ∞ only if m = n. Notice that if m = n and x → ∞, h m,n (a, x) will either not admit constant other than zero or does not exist (involves ∞ or −∞). Furthermore, Equation (41) corresponds to an assumption we made at Section 4.2. The possibilities of choosing m = n equals to some constant is endless. However, the authors [29] have noted that m = n = 2 and m = n = 4 do not perform well for the fractional Riccati equation case. Fortunately, the case m = n = 3 does produce excellent quality of solution for the fractional Riccati equation, at least for the H → 0 solution. Therefore, we will work out the case for multipoint Padé approximation method for m = n = 3 in this section. When we mention multipoint Padé approximation method, we intend to refer it as multipoint Padé approximation method using m = n = 3 for the rest of the discussion in this paper. We follow closely from the presentation of [29] on this section. Remember that we have set y = x α for the rest of our discussion. Suppose that from Equation (29), we obtain the small time expansion in the form of and from Equation (36), we have the series expansion for small time on the solution of fractional Riccati equation as The Padé (3, 3) approximation can be written as h (3,3) (y) = p 1 y + p 2 y 2 + p 3 y 3 1 + q 1 y + q 2 y 2 + q 3 y 3 (44) with q 0 set to 1. Matching the coefficients to Equations (29) and (36) using Padé approximation method, we have the system of equations as Accordingly, we obtain the constants p 1 , p 2 , p 3 , q 1 , q 2 , q 3 as Note that the solutions p 1 , p 2 , p 3 , q 1 , q 2 , q 3 can be easily obtained using solve function in MATLAB in less than one second. Now that we have the required multipoint Padé method, it is time we test it against the general approach (fractional Adams method).

Numerical Experiment And Performances
In the previous sections, we have demonstrated the fractional Adams method and the multipoint Padé approximation method. Therefore, it is time that we conduct some numerical experiment on the fractional Riccati equation. In particular, we have decided to compare the result of the fractional Riccati equation using fractional Adams method and the multipoint Padé approximation on a different scale of the Hurst parameter H.
The parameters used in our numerical studies are as following unless specifically stated otherwise: where N belongs to the time steps of fractional Adams method. Although we have given the ν value, note that if we use x in x α = νt α as a substitute of ν and t for the fractional Adams method and multipoint Padé method, we will not directly be affected by the value ν is computation for the plotting of D α h(a, x) and x. For perspective, we first display the Re(D α h(a, x)) for H = 0.05 and H = 0.499 in Figure 1. Notice that the solution for the real part of the D α h(a, x) performs extremely well for small H and moderate for H → 0.5. The following observation is taken from [29]. In particular, the deterioration effect is deeply enlarged when we observe for the imaginary part of D α h(a, x) as seen in Figure 2. This in turn would provide an erroneous result for the option price if not fixed. In general, the performance of multipoint Padé approximation method should not be too surprising as we have matched on only two points which are x → 0 and x → ∞. Supposedly, we should also understand the importance of smaller u values in φ T (u − i/2) from Equation (28), smaller u values have higher contribution to the computation of the call option price due to the 1/(u 2 + 1/4) term in the integral. In the subsequent part of this section, we will evaluate D α h(u − i/2, x) in a more robust setting which concerns the choices of u values at an overall scale of Hurst parameter H. We first introduce some of the notations for the purpose of computing errors. We start with the Euclidean Norm Distance of the fractional Riccati equation's solution: (48) and the percentage error of the fractional Riccati equation's solution: where u is our parameter of this study, i is the imaginary number,ĥ(u − i/2, x j ) is the multipoint Padé (3,3) approximation at time x j and h(u − i/2, x j ) is the solution at time x j generated from the fractional Adams method using parameters in Equation (46) with a slight modification of N = 3000 such that x N = 3 (we generate step size of ∆ = 3/3000 = 0.001). We find that N = 3000 is a sufficient time steps for fractional Adams method with reasonable computing time. D α h is from Equation (27) and note that α = H + 1/2 with H as the Hurst parameter. Furthermore, Re(·) and Im(·) denote the real part and imaginary part of the solution. Equations (47) and (48) give an overall error computation of fractional Riccati equation's solution using the multipoint Padé approximation method at different time x j up to x N = 3, whereas Equations (49) and (50). Furthermore, in order to evaluate the quality of D αĥ (u − i/2, x) solution over different values of the Hurst parameter H, we have decided to iterate through different H for our study. Suppose we introduce χ Re , Ψ Re , χ Im , Ψ Im to achieve our final numerical computation result as follows:   [8,20]. As such, multipoint Padé approximation method would be able to perform extremely well with low computing cost. Out of the ordinary, Ψ Re and Ψ Im appear to be large across the values of u in Table 1, and upon investigation this is due to the fact that Re(D α h(u − i/2, x)) and Im(D α h(u − i/2, x)) tend to approach zero as x → 3 during the computation, therefore incurring large percentage error. This is further verified by the consistent low χ Re and χ Im values. In short, the multipoint Padé approximation method performs similarly to the fractional Adams method in low computational errors manner while still retaining the low computational cost advantage.

Concluding Remarks
First of all, we have provided the brief history and literature review for the general stochastic volatility model, rough volatility, rough Heston model and ways to obtain solution for the rough Heston model. In particular, a large part of the literature review is dedicated to numerical methods of solving for the fractional Riccati equation where it can be used to obtain the option prices. Note that the introductory level of the fractional calculus and Mittag-Leffler function are located in the appendix for interested readers. They are subsequently used in describing the dynamics of rough Heston model and represented as the solution of fractional Riccati equation. Later, the fractional Adams-Bashforth-Moulton method and its error analysis are provided. The fractional Adams method still remains a classical and frequently used method of solving any fractional ordinary differential equation. Unfortunately, the fractional Adams method typically requires a large number of time steps which in turn requires large computational cost. But before that, we have introduced some recent advancements in the option pricing theory. Classical Black-Scholes option pricing models, classical Heston model, and the newly rough Heston models have been discussed subsequently. In addition, we gave a slight touch on the implied volatility topic in case readers are keen on testing the models on real data. Furthermore, the characteristic functions and their inversions to the option prices were given. Since the multipoint Padé approximation method requires some expansion on multiple points for matching purposes, we introduce some literature on the small and large time expansion of the solution of the fractional Riccati equation which is being used in the characteristic function of rough Heston model. We later then discuss the required multipoint Padé method (3,3), as suggested by some authors, in a great effort of using the small and large time expansion to obtain the solution on fractional Riccati equation. Finally, we test the multipoint Padé (3, 3) method on a setting and compare it with the fractional Adams method. Lastly, we have further investigated the possible outcomes on u values in D αĥ (u − i/2, x), which is the key component in Equation (28). Consequently, the errors computation of different u values in D αĥ (u − i/2, x) further agrees that multipoint Padé (3, 3) approximation method performs extremely consistent with the fractional Adams method. This enables practitioners to price their options with reasonable confidence such that the solution to fractional Riccati equation obtained using multipoint Padé (3, 3) approximation method will be identical to the ones produced by the fractional Adams method, while enjoying the benefit of low computational cost. Accordingly, what we have found coincides with the previous authors such that multipoint Padé method (3,3) performed extremely well with H → 0 but not the H → 0.5 case. Above all, we wish to make a minor recommendation on this issue-since the solution for the classical Riccati equation or fractional Riccati equation with H = 0.5 exists, in our opinion, we believe it is possible for an attempt of using a hybrid model between the existing multipoint Padé method and the exact solution for the classical Riccati equation. An attempt to match it at H → 0 and H → 0.5 would result in a more robust setting and applicable to altering volatility behavior.

Appendix B. Mittag-Leffler Function
Following on to the discussion Mittag-Leffler function which is a crucial part in Section 4.2 for the derivation of the large time expansion on solution for fractional Riccati equation in a linearisation form. Again from [40], we define the Mittag-Leffler function of one parameter E α (z) for order α ∈ (0, 1] as with z ∈ C. An obvious deduction of E α (Ax α ) is as following A k x kα Γ(kα + 1) (A5) with A ∈ C as an arbitrary constant. We note that Equation (A5) will be used as a direct application to Equation (35). Furthermore, by the relationship of Equation (A3), we obtain A k x kα Γ(kα + 1) We now discuss the asymptotic expansion of Mittag-Leffler function from a review article of [50] originally taken from [51]. The rest of the discussion in this section will focus on the asymptotic expansion and its error, subsequently this discussion is important for the understanding of the large time expansion and its error behaviour in Section 4.2. From Equation (6.5) of [50], we let α ∈ (0, 1] and µ ∈ R such that πα 2 < µ < πα (A7) Then for N ∈ N, N = 1, we have the asymptotic expansion as From Lemma B.2. and Corollary B.1. of [29], we have Lemma A1. For α ∈ (0, 1]. Let a = u + iy with R ≥0 and i as the imaginary number, y ∈ −1/(1 − ρ 2 ), 0 and let A = a(a + i) − ae 2 a 2 , Then for any x ∈ R >0 | arg(−Ax α )| ∈ 3π 4 , π (A9) By setting µ = 3 4 πα in Equation (A8), it follows from Lemma A1 that Corollary A1. For α ∈ (0, 1]. Let a = u + iy with R ≥0 , y ∈ −1/(1 − ρ 2 ), 0 and let A = a(a + i) − ae 2 a 2 , Then for any x ∈ R >0 and N ∈ N