Series Expansion and Fourth-Order Global Padé Approximation for a Rough Heston Solution

: The rough Heston model has recently been shown to be extremely consistent with the observed empirical data in the ﬁnancial market. However, the shortcoming of the model is that the conventional numerical method to compute option prices under it requires great computational effort due to the presence of the fractional Riccati equation in its characteristic function. In this study, we contribute by providing an efﬁcient method while still retaining the quality of the solution under varying Hurst parameter for the fractional Riccati equations in two ways. First, we show that under the Laplace–Adomian-decomposition method, the inﬁnite series expansion of the fractional Riccati equation’s solution corresponds to the existing expansion method from previous work for at least up to the ﬁfth order. Then, we show that the fourth-order Padé approximants can be used to construct an extremely accurate global approximation to the fractional Riccati equation in an unexpected way. The pointwise approximation error of the global Padé approximation to the fractional Riccati equation is also provided. Unlike the existing work of third-order global Padé approximation to the fractional Riccati equation, our work extends the availability of Hurst parameter range without incurring huge errors. Finally, numerical comparisons were conducted to verify that our methods are indeed accurate and better than the existing method for computing both the fractional Riccati equation’s solution and option prices under the rough Heston model.


Introduction
Financial option pricing theory heavily relies on the development of the Black-Scholes model [1], which has several assumptions. One of the unrealistic assumptions assumes that the volatility is constant or non-fluctuating. Due to the mismatch between that expectation and real market prices, it caused a market crash on 19 October 1987-the world's stock market lost more than 20% of its value within a few hours. Due to the event, most investors have realized a lack of volatility that moves stochastically is a major problem. More models have been proposed to incorporate the observed features present in the current stock market.
One of the famous models is the Merton jump-diffusion model [2], which incorporates the stock price jump effect displayed in the financial market, especially when there is an update of the market or a company's condition. Subsequently, the stochastic volatility models such as the Hull-White model [3] and Heston model [4] were proposed. The proposed models suggest that volatility should move stochastically instead of remaining constant. Evidently, the stochastic volatility models have played an important role to match the skew or smile effect displayed by the empirical implied volatility [5]. The skew is the change of implied volatility with respect to the change of strike price, An interesting component of the approximation is that the method uses a n-dimensional classical Riccati equation to approximate the fractional Riccati equation. The authors managed to prove that when n → ∞, the approximation method becomes exact for the fractional Riccati equation. Specifically, the highest dimension of n = 500 was considered in the numerical study. A similar time-efficient method was also proposed by [23] with the introduction of a new model called the lifted Heston model which can be both a classical Heston model or a rough Heston model depending on the parameter n.
There are also methods proposed by [24,25]. The former method approximates the option prices under rough Heston model by substituting a scaled "volatility of volatility" term to the classical Heston model. This directly reduces the computational cost down to the same cost of classical Heston model, which is very appealing to compute for a reasonable accuracy of the call option prices. A hybrid method is proposed as the latter method; it contains Richardson-Romberg extrapolation method coupled with a series expansion of the fractional Riccati equation for approximating the solution to the fractional Riccati equation. The hybrid method is both efficient and flexible as far as the numerical tests have indicated.
In fact, the fractional Riccati equation has been studied for a long time, long before it was discovered that it would be relevant in option pricing theory. The most notable methods are the variational iteration method [26] and the Adomian decomposition method [27]. Furthermore, the homotopy analysis method [28] and the homotopy perturbation method [29] are also popular choices in approximating the fractional Riccati equation's solution. Part of our work also relies on the Laplace-Adomian-decomposition method (LADM) which can be found in [30,31]. Other methods such as [32,33] are also employed to approximate for the solution of the fractional Riccati equation.
Moreover, decomposition formula for a rough volatility model [34] has been proposed and derived based on the work of [35]. The decomposition formulas involves using semi-closed form method to efficiently price European options under the rough volatility models. In particular, the rBergomi model has been successfully tested using a hybrid calibration scheme proposed in [20]. The result of the study is relatively consistent with AAPL options in the financial market. An interesting study done by [36] managed to derive a closed-form option pricing approximation formula for the fractional Heston model (the study has also focused on decomposing option prices as sums of the Black-Scholes formula).
The conventional algorithm to solve fractional ordinary differential equation is called fractional Adams-Bashforth-Moulton method. It is also normally known as the fractional Adams method; error analysis of the method is provided in [37]. Along with the introduction of the characteristic function of the rough Heston model by [17], the authors utilized fractional Adams to obtain the solution for the call option in the numerical experiment section. While it is noted that the fractional Adams method generally provides extremely consistent results with a large number of time steps, the time complexity of the method is at O(N 2 ), where N is the number of time steps of the fractional Adams method. Nevertheless, it has became a standard method to compare consistency on the computation of call option prices for different methods, such as in [19,24]. By the way, the time complexity O(N 2 ) is not the final computational cost of the option pricing formula, as the time complexity is only for the characteristic function. To evaluate the call option pricing, it is necessary to use the inversion of characteristic function formula from [7,38,39]. Therefore, the total computational cost is O(n a N 2 ) to evaluate a single call option price where n a is the number of space steps in the inversion of characteristic function formula. Similarly to other papers, we ought to use fractional Adams method to compare our hybrid method later on.
We move on to the discussion of multipoint Padé approximation method from [19]. The multipoint Padé approximation method matches two extreme points or the asymptotic solutions at t → 0 and t → ∞ using the series expansion solutions of the fractional Riccati equation. Due to the investigation of rational approximations on the power series by Georg Frobenis, it led Henri Padé to discover the Padé approximation method in the 1890s. Noticeably, the Padé approximation method has been used in many areas, such as physics, e.g., quantum resonances [40,41], the kinetic electron model [42] and many more noted in [43]. We utilize the literature given in [19]. The author of [44] provided the theory of convergence on the Padé approximants. In addition, Montessus de Ballore's theorem states that for a fixed j, the sequence of { f (i,j) } ∞ i=1 will converge uniformly to f in the compact subsets, omitting the poles only if the f is analytical in a ball centered at zero except when i is not a multiple of j.
Coincidentally, it is required for us to use diagonal Padé approximants for our multipoint Padé method. Particularly, diagonal Padé approximants refer to Padé approximants whenever i = j. Interestingly, the diagonal Padé approximant { f (j,j) } ∞ j=1 is conjectured to converge uniformly to f under some certain conditions as noted in [45]. However, the conjecture was disproved using a counterexample [44]. Nevertheless, a theorem called Nutall-Pommerenke's theorem implies the existence of the almost universal convergence for a subsequence of { f (j,j) } ∞ j=1 . As of this condition, the authors of [19] took the initiative to try out different n values of the diagonal Padé approximant on the fractional Riccati equation. See [46] for additional details on the methods.
The multipoint approximation method for transcendental functions was originally used in [47] as a global rational approximant. As the name suggests, the main aim of the method is to utilize available points of their expansion series function and match them using the (i, j) Padé approximant. Interestingly, the Padé approximation method and multipoint Padé approximation method have been used to enhance the approximation at a longer interval [48,49]. In addition, the approximations to the Mittag-Leffler function and generalized Mittag-Leffler function using the Padé approximation method were proposed by [50,51] respectively. On a side note, the uniform convergence of the Mittag-Leffler function using Padé approximation function is given in [52]. It is noted that the Padé approximation method performs better than the corresponding truncated Taylor series as usual, but it is not necessarily compatible for large t. This is related to the study from [19] where the authors have utilized global rational approximants in obtaining the fractional Riccati equation's solution by matching two series asymptotic expansions' functions, as the expansion functions are similar to the Mittag-Leffler function and its asymptotic function. Just recently, fourth-order global Padé approximants have been successfully applied to the generalized Mittag-Leffler function by the authors of [53]. The authors managed to approximate the generalized Mittag-Leffler in a highly accurate manner as compared to the second-order global Padé approximants in [51]. Motivated by the study, we employ the methods from [53] in approximating the fractional Riccati equation's solution, and subsequently we were able compute option prices under the rough Heston model in an extremely efficient and accurate manner. This paper is organized as follows. Firstly, Riemann-Liouville fractional calculus will be introduced in Appendix A. Section 2 serves as an introduction to the rough Heston model. In particular, it contains the discussion of the rough Heston's characteristic function and the evaluation of call option prices in terms of the characteristic function. In Section 3, we present a step-by-step approach using the Laplace-Adomian-decomposition method (LADM) to solve the generalized fractional Riccati equation with constant coefficients (mainly motivated by [30]). The details of the LADM on fractional Riccati equation are placed in Appendix B. Section 4 contains some existing work on fractional Adams-Bashforth-Moulton method, and both the small time expansion and large time expansion of the fractional Riccati equation under the rough Heston model by [19,37,54]. Global Padé approximants of third (work of [19]) and fourth order are presented to solve the fractional Riccati equation under the rough Heston model in Section 5. In Section 6, we show the graphical and numerical results of the performance of the fractional Adams method, the third-order Padé approximant and the fourth-order Padé approximants to solve for the fractional Riccati equation, and subsequently use the methods to compute for the option prices. Finally, we provide an overall summary to our work and some recommendations for future research in Section 7.

The Rough Heston Model
Using the ideas from [10,18], the authors from [17] were able to develop the rough version of the Heston model also now known as the rough Heston model. The rough Heston model composes several features of high frequency data in the stock market, e.g., fat tails in distribution, time-varying volatility and the leverage effect. As we have stated from the beginning, rough volatility models are introduced for the purpose of matching the short-time skew of an order of O(T −1/2 ). The numerical result shown in [17] was able to confirm the consistency of the rough Heston model against the empirical financial data.

The Model
We have decided not to introduce the base version of the rough Heston model from [17], instead we will introduce a concise version of the rough Heston model from [55]. How the model from [55] differs from the base model from [17] is that under the perfect hedging scheme, the same authors were able to remove a drift integral component and a part of fluctuating component from the variance dynamics by substituting a forward variance curve with the volatility component. Computation-wise, the authors have successfully removed an additional two parameters in an integral. From [55], we show the stock return dynamics as and the variance dynamics as where V u is the variance at time u, ξ t (u) is the forward variance curve on time u observed at time t with u ≥ t, ρ ∈ [−1, 1] is the correlation between the volatility movement and the stock return, ν is the volatility of the volatility parameter and B t and B ⊥ t are independent Brownian motions. The parameter H is known as the Hurst parameter or Hurst exponent. Hurst parameter H is used to control the roughness or smoothness of the fBm behavior in the integral of Equation (2) with a range of H ∈ (0, 1) \ 1/2. To give a perspective of how it works, the Hurst parameter in the range of 1/2 < H < 1 would cause the denominator term of integral in Equation (2) to favor long dependency on the previous volatility, whereas 0 < H < 1/2 would have the effect of favoring and enlarging the effect on volatility closer to the evaluation time at time u. In other words, the fBm displays a long dependency effect when 1/2 < H < 1 and it displays an anti-persistency effect when 0 < H < 1/2. Specifically, we will limit the Hurst parameter range to 0 < H < 1/2 in Equation (2) since we are using the rough Heston model. By the way, the forward variance curve ξ t (u) can be obtained by differentiating the variance swap curve. More details can be found in [56].

Characteristic Function of the Rough Heston Model
From [54], using the forest expansion method, the characteristic function of the rough Heston model is given as where h(a, ·) is the unique solution of the following fractional Riccati equation: The D α (·) and I 1−α are the Riemann-Liouville fractional derivative and integration, respectively, as introduced in Appendix A. Meanwhile, X t is the log spot-price log(S t ), α is the equivalent to 1/2 < α = H + 1/2 < 1 and other parameters ρ, ν, ξ t (u) have the same meanings/definitions as in Equations (1) and (2).
To obtain the call option price, it is necessary to use the inversion of characteristic function formula derived by [7,38,39]. The call option prices can be evaluated using where S is the stock price, K is the strike price, T is the time to maturity, k is log( K S ) and φ T t (a) is the characteristic function from Equation (3). We have made further assumptions of risk-free interest rate r and the dividend yield q being zero, and t = 0, for simplicity. Before evaluating Equation (5), we have to solve the fractional Riccati equation (Equation (4)) first. As mentioned before, the standard method for solving a fractional ordinary differential equation, such as the fractional Riccati equation (Equation (4)), is computing through an iterative method called the fractional Adams-Bashforth-Moulton method, which we will be discussing further in Section 4.

Analytical Approximation
In this section, we present the Laplace-Adomian-decomposition method (LADM) for obtaining the analytical approximation in terms of infinite series for the generalized fractional Riccati equation's solution with constant coefficients. The method will be applied in a step-by-step manner to ease readers. Readers can go through [30] for the LADM on non-constant coefficients of the fractional Riccati equation.

Preliminaries
The generalized fractional Riccati equation is denoted as where W, R, U ∈ C, α ∈ (0.5, 1), f (t) is the fractional Riccati equation's unique solution and D α (·) is the Riemann-Liouville fractional derivative as introduced in Appendix A. The Laplace transform L(·) of the Riemann-Liouville fractional derivative can be obtained as The generalized Mittag-Leffler function plays an important role in the realm of fractional calculus. We denote the following generalized Mittag-Leffler function as where α, β > 0, z ∈ C and Γ(·) is the Gamma function. Following the work of [30], we denote the kth derivative of the Mittag-Leffler function as Based on the work of [57], we treat one of its results as a lemma: Lemma 1. Suppose we consider the following function: Then, the Laplace transform of the function ε k (t, a, α, β) can be evaluated as and the inverse Laplace transform relationship can be established as Proof of Lemma 1. See pages 6 and 7 of [57].
The next lemma can be further developed from Lemma 1 and from page 210 of [58].

Lemma 2.
For R ∈ C, p > −1 and α > 0, we can obtain the following inverse Laplace transform result: Proof of Lemma 2. The following relationship (see page 210 in [58]) is well known as: Let k = 0 and β = p + α + 1 in Lemma 1; then we can obtain the following result:

Laplace-Adomian-Decomposition Method
In this section, the Laplace-Adomian-decomposition method (LADM) is used to provide a series expansion on the fractional Riccati equation's solution in terms of h(a, t). We provide a summary of the method over here. The details of the method is provided in Appendix B.
The general idea is to apply Laplace transform to Equation (11) which will yield us an easier-looking equation. Then, we employ the idea of the Adomian decomposition method to the solution; i.e., the fractional Riccati equation's solution can be decomposed into an infinite series and the non-linear term can be decomposed into Adomian polynomials. Specifically, the Adomian polynomials can be computed using methods of [31,59,60]. Substituting the series expansions of the solution and the non-linear solution, we then compare and match the left and right side of the equation. In the next step, we can apply the inverse Laplace transform to obtain the terms of the solution. Lemmas 1 and 2 are used to speed up the derivation of the solutions. Repeat the process to construct more terms for the solution. Finally, we substitute the coefficients of the fractional Riccati equation (Equation (4)) into the series expansion. We approach the problem with four distinct steps that are listed in Appendix B. Specifically, the fractional Riccati equation's solution h(a, t) from Equation (4) obtained using LADM is as follows:

The Fractional Adams-Bashforth-Moulton Method and Series Expansion Solutions for the Fractional Riccati Equation
In this section, we present some existing work on the fractional Adams-Bashforth-Moulton method and the series expansion solution for the fractional Riccati equation. The aim is to let the readers familiarize themselves with some of the recent advancements on the computation for the rough Heston model. This section is crucial to developing Section 5 later and will be used in the numerical experiment in Section 6.

The Fractional Adams-Bashforth-Moulton
To solve the fractional Riccati Equation (4), we employ the fractional Adams-Bashforth-Moulton method originally introduced by [61,62] and further explored in [37]. In addition, the pseudocode of the fractional Adams method can be found in [63]. For a reference, we refer to the text [64] for a clearer presentation of the method. In particular, the fractional Adams method is also known as the PECE (predict-evaluate-correct-evaluate) method or the predictor corrector method. In particular, the product rectangle method and the product trapezoidal method from quadrature theory are used to construct the fractional Adams method. Using fractional order of α = H + 1/2, the fractional Adams method is as follows: where a i,k and b i,k are as follows: and Note that h k is the final solution or h k = h(a, t k ) for the fractional Adams method. Notice that the initial condition of h(a, 0) in Equation (20) indicates that the first sums in Equations (18) and (19) are equal to Suppose that we set the total number of time steps as N and the time to maturity as T; then the step size is ∆ = T/N. Moreover, it is apparent that the time complexity of this algorithm is O(N 2 ), which is extremely computational expensive considering it is only for obtaining the rough Heston model's characteristic function.
Originally from [65], the convergence analysis has been done on the fractional Adams method. We follow the short literature from [17] for the result of convergence analysis. Suppose that h(t j ) is the real solution and h j is the solution obtained by using the fractional Adams method. Consequently, for t > 0 and a ∈ R, the error is as follows: and max for any ε > 0.

Expansions for the Solution of the Fractional Riccati Equation
Previously, we have discussed a standard method for obtaining the fractional Riccati equation's solution. Depending on number of time steps N and fractional order α, the error of fractional Adams can be minimized to an extreme precision, as the Equations (23) and (24) have suggested. However, of course, there would be a sacrifice on efficiency in doing so.
Here we introduce work of [19,54] for their small and large time infinite series expansion of the fractional Riccati equation's solution. The purpose of introducing such expansions is that in Section 5, these expansions will be put to use to solve for the terms of global Padé approximants.

Small Time Series Expansion
The solution as an infinite series expansion of the fractional Riccati Equation (4) can be obtained through the forest expansion method by [54]. Specifically, we have where β k (a) is where 1 is the indicator function; ρ and ν have the same meaning from the rough Heston model in Equations (1) and (2); and α = H + 1/2 where H is the Hurst parameter. Noticeably, Equation (25) can be used to construct a small time series expansion for the fractional Riccati equation's solution, just as the authors of [19] did. (25), we obtained an analytical approximation solution to the fractional Riccati equation by using the Laplace-Adomian-decomposition method (LADM). We have verified that the solution (16) obtained using LADM and Equation (25) have the same order of coefficients up to at least fifth order (t 5α ). Since the coefficients of the solution from LADM are the same as the Equation (25) up to the fifth order, we will make use of Equation (25) throughout the rest of the paper, as it is easier to compute for the readers.

Large Time Series Expansion
The large time series expansion of the fractional Riccati equation's solution was derived by [19]. Under x α = νt α , the fractional Riccati equation can be re-expressed as where 1 is the Mittag-Leffler function in Equation (8). Then, for x ∈ R ≥0 and a ∈ A, h ∞ (a, x) solves the Equation (27) up to an error term of O(|Ax α | −2 ) as x → ∞.
Subsequently, Proposition 1 motivated the same authors to provide an ansatz for the asymptotic expansion of h(a, t) as t → ∞ is denoted as where the coefficient γ 0 , γ 1 , γ 2 are denoted as along with the general recursion formula of γ k as where terms r − and A are defined in Equation (28). Note that the authors of [19] used h(a, x) to construct their large time series expansion of the fractional Riccati equation's solution, whereas we have merely transformed it back to h(a, t) in Equation (29) for consistency's sake.

The Global Padé Approximation
The section was directly inspired by [53] for their work on the global approximation to the generalized Mittag-Leffler function. Here, we develop the global Padé approximation to the fractional Riccati equation's solution that has a time complexity of O(1); thus, it is significantly more efficient to compute for the fractional Riccati equation's solution as compared to the fractional Adams method that has a time complexity of O(N 2 ), where N is the number of time steps. Consequently, the fractional Riccati equation's solution obtained by the approximation method can be applied to efficiently compute option prices under rough Heston model using Equation (5). The pointwise approximation errors for the global Padé approximation method on t → 0 and t → ∞ are also included. Existing work on the third-order global Padé approximation by [19] is also mentioned in this section.

Definition
We begin by lettingt = t α ; then we define a function such that the asymptotic expansion forĥ(a,t) → 1 ast → ∞. Specifically, the function admits the following behavior:ĥ where Since when n = 1, then b(t −1 ) = 1, we should always take n > 1 such that we can take advantage of the asymptotic expansions in the global approximation. Subsequently, we move on to the global Padé approximation and its approximation error. Firstly, the x-th order Padé approximation is denoted as where x ∈ N + . Then we denote the global Padé approximation toĥ(a,t) aŝ such that m + n is odd and x := m+n−1 2 ≥ 1. Due to the transformation in Equation (32),ĥ (m,n) (a,t) → 1 ast → ∞ with the Padé approximation approaching limˆt →∞ p(t)/q(t) = p x /q x , we can thus set Using similar argument, we set p 0 = 0 We are now left with 2x − 1 coefficients to solve. We can obtain the coefficients by solving the system of linear equations from the expansion Equations (34) and (35) as In order to yield h (m,n) (a,t) instead ofĥ (m,n) (a,t), a casual inverse transformation gives us where we only consider n > 1, m ≥ 2 and m + n to be odd. Equation (42) is the x-th order global Padé approximation (m, n) to the fractional Riccati Equation (4). For q(t) = 0, we can deduce from Equations (40) and (41):

Approximation Errors
Suppose now we consider the pointwise approximation error between the global Padé approximation method and the fractional Riccati equation's solution fort → 0 andt → ∞. From Equations (33) and (43), we obtain the pointwise approximation error ast → 0 e (m,n) 0 Similarly, we can obtain the following pointwise approximation error fort → ∞ e (m,n) ∞ = O t −n Remark 2. The approximation errors (44) indicate that we should indeed pick m > n − 1 for an accurate small time approximation for the solution of the fractional Riccati equation. Besides that, the approximation error for large time (45) indicates that the approximation error can be arbitrary small by taking a large n value.

Third-Order Global Padé Approximant
The second-order, third-order and fourth order global Padé approximants have been introduced in [19], but the authors did not consider different m and n values in solving the coefficients p 1 , ..., p x and q 0 , ..., q x with x, indicating the x-th order global Padé approximant. The work of [19] indicates that the third-order global Padé approximant with m = 4 and n = 3 performed the best overall as compared to the second-order and fourth-order. Although the authors did not specify which m and n values they used specifically for the fourth-order global Padé approximant, our numerical result shows that they used m = 5 and n = 4. As compared to our work, the authors of [19] have presented their work in a different manner. Their third-order Padé approximant is as follows: Note that the authors set q 0 = 1 and did not set p 3 = q 3 = 1. Specifically, the fractional Riccati equation's solution as the small time series expansion (25) in a simplified form can be denoted as where b 1 , b 2 and b 3 are to match the coefficients of Equation (25). Additionally, next, the fractional Riccati equation's solution as the large time series expansion (29) in concise form can be denoted as where u 0 , u 1 and u 2 are to match the coefficients of Equation (29). Matching the terms of Equations (46), (25) and (29) just like in Equations (40) and (41), we obtain the following sets of equations to be solved: Solving Equations (49), we obtain the following constants as In order to separate our work from the work of [19], we refer to the third-order global Padé approximant in Equation (46) as GPadé.

Fourth-Order Global Padé Approximant
In this subsection, we present a highly accurate fourth-order Global Padé approximant for the fractional Riccati equation's solution. Setting x = 4, the fourth-order Global Padé approximant takes the form of We now present the system of equations in terms of matrix for (m, n) = (5, 4), (6, 3), (7,2). Notice that with the the combinations of m and n values all correspond to m + n = 9, which is an odd number. Before that, we define two notations that will be used to simplify the matrices in Sections 5.4.1-5.4.3: and where n ∈ N 0 . We then utilize the Equations (40) and (41)  Finally, we are able to solve for the terms p 1 , p 2 , p 3 , q 0 , q 1 , q 2 , q 3 by inversing the matrices.

Numerical Experiments
In this section, we focus on the numerical performance of the previous work by [19] (GPadé), our fourth-order global Padé approximation and the fractional Adams method introduced in Section 4.1. First, we compare the performance of the imaginary part of the solution from the fractional Riccati equation in terms of D α h(3 − i/2, t) as most of the imaginary part poses the most noticeable differences [19]. Accordingly, we will move on to the comparison of option prices obtained using GPadé [19], the fractional Adams method [61] and our newly proposed method-fourth order global Padé approximants. We will emphasize the errors of implied volatilities as a meaningful comparison.

Comparisons of Methods on the Fractional Riccati Equation
We consider the following parameters for our numerical experiments: where ν, ρ are from Equations (1) and (2), N is the number of time steps for the fractional Adams method and T is the evaluated time period. Specifically, we let h(a, t) to be the solution obtained by the fractional Adams method as it converges when N → ∞, as indicated from Equations (23) and (24). Furthermore, we denote the maximum absolute error: where H is the Hurst parameter, a = 3 − i/2 and (·) is the imaginary part of the solution. Note that we have included GPadé's solution for h (m,n) (a, t) to be considered in Equation (58). The absolute error is denoted as for the purpose of making comparisons on the solution of the fractional Riccati equation using different methods. Figure 1 shows the comparison of the solutions (D α h(3 − i/2, t)) by using the fractional Adams method, the GPadé method and the fourth-order global Padé approximation methods for (m, n) = (5, 4), (6, 3), (7, 2) on different values of the Hurst parameter H. Evidently, we can directly observe from Figures 1 and 2 and Table 1 that fourth-order global Padé approximants h (m,n) (a, t) for (m, n) = (6, 3), (7, 2) approximate the solution of the fractional Riccati equation in an extremely accurate manner. The key takeaway is that as compared to the GPadé method, the fourth-order global Padé approximants h (m,n) (a, t) for (m, n) = (6, 3), (7, 2) remain well-behaved even when the Hurst parameter approaches 0.5. GPadé [19] h (5,4) (a, t) h (6,3) (a, t) h (7,2)    Plots of E (t) for third-order GPadé approximant (dotted dark green) [19] and fourth-order h (m,n) (a, t) for (m, n) = (5, 4), (6, 3), (7, 2) (dashed-dotted dark blue, dotted mustard and dotted dark purple lines respectively) and a = 3 − i/2.

Option Pricing
In this section, we utilize the methods we have developed throughout the paper to compute the option prices under the rough Heston model. The parameters used in the option pricing are: where ν, ρ are from Equation (2), N is the number of time steps in fractional Adams method and ξ 0 (t) is the flat forward variance curve. We denote the average absolute error of the implied volatility as where k i ∈ [−0.4, 0.4], y is the number of times implied volatility is sampled, σ(T, H, k i ) is the implied volatility at log-moneyness k i obtained using the fractional Adams method and σ (m,n) (T, H, k i ) is the implied volatility at log-moneyness k i being approximated using the GPadé method and fourth-order Padé approximants. From Table 2 and Figure 3, it is obvious that the fourth-order global Padé approximants performed better in all the circumstances that we evaluated. Notably, the use of the third-order global Padé approximant in computing option prices by [19] deteriorates the quality of the solution as the Hurst parameter H approaches 0.5, whereas fourth-order global Padé approximants in computing option prices performed rather well and incurred lower errors. Fourth-order Padé approximants of h (6,3) (a, t) and h (7,2) (a, t) seem to provide the best overall performance as compared to h (5,4) (a, t) in computing the fractional Riccati equation's solution.
It is important to realize that in the setting of option pricing under the rough Heston model, the application of the fractional Adams method to the fractional Riccati equation is undesirable in the real life calibration and computation due to its excessive computational cost needed (total time complexity of O(n a N 2 ) where n a is the number of space steps for inversion of characteristic function (5), and N is the number of time steps for the fractional Adams method, (17)- (22)). For this reason, we have proposed and demonstrated the use of fourth-order global Padé approximants to approximate the fractional Riccati equation's solution, and directly applied it to compute option prices under the rough Heston model with a total time complexity of just O(n a ). By considering efficiency, accuracy and robustness (Hurst parameter H), the proposed fourth-order global Padé approximation method is definitely preferable and superior in real life practice as compared to the existing work of the fractional Adams method and the third-order global Padé approximation [19]. Table 2. Average absolute error of the implied volatility in terms of 10 −3 .

Summary and Future Research
First off, we have provided an extensive literature review relevant to the fractional Riccati equation and rough Heston model. That included the general stochastic volatility models, rough volatility, the Hurst parameter and methods for obtaining option prices from the rough Heston model. This study heavily focuses on the numerical methods for computing the fractional Riccati equation's solution. One of the key contributions from this paper is that the infinite series expansion of the fractional Riccati equation obtained using the Laplace-Adomian-decomposition method corresponds to the small time expansion of the fractional Riccati equation by [54] up to at least fifth order of the solution (t 5α ). Here we have implemented a fourth-order global Padé approximation method for the fractional Riccati equation; it includes the systems of equations in terms of the matrix to be solved in order to obtain the coefficients for fourth-order Padé approximants. The pointwise approximation errors on the global Padé approximation method for t → 0 and t → ∞ are also included, whereby it provides a sense of the number of terms (m and n) that we should pick from the small time and the large time series expansion solution of the fractional Riccati equation to solve for the coefficients of the fourth-order global Padé approximants. From the numerical computations of the fractional Riccati equation's solution using different methods, we have found that indeed the fourth-order global Padé approximants performed the most consistently and robustly of any Hurst parameter on the solution of the fractional Riccati equation, unlike the third-order Padé approximant, for which it yields erroneous solutions as the Hurst parameters get closer to 0.5. Together with the application of the fractional Riccati equation, we computed option prices under different methods, and once again the fourth-order global Padé approximants performed the best as compared to the third-order global Padé approximant from the previous study by [19]. In addition to the accurate computation of the fractional Riccati equation's solution, we note that the fourth-order global Padé approximation method possesses a time complexity of only O(1) (same as the third-order global Padé approximant) in solving the fractional Riccati equation as compared to the fractional Adams method that possesses a time complexity of O(N 2 ) where N is the number of time steps. A key limitation of this global Padé approximation method is that while it guarantees pointwise convergence (t → 0 and t → ∞) on the fractional Riccati equation's solution, it does not provide any sense of convergence at elsewhere besides the two points.
Although we are able to reduce the errors of the fractional Riccati equation's solution significantly as compared to the previous third-order Padé approximant, we wish to make three recommendations for future work: (1) consider higher order global Padé approximants such as fifth or sixth-order global Padé approximants to obtain for the fractional Riccati equation's solution and subsequently apply them to the computation of option prices under the rough Heston model; (2) it would be interesting to see that this fourth order global Padé approximant is being applied to a general setting of the fractional Riccati equation, since they are applicable in many fields, such as physics, engineering and general mathematics; (3) while the third-order and fourth-order global Padé approximants have the same time complexity O(1) to compute for the fractional Riccati equation's solution, we can further consider that by how much additional computational time is needed to switch from the third-order global Padé approximant to a fourth-order or perhaps fifth-order global Padé approximant given the advantage of accuracy and robustness in the Hurst parameter. Overall, it can be concluded that together with the small time and large time series expansions of the fractional Riccati equation's solution, fourth-order global Padé approximants can be used to efficiently and accurately compute the option prices under varying Hurst parameter H. where A n := A n (h 0 , h 1 , ..., h n ) are the Adomian polynomials and can be computed using methods stated from [59,60]. The Adomian polynomials can be computed by the formula As for f (h) := h 2 (t), the first 12 Adomian polynomials A 0 , A 1 , ..., A 11 are computed by [31] as (A9) Step 3: By substituting Equations (A6) and (A7) from Step 2 into Equation (A5), we obtain the following: By matching both sides of Equation (A10) as in the usual Adomian decomposition method, we get . . .

(A16)
By performing an inverse Laplace transform on Equation (A12) and utilizing Lemma 2, we obtain h 1 as follows.