Multi-factor approximation of rough volatility models

Rough volatility models are very appealing because of their remarkable fit of both historical and implied volatilities. However, due to the non-Markovian and non-semimartingale nature of the volatility process, there is no simple way to simulate efficiently such models, which makes risk management of derivatives an intricate task. In this paper, we design tractable multi-factor stochastic volatility models approximating rough volatility models and enjoying a Markovian structure. Furthermore, we apply our procedure to the specific case of the rough Heston model. This in turn enables us to derive a numerical method for solving fractional Riccati equations appearing in the characteristic function of the log-price in this setting.


Introduction
Empirical studies of a very wide range of assets volatility time-series in [15] have shown that the dynamics of the log-volatility are close to that of a fractional Brownian motion W H with a small Hurst parameter H of order 0.1. Recall that a fractional Brownian motion W H can be built from a two-sided Brownian motion thanks to the Mandelbrot-van Ness representation The fractional kernel (t − s) H− 1 2 is behind the H − ε Hölder regularity of the volatility for any ε > 0. For small values of the Hurst parameter H, as observed empirically, stochastic volatility models involving the fractional kernel are called rough volatility models.
Aside from modeling historical volatility dynamics, rough volatility models reproduce accurately with very few parameters the behavior of the implied volatility surface, see [3,10], especially the at-the-money skew, see [14]. Moreover, microstructural foundations of rough volatility are studied in [12,19].
In this paper, we are interested in a class of rough volatility models where the dynamics of the asset price S and its stochastic variance V are given by for all t ∈ [0, T ], on some filtered probability space (Ω, F, F, P). Here T is a positive time horizon, the parameters λ and V 0 are non-negative, H ∈ (0, 1/2) is the Hurst parameter, σ is a continuous function and W = ρB + 1 − ρ 2 B ⊥ with (B, B ⊥ ) a two-dimensional F-Brownian motion and ρ ∈ [−1, 1]. Moreover, θ is a deterministic mean reversion level allowed to be time-dependent to fit the market forward variance curve (E[V t ]) t≤T as explained in Section 2 and in [13]. Under some general assumptions, we establish in Section 2 the existence of a weak non-negative solution to the fractional stochastic integral equation in (1.2) exhibiting H − ε Hölder regularity for any ε > 0. Hence, this class of models is a natural rough extension of classical stochastic volatility models where the fractional kernel is introduced in the drift and stochastic part of the variance process V . Indeed, when H = 1/2, we recover classical stochastic volatility models where the variance process is a standard diffusion.
Despite the fit to the historical and implied volatility, some difficulties are encountered in practice for the simulation of rough volatility models and for pricing and hedging derivatives with them. In fact, due to the introduction of the fractional kernel, we lose the Markovian and semimartingale structure. In order to overcome theses difficulties, we approximate these models by simpler ones that we can use in practice.
In [11,12,13], the rough Heston model (which corresponds to the case of σ(x) = ν √ x) is built as a limit of microscopic Hawkes-based price models. This allowed the understanding of the microstructural foundations of rough volatility and also led to the formula of the characteristic function of the log-price. Hence, the Hawkes approximation enabled us to solve the pricing and hedging under the rough Heston model. However, this approach is specific to the rough Heston case and can not be extended to an arbitrary rough volatility model of the form (1.1)-(1.2).
Inspired by the works of [1,4,5,17,21], we provide a natural Markovian approximation for the class of rough volatility models (1.1)-(1.2). The main idea is to write the fractional kernel We then approximate µ by a finite sum of Dirac measures µ n = n i=1 c n i δ γ n i with positive weights (c n i ) 1≤i≤n and mean reversions (γ n i ) 1≤i≤n , for n ≥ 1. This in turn yields an approxi-mation of the fractional kernel by a sequence of smoothed kernels (K n ) n≥1 given by c n i e −γ n i t , n ≥ 1.
This leads to a multi-factor stochastic volatility model (S n , V n ) = (S n t , V n t ) t≤T , which is Markovian with respect to the spot price and n variance factors (V n,i ) 1≤i≤n and is defined as follows where dV n,i t = (−γ n i V n,i t − λV n t )dt + σ(V n t )dB t , and g n (t) = V 0 + t 0 K n (t − s)θ(s)ds with the initial conditions S n 0 = S 0 and V n,i 0 = 0. Note that the factors (V n,i ) 1≤i≤n share the same dynamics except that they mean revert at different speeds (γ n i ) 1≤i≤n . Relying on existence results of stochastic Volterra equations in [1,2], we provide in Theorem 3.1 the strong existence and uniqueness of the model (S n , V n ), under some general conditions. Thus the approximation (1.4) is uniquely well-defined. We can therefore deal with simulation, pricing and hedging problems under these multi-factor models by using standard methods developed for stochastic volatility models. Theorem 3.5, which is the main result of this paper, establishes the convergence of the multifactor approximation sequence (S n , V n ) n≥1 to the rough volatility model (S, V ) in (1.1)-(1.2) when the number of factors n goes to infinity, under a suitable choice of the weights and mean reversions (c n i , γ n i ) 1≤i≤n . This convergence is obtained from a general result about stability of stochastic Volterra equations derived in Section 3.4.
In [2,11,13], the characteristic function of the log-price for the specific case of the rough Heston model is obtained in terms of a solution of a fractional Riccati equation. We highlight in Section 4.1 that the corresponding multi-factor approximation (1.4) inherits a similar affine structure as in the rough Heston model. More precisely, it displays the same characteristic function formula involving a n-dimensional classical Riccati equation instead of the fractional one. This suggests solving numerically the fractional Riccati equation by approximating it through a n-dimensional classical Riccati equation with large n, see Theorem 4.1. In Section 4.2, we discuss the accuracy and complexity of this numerical method and compare it to the Adams scheme, see [7,8,9,11].
The paper is organized as follows. In Section 2, we define the class of rough volatility models (1.1)-(1.2) and discuss the existence of such models. Then, in Section 3, we build a sequence of multi-factor stochastic volatility models of the form of (1.4) and show its convergence to a rough volatility model. By applying this approximation to the specific case of the rough Heston model, we obtain a numerical method for computing solutions of fractional Riccati equations that is discussed in Section 4. Finally, some proofs are relegated to Section 5 and some useful technical results are given in an Appendix.

A definition of rough volatility models
We provide in this section the precise definition of rough volatility models given by (1.1)-(1.2). We discuss the existence of such models and more precisely of a non-negative solution of the fractional stochastic integral equation (1.2). The existence of an unconstrained weak solution V = (V t ) t≤T is guaranteed by Corollary B.2 in the Appendix when σ is a continuous function with linear growth and θ satisfies the condition Furthermore, the paths of V are Hölder continuous of any order strictly less than H and Moreover using Theorem B.4 together with Remarks B.5 and B.6 in the Appendix 1 , the existence of a non-negative continuous process V satisfying (1.2) is obtained under the additional conditions of non-negativity of V 0 and θ and σ(0) = 0. We can therefore introduce the following class of rough volatility models.
As done in [13], we allow the mean reversion level θ to be time dependent in order to be consistent with the market forward variance curve. More precisely, the following result shows that the mean reversion level θ can be written as a functional of the forward variance curve Proposition 2.2. Let (S, V ) be a rough volatility model given by Definition 2.1. Then, t≤T is linked to θ by the following formula

3)
where α = H + 1/2 and E α (x) = k≥0 x k Γ(α(k+1)) is the Mittag-Leffler function. Moreover, (E[V t ]) t≤T admits a fractional derivative 2 of order α at each time t ∈ (0, T ] and 2 Recall that the fractional derivative of order α ∈ (0, 1) of a function f is given by Proof. Thanks to (2.2) together with Fubini theorem, t → E[V t ] solves the following fractional linear integral equation

Multi-factor approximation of rough volatility models
Thanks to the small Hölder regularity of the variance process, models of Definition 2.1 are able to reproduce the rough behavior of the volatility observed in a wide range of assets. However, the fractional kernel forces the variance process to leave both the semimartingale and Markovian worlds, which makes numerical approximation procedures a difficult and challenging task in practice. The aim of this section is to construct a tractable and satisfactory Markovian approximation of any rough volatility model (S, V ) of Definition 2.1. Because S is entirely determined by ( , it suffices to construct a suitable approximation of the variance process V . This is done by smoothing the fractional kernel. More precisely, denoting by K(t) = t H− 1 2 Γ(H+1/2) , the fractional stochastic integral equation (1.2) reads which is a stochastic Volterra equation. Approximating the fractional kernel K by a sequence of smooth kernels (K n ) n≥1 , one would expect the convergence of the following corresponding sequence of stochastic Volterra equations to the fractional one.
The argument of this section runs as follows. First, exploiting the identity (1.3), we construct a family of potential candidates for (K n , V n ) n≥1 in Section 3.1 such that V n enjoys a Markovian structure. Second, we provide convergence conditions of (K n ) n≥1 to K in L 2 ([0, T ], R) in Section 3.2. Finally, the approximation result for the rough volatility model (S, V ) is established in Section 3.3 relying on an abstract stability result of stochastic Volterra equations postponed to Section 3.4 for sake of exposition.

Construction of the approximation
In [4,17,21], a Markovian representation of the fractional Brownian motion of Riemann-Liouville type is provided by writing the fractional kernel K(t) = t H− 1 2 Γ(H+1/2) as a Laplace transform of a non-negative measure µ as in (1.3). This representation is extended in [1] for the Volterra square-root process. Adopting the same approach, we establish a similar representation for any solution of the fractional stochastic integral equation (1.2) in terms of an infinite dimensional system of processes sharing the same Brownian motion and mean reverting at different speeds. Indeed by using the linear growth of σ together with the stochastic Fubini theorem, see [24], we obtain that Inspired by [4,5], we approximate the measure µ by a weighted sum of Dirac measures leading to the following approximation V n = (V n t ) t≤T of the variance process V and The choice of the positive weights (c n i ) 1≤i≤n and mean reversions (γ n i ) 1≤i≤n , which is crucial for the accuracy of the approximation, is studied in Section 3.2 below. Before proving the convergence of (V n ) n≥1 , we shall first discuss the existence and uniqueness of such processes. This is done by rewriting the stochastic equation (3.2) as a stochastic Volterra equation of the form The existence of a continuous non-negative weak solution V n is ensured by Theorem B.4 together with Remarks B.5 and B.6 in the Appendix 3 , because θ and V 0 are non-negative and σ(0) = 0. Moreover, pathwise uniqueness of solutions to (3.5) follows by adapting the standard arugments of [25], provided a suitable Hölder continuity of σ, see Proposition B.3 in the Appendix. Note that this extension is made possible due to the smoothness of the kernel K n . For instance, this approach fails for the fractional kernel because of the singularity at zero. This leads us to the following result which establishes the strong existence and uniqueness of a non-negative solution of (3.5) and equivalently of (3.2).
Theorem 3.1. Assume that θ : [0, T ] → R is a deterministic non-negative function satisfying (2.1) and that σ : R → R is η-Hölder continuous with σ(0) = 0 and η ∈ [1/2, 1]. Then, there exists a unique strong non-negative solution V n = (V n t ) t≤T to the stochastic Volterra equation Due to the uniqueness of (3.2), we obtain that V n is a Markovian process according to n state variables (V n,i ) 1≤i≤n that we call the factors of V n . Moreover, V n being non-negative, it can model a variance process. This leads to the following definition of multi-factor stochastic volatility models.
Definition 3.2. (Multi-factor stochastic volatility models). We define the following sequence of multi-factor stochastic volatility models (S n , V n ) = (S n t , V n t ) t≤T as the unique R × R +valued strong solution of with dV n,i t = (−γ n i V n,i t − λV n t )dt + σ(V n t )dB t , V n,i 0 = 0, S n 0 = S 0 > 0, on a filtered probability space (Ω, F, P, F), where F is the canonical filtration a two-dimensional Brownian motion (W, W ⊥ ) and B = ρW + 1 − ρ 2 W ⊥ with ρ ∈ [−1, 1]. Here, the weights (c n i ) 1≤i≤n and mean reversions (γ n i ) 1≤i≤n are positive, σ : R → R is η-Hölder continuous such that σ(0) = 0, η ∈ [1/2, 1] and g n is given by (3.3), that is with a non-negative initial variance V 0 , a kernel K n defined as in (3.4) and a non-negative deterministic function θ : [0, T ] → R satisfying (2.1).
Note that the strong existence and uniqueness of (S n , V n ) follows from Theorem 3.1. This model is Markovian with n + 1 state variables which are the spot price S n and the factors of the variance process V n,i for i ∈ {1, . . . , n}.

An approximation of the fractional kernel
Relying on (3.5), we can see the process V n as an approximation of V , solution of (1.2), obtained by smoothing the fractional kernel K(t) = t H− 1 2 Γ(H+1/2) into K n (t) = n i=1 c n i e −γ n i t . Intuitively, we need to choose K n close to K when n goes to infinity, so that (V n ) n≥1 converges to V . Inspired by [5], we give in this section a condition on the weights (c n i ) 1≤i≤n and mean reversion terms 0 < γ n 1 < ... < γ n n so that the following convergence holds as n goes to infinity, where · 2,T is the usual L 2 ([0, T ], R) norm. Let (η n i ) 0≤i≤n be auxiliary mean reversion terms such that η n 0 = 0 and η n i−1 ≤ γ n i ≤ η n i for i ∈ {1, . . . , n}. Writing K as the Laplace transform of µ as in (1.3), we obtain that We start by dealing with the first term,

Moreover by choosing
and using the Taylor-Lagrange inequality up to the second order, we obtain This leads to the following inequality n is a function of the auxiliary mean reversions defined by Hence, we obtain the convergence of K n to the fractional kernel under the following choice of weights and mean reversions.
Assumption 3.1. We assume that the weights and mean reversions are given by (3.6) such that η n 0 = 0 < η n 1 < . . . < η n n and as n goes to infinity. There exists several choices of auxiliary factors such that condition (3.9) is met. For instance, assume that η n i = iπ n for each i ∈ {0, . . . , n} such that π n > 0. It follows from as n tends to infinity. In this case, This upper bound is minimal for and where C H is a positive constant that can be computed explicitly and that depends only on the Hurst parameter H ∈ (0, 1/2). Remark 3.4. Note that the kernel approximation in Proposition 3.3 can be easily extended to any kernel of the form where µ is a non-negative measure such that

Convergence result
We assume now that the weights and mean reversions of the multi-factor stochastic volatility model (S n , V n ) satisfy Assumption 3.1. Thanks to Proposition 3.3, the smoothed kernel K n is close to the fractional one for large n. Because V n satisfies the stochastic Volterra equation (3.5), V n has to be close to V and thus by passing to the limit, (S n , V n ) n≥1 should converge to the rough volatility model (S, V ) of Definition 2.1 as n goes large. This is the object of the next theorem, which is the main result of this paper.
Theorem 3.5. Let (S n , V n ) n≥1 be a sequence of multi-factor stochastic volatility models given by Definition 3.2. Then, under Assumption 3.1, the family (S n , V n ) n≥1 is tight for the uniform topology and any point limit (S, V ) is a rough volatility model given by Definition 2.1.
Theorem 3.5 states the convergence in law of (S n , V n ) n≥1 whenever the fractional stochastic integral equation (1.2) admits a unique weak solution. In order to prove Theorem 3.5, whose proof is in Section 5.2 below, a more general stability result for d-dimensional stochastic Volterra equations is established in the next subsection.

Stability of stochastic Volterra equations
As mentioned above, Theorem 3.5 relies on the study of the stability of more general ddimensional stochastic Volterra equations of the form for any t, h ≥ 0 with t + h ≤ T and for some positive constants C and γ, to guarantee the weak existence of a continuous solution X of (3.11).
More precisely, we consider a sequence X n = (X n t ) t≤T of continuous weak solutions to the stochastic Volterra equation (3.11) with a kernel K n ∈ L 2 ([0, T ], R d×d ) admitting a resolvent of the first kind, on some filtered probability space (Ω n , F n , F n , P n ), with g n : [0, T ] → R d and K n satisfying (3.12) for every n ≥ 1. The stability of (3.11) means the convergence in law of the family of solutions (X n ) n≥1 to a limiting process X which is a solution to (3.11), when (K n , g n ) is close to (K, g) as n goes large.
This convergence is established by verifying first the Kolmogorov tightness criterion for the sequence (X n ) n≥1 . It is obtained when g n and K n satisfy (3.12) uniformly in n in the following sense.
Assumption 3.2. There exists positive constants γ and C such that The following result, whose proof is postponed to Section 5.1 below, states the convergence of (X n ) n≥1 to a solution of (3.11).
for any t ∈ [0, T ] as n goes to infinity. Then, under Assumption 3.2, the sequence (X n ) n≥1 is tight for the uniform topology and any point limit X is a solution of the stochastic Volterra equation (3.11).

The particular case of the rough Heston model
The rough Heston model introduced in [11,13] is a particular case of the class of rough volatility models of Definition 2.1, with σ(x) = ν √ x for some positive parameter ν, that is Γ(H+1/2) denotes the fractional kernel and g is given by (3.1). Aside from reproducing accurately the historical and implied volatility, the rough Heston model displays a closed formula for the characteristic function of the log-price in terms of a solution to a fractional Riccati equation allowing to fast pricing and calibration, see [10]. More precisely, it is shown in [1,11,13] that where ψ(·, z) is the unique continuous solution of the fractional Riccati equation . Unlike the classical case H = 1/2, (4.2) does not exhibit an explicit solution. However, it can be solved numerically through the Adam scheme developed in [7,8,9,11] for instance. In this section, we show that the multi-factor approximation applied to the rough Heston model gives rise to another natural numerical scheme for solving the fractional Riccati equation. Furthermore, we will establish the convergence of this scheme with explicit errors.

Multi-factor scheme for the fractional Riccati equation
We consider the multi-factor approximation (S n , V n ) of Definition 3.2 with σ(x) = ν √ x, where the number of factors n is large, that is Recall that g n is given by (3.3) and it converges pointwise to g as n goes large, see Lemma 5.1.
We write the dynamics of (S n , V n ) in terms of a Volterra Heston model with the smoothed kernel K n given by (3.4) as follows In [1,2], the characteristic function formula of the log-price (4.1) is extended to the general class of Volterra Heston models. In particular, is given by where ψ n (·, z) is the unique continuous solution of the Riccati Volterra equation Thanks to the weak uniqueness of the rough Heston model, established in several works [1,2,22], and to Theorem 3.5, (S n , V n ) n≥1 converges in law for the uniform topology to (S, V ) when n tends to infinity. In particular, L n (t, z) converges pointwise to L(t, z). Therefore, we expect ψ n (·, z) to be close to the solution of the fractional Riccati equation (4.2). This is the object of the next theorem, whose proof is reported to Section 5.3 below.
Relying on the L 1 -convergence of (K n ) n≥1 to K under Assumption 3.1, see Proposition 3.3, we have the uniform convergence of (ψ n (·, z)) n≥1 to ψ(·, z) on [0, T ]. Hence, Theorem 4.1 suggests a new numerical method for the computation of the fractional Riccati solution (4.2) where an explicit error is given. Indeed, set Then, and (ψ n,i (·, z)) 1≤i≤n solves the following n-dimensional system of ordinary Riccati equations Hence, (4.5) can be solved numerically by usual finite difference methods leading to ψ n (·, z) as an approximation of the fractional Riccati solution.

Numerical illustrations
In this section, we consider a rough Heston model with the following parameters We discuss the accuracy of the multi-factor approximation sequence (S n , V n ) n≥1 as well as the corresponding Riccati Volterra solution (ψ n (·, z)) n≥1 , for different choices of the weights (c n i ) 1≤i≤n and mean reversions (γ n i ) 1≤i≤n . This is achieved by first computing, for different number of factors n, the implied volatility σ n (k, T ) of maturity T and log-moneyness k by a Fourier inversion of the characteristic function formula (4.3), see [6,20] for instance. In a second step, we compare σ n (k, T ) to the implied volatility σ(k, T ) of the rough Heston model. We also compare the Riccati Volterra solution ψ n (T, z) to the fractional one ψ(T, z).
Note that the Riccati Volterra solution ψ n (·, z) is computed by solving numerically the ndimensional Riccati equation (4.5) with a classical finite difference scheme. The complexity of such scheme is O(n × n ∆t ), where n ∆t is the number of time steps applied for the scheme, while the complexity of the Adam scheme used for the computation of ψ(·, z) is O(n 2 ∆t ). In the following numerical illustrations, we fix n ∆t = 200.
In order to guarantee the convergence, the weights and mean reversions have to satisfy Assumption 3.1 and in particular they should be of the form (3.6) in terms of auxiliary mean reversions (η n i ) 0≤i≤n satisfying (3.9). For instance, one can fix where π n is defined by (3.10), as previously done in Section 3.2. For this particular choice, Figure 1 shows a decrease of the relative error ψ n (T,ib)−ψ(T,ib) towards zero for different values of b. We also observe in the Figure 2 below that the implied volatility σ n (k, T ) of the multifactor approximation is close to σ(k, T ) for a number of factors n ≥ 20. Notice that the approximation is more accurate around the money. In order to obtain a more accurate convergence, we can minimize the upper bound f (3.8). Hence, we choose (η n i ) 0≤i≤n to be a solution of the constrained minimization problem inf | is smaller under the choice of factors (4.7). Indeed the Volterra approximation ψ n (T, ib) is now closer to the fractional Riccati solution ψ(T, ib) especially for small number of factors. However, when n is large, the accuracy of the approximation seems to be close to the one under (4.6). For instance when n = 500, the relative error is around 1% under both (4.6) and (4.7). In the same way, we observe in Figure 4 that the accuracy of the implied volatility approximation σ n (k, T ) is more satisfactory under (4.7) especially for a small number of factors.
Theorem 4.1 states that the convergence of ψ n (·, z) depends actually on the L 1 ([0, T ], R)-error between K n and K. Similarly to the computations of Section 3.2, we may show that, This leads to choosing (η n i ) 0≤i≤n as a solution of the constrained minimization problem It is easy to show that such auxiliary mean-reversions (η n i ) 0≤i≤n satisfy (3.9) and thus Assumption 3.1 is met.

Upper bound for call prices error
Using a Fourier transform method, we can also provide an error between the price of the call C n (k, T ) = E[(S n T − S 0 e k ) + ] in the multi-factor model and the price of the same call C(k, T ) = E[(S T − S 0 e k ) + ] in the rough Heston model. However, for technical reasons, this bound is obtained for a modification of the multi-factor approximation (S n , V n ) n≥1 of Definition 3.2 where the function g n given initially by (3.3) is updated into where K n is the smoothed approximation (3.4) of the fractional kernel. Note that the strong existence and uniqueness of V n is still directly obtained from Proposition B.3 and its nonnegativity from Theorem B.4 together with Remarks B.5 and B.6 in the Appendix 4 . Although for g n satisfying (4.9), (V n ) n≥1 can not be tight 5 , the corresponding spot price (S n ) n≥1 converges as shown in the following proposition.
Note that the characteristic function (4.3) still holds. Using Theorem 4.1 together with a Fourier transform method, we obtain an explicit error for the call prices. We refer to Section 5.5 below for the proof.
Proposition 4.3. Let C(k, T ) be the price of the call in the rough Heston model with maturity T > 0 and log-moneyness k ∈ R. We denote by C n (k, T ) the price of the call in the multifactor Heston model of Definition 3.2 such that g n is given by (4.9). If |ρ| < 1, then there exists a positive constant c > 0 such that

Proofs
In this section, we use the convolution notations together with the resolvent definitions of Appendix A. We denote by c any positive constant independent of the variables t, h and n and that may vary from line to line. For any h ∈ R, we will use the notation ∆ h to denote the semigroup operator of right shifts defined by ∆ h f : t → f (h + t) for any function f .
We first prove Theorem 3.6, which is the building block of Theorem 3.5. Then, we turn to the proofs of the results contained in Section 4, which concern the particular case of the rough Heston model.

Proof of Theorem 3.6
Tightness of (X n ) n≥1 : We first show that, for any p ≥ 2, Thanks to Proposition B.1, we already have Using the linear growth of (b, σ) and (5.2) together with Jensen and BDG inequalities, we get Relying on Assumption 3.2 and the convergence of (g n (0), T 0 |K n (s)| 2 ds) n≥1 , sup t≤T |g n (t)| p and T 0 |K n (s)| 2 ds are uniformly bounded in n. This leads to By the Grönwall type inequality in Lemma A.4 in the Appendix, we deduce that where E n c ∈ L 1 ([0, T ], R) is the canonical resolvent of |K n | 2 with parameter c, defined in Appendix A.3, and the last inequality follows from the fact that We now show that (X n ) n≥1 exhibits the Kolmogorov tightness criterion. In fact, using again the linear growth of (b, σ) and (5.1) together with Jensen and BDG inequalities, we obtain, for any p ≥ 2 and t, h ≥ 0 such that t + h ≤ T , Hence, Assumption 3.2 leads to and therefore to the tightness of (X n ) n≥1 for the uniform topology.
Convergence of (X n ) n≥1 : Let M n t = t 0 σ(X n s )dW n s . As M n t = t 0 σσ * (X n s )ds, ( M n ) n≥1 is tight and consequently we get the tightness of (M n ) n≥1 from [18, Theorem VI-4.13]. Let (X, M ) = (X t , M t ) t≤T be a possible limit point of (X n , M n ) n≥1 . Thanks to [18, Theorem VI-6.26], M is a local martingale and necessarily Moreover, setting Y n t = t 0 b(X n s )ds + M n t , the assoicativity property (A.1) in the Appendix yields (L * X n ) t = (L * g n )(t) + L * (K n − K) * dY n t + Y n t , where L is the resolvent of the first kind of K defined in Appendix A.2. By the Skorokhod representation theorem, we construct a probability space supporting a sequence of copies of (X n , M n ) n≥1 that converges uniformly on [0, T ], along a subsequence, to a copy of (X, M ) almost surely, as n goes to infinity. We maintain the same notations for these copies. Hence, we have sup almost surely, as n goes to infinity. Relying on the continuity and linear growth of b together with the dominated convergence theorem, it is easy to obtain for any t ∈ [0, T ] almost surely as n goes to infinity. Moreover for each t ∈ [0, T ] by the uniform boundedness of g n in n and t and the dominated convergence theorem. Finally thanks to the Jensen inequality, From (5.1) and the linear growth of (b, σ), we deduce which goes to zero when n is large. Consequently, we send n to infinity in (5.3) and obtain the following almost surely equality, for each t ∈ [0, T ], Recall also that M = · 0 σσ * (X s )ds. Hence, by [23, Theorem V-3.9], there exists a mdimensional Brownian motion W such that The processes in (5.4) being continuous, we deduce that, almost surely, We convolve by K and use the associativity property (A.1) in the Appendix to get that, almost surely, Finally it is easy to see that the processes above are differentiable and we conclude that X is solution of the stochastic Volterra equation (3.11), by taking the derivative.

Proof of Theorem 3.5
Theorem 3.5 is easily obtained once we prove the tightness of (V n ) n≥1 for the uniform topology and that any limit point V is solution of the fractional stochastic integral equation (1.2). This is a direct consequence of Theorem 3.6, by setting d = m = 1, g and g n respectively as in (3.1) and (3.3), b(x) = −λx, K being the fractional kernel and K n (t) = n i=1 c n i e −γ n i t its smoothed approximation. Under Assumption 3.1, (K n ) n≥1 converges in L 2 ([0, T ], R) to the fractional kernel, see Proposition 3.3. Hence, it is left to show the pointwise convergence of (g n ) n≥1 to g on [0, T ] and that (K n , g n ) n≥1 satisfies Assumption 3.2. converges to zero as n goes large, for some ε > 0 and K n given by (3.4). Using the representation of K as the Laplace transform of µ as in (1.3), we obtain that (5.5) is bounded by The first term in (5.6) converges to zero for large n by the dominated convergence theorem because η n n tends to infinity, see Assumption 3.1. Using the Taylor-Lagrange inequality (3.7), the second term in (5.6) is dominated by which goes to zero thanks to Assumption 3.1.

Lemma 5.2 (K n satisfying Assumption 3.2). Under Assumption 3.1, there exists
where K n is defined by (3.4).

Proof. We start by proving that for any
In fact we know that this inequality is satisfied for K(t) = t H− 1 2 Γ(H+1/2) . Thus it is enough to prove where · 2,h stands for the usual L 2 ([0, h], R) norm. Relying on the Laplace transform representation of K given by (1.3), we obtain We start by bounding the first term, As in Section 3.2, we use the Taylor-Lagrange inequality (3.7) to get Using the boundedness of n i=1 η n i η n i−1 (γ − γ n i ) 2 µ(dγ) n≥1 from Assumption 3.1, we deduce (5.7). We now prove In the same way, it is enough to show Similarly to the previous computations, we get Moreover, fix h, t > 0 and set χ(γ) = e −γt − e −γ(t+h) . The second derivative reads The Taylor-Lagrange formula, up to the second order, leads to   Proof. Because θ satisfies (2.1), it is enough to prove that, for each fixed ε > 0, there exists C > 0 such that for any t, h ≥ 0 with t + h ≤ T . (5.10) being satisfied for the fractional kernel, it is enough to establish In the proof of Lemma 5.1, it is shown that

The first term is dominated by
where B is the usual Beta function. Moreover thanks to (3.7) and Assumption 3.1, we get By similar computations as previously and using (5.9), we get that The first term being bounded by Assumption 3.1 leads to (5.11).

Proof of Theorem 4.1
Uniform boundedness : We start by showing the uniform boundedness of the unique continuous solutions (ψ n (·, a + ib)) n≥1 of (4.4).
Proposition 5.4. For a fixed T > 0, there exists C > 0 such that for any a ∈ [0, 1] and b ∈ R.
Proof. Let z = a + ib and start by noticing that (ψ n (·, z)) is non-positive because it solves the following linear Volterra equation with continuous coefficients is continuous non-positive, see Theorem C.1. In the same way (ψ(·, z)) is also non-positive. Moreover, observe that ψ n (·, z) solves the following linear Volterra equation with continuous coefficients Therefore, Corollary C.4 leads to where E n ν−λ denotes the canonical resolvent of K n with parameter ν − λ, see Appendix A.3. This resolvent converges in L 1 ([0, T ], R) because K n converges in L 1 ([0, T ], R) to K, see [16,Theorem 2.3.1]. Hence, ( T 0 E n ν−λ (s)ds) n≥1 is bounded, which ends the proof.

Proof of Proposition 4.2
We consider for each n ≥ 1, (S n , V n ) defined by the multi-factor Heston model in Definition Tightness of ( · 0 V n s ds, Thus t → E[V n t ] solves the following linear Volterra equation with K n given by (3.4). Theorem A.3 in the Appendix leads to and then by Fubini theorem again where E n λ is the canonical resolvent of K n with parameter λ, defined in Appendix A.3. Because (K n ) n≥1 converges to the fractional kernel K in L 1 ([0, T ], R), we obtain the convergence of E n λ in L 1 ([0, T ], R) to the canonical resolvent of K with parameter λ, see [16,Theorem 2.3.1]. In particular thanks to Corollary C.2 in the Appendix, t 0 E n λ (s)ds is uniformly bounded in t ∈ [0, T ] and n ≥ 1. This leads to (5.13) and then to the tightness of ( · 0 V n s ds, V n s dB s ) n≥1 by [18,.
Convergence of (S n , V n s dB s . Denote by (X, M 1 , M 2 ) a limit in law for the uniform topology of a subsequence of the tight family ( · 0 V n s ds, M n,1 , M n,2 ) n≥1 . An application of stochastic Fubini theorem, see [24], yields which converges in law for the uniform topology to zero thanks to the convergence of ( · 0 V n s ds) n≥1 together with Proposition 3.3. Finally, which goes to zero thanks to (5.13) and Proposition 3.3. Hence, by passing to the limit in (5.14), we obtain for any t ∈ [0, T ], almost surely. The processes being continuous, the equality holds on [0, T ]. Then, by the stochastic Fubini theorem, we deduce that X = · 0 V s ds, where V is a continuous process defined by In particular V is solution of the fractional stochastic integral equation in Definition 2.1 with σ(x) = ν √ x. Because S n = exp(M n,1 − 1 2 · 0 V n s ds), we deduce the convergence of (S n , · 0 V n s ds) n≥1 to the limit point (S, · 0 V s ds) that displays the rough-Heston dynamics of Definition 2.1. The uniqueness of such dynamics, see [1,2,22], enables us to conclude that (S n , V n ) n≥1 admits a unique limit point and hence converges to the rough Heston dynamics.

Proof of Proposition 4.3
We use the Lewis Fourier inversion method, see [20], to write Hence, Because L(T, z) and L n (T, z) satisfy respectively the formulas (4.1) and (4.3) with g and g n given by and ψ(·, z) and ψ n (·, z) solve respectively (4.2) and (4.4), we use the Fubini theorem to deduce that and with z = 1/2 + ib. Therefore, relying on the local Lipschitz property of the exponential function, it suffices to find an upper bound for (ψ n (·, z)) in order to get an error for the price of the call from (5.15). This is the object of the next paragraph.
where E n η is the canonical resolvent of K n with parameter η defined in Appendix A.3.
Proof. The claimed monotonicity of b → φ n η (t 0 , b) is directly obtained from Theorem C.1. Consider now h, b 0 > 0. It is easy to see that ∆ h φ n η (·, b 0 ) solves the following Volterra equation ) being solution of the following linear Volterra integral equation with continuous coefficients, we deduce its non-negativity using again Theorem C.1. Thus, t ∈ R + → φ n η (t, b 0 ) is nonincreasing and consequently sup s∈[0,t] |φ η (s, b)| = |φ n η (t, b 0 )| as φ n η (0, b) = 0. Hence, Theorem A.3 leads to We end the proof by solving this inequality of second order in φ n η (t, b) and using that φ n η is non-positive. Notice that t 0 E n η (s)ds > 0 for each t > 0, see Corollary C.2.
Finally from Lemma 2.4 in [2] together with the Kolmogorov continuity theorem, we can show that there exists a unique version of (K * dM t ) t≥0 that is continuous whenever b and σ are locally bounded. In this paper, we will always work with this continuous version.
Note that the convolution notation could be easily extended for matrix-valued K and L. In this case, the associativity properties exposed above hold.

A.2 Resolvent of the first kind
We define the resolvent of the first kind of a d × d-matrix valued kernel K, as the R d×d -valued measure L on R + of locally bounded variation such that where id stands for the identity matrix, see [16,Definition 5.5.1]. The resolvent of the first kind does not always exist. In the case of the fractional kernel K(t) = t H− 1 2 Γ(H+1/2) the resolvent of the first kind exists and is given by for any H ∈ (0, 1/2). If K is non-negative, non-increasing and not identically equal to zero on R + , the existence of a resolvent of the first kind is guaranteed by [16,Theorem 5.5.5].
The following result shown in [2, Lemma 2.6], is stated here for d = 1 but is true for any dimension d ≥ 1.
Lemma A.1. Assume that K ∈ L 1 loc (R + , R) admits a resolvent of first kind L. For any F ∈ L 1 loc (R + , R) such that F * L is right-continuous and of locally bounded variation one has Here, df denotes the measure such that f (t) = f (0) + [0,t] df (s), for all t ≥ 0, for any right-continuous function of locally bounde variation f on R + .
Remark A.2. The previous lemma will be used with F = ∆ h K, for fixed h > 0. If K is continuous on (0, ∞), then ∆ h K * L is right-continuous. Moreover, if K is non-negative and L non-increasing in the sense that s → L([s, s + t]) is non-increasing for all t ≥ 0, then ∆ h K * L is non-decreasing since In particular, ∆ h K * L is of locally bounded variation.

A.3 Resolvent of the second kind
We consider a kernel K ∈ L 1 loc (R + , R) and define the resolvent of the second kind of K as the unique function R K ∈ L 1 loc (R + , R) such that For λ ∈ R, we define the canonical resolvent of K with parameter λ as the unique solution This means that E λ = −R −λK /λ, when λ = 0 and E 0 = K. The existence and uniqueness of R K and E λ is ensured by [16,Theorem 2.3.1] together with the continuity of K → E λ (K) in the topology of L 1 loc (R + , R). Moreover, if K ∈ L 2 loc (R + , R) so does E λ due to [16,Theorem 2.3.5].
We recall [16,Theorem 2.3.5] regarding the existence and uniqueness of a solution of linear Volterra integral equations in L 1 loc (R + , R).
Theorem A.3. Let f ∈ L 1 loc (R + , R). The integral equation admits a unique solution x ∈ L 1 loc (R + , R) given by When K and λ are positive, E λ is also positive, see [16,Proposition 9.8.1]. In that case, we have a Grönwall type inequality given by [16,Lemma 9.8.2]. Then, Note that the definition of the resolvent of the second kind and canonical resolvent can be extended for matrix-valued kernels. In that case, Theorem A.3 still holds.

B Some existence results for stochastic Volterra equations
We collect in this Appendix existence results for general stochastic Volterra equations as introduced in [2]. We refer to [1,2] for the proofs. We fix T > 0 and consider the ddimensional stochastic Volterra equation is a kernel admitting a resolvent of the first kind L, g : [0, T ] → R d is a continuous function and B is a m-dimensional Brownian motion on a filtered probability space (Ω, F, F, P). In order to prove the weak existence of continuous solutions to (B.1), the following regularity assumption is needed.
Assumption B.1. There exists γ > 0 and C > 0 such that for any t, h ≥ 0 with t + h ≤ T , The following existence result can be found in [ and admits Hölder continuous paths on [0, T ] of any order strictly less than γ.
In particular, for the fractional kernel, Proposition B.1 yields the following result.
Corollary B.2. Fix H ∈ (0, 1/2) and θ : [0, T ] → R satisfying The fractional stochastic integral equation admits a weak continuous solution X = (X t ) t≤T for any X 0 ∈ R. Moreover X satisfies (B.2) and admits Hölder continuous paths on [0, T ] of any order strictly less than H.
Proof. It is enough to notice that the fractional stochastic integral equation is a particular case of (B.1) with d = m = 1, K(t) = t H− 1 2 Γ(H+1/2) the fractional kernel, which admits a resolvent of the first kind, see Section A.2, and As t → t 1/2+ε θ(t) is bounded on [0, T ], we may show that g is H − ε Hölder continuous for any ε > 0. Hence, Assumption B.1 is satisfied and the claimed result is directly obtained from Proposition B.1.
We now establish the strong existence and uniqueness of (B.1) in the particular case of smooth kernels. This is done by extending the Yamada-Watanabe pathwise uniqueness proof in [25]. Proposition B.3. Fix m = d = 1 and assume that g is Hölder continuous, K ∈ C 1 ([0, T ], R) admitting a resolvent of the first kind and that there exists C > 0 and η ∈ [1/2, 1] such that for any x, y ∈ R, Then, the stochastic Volterra equation (B.1) admits a unique strong continuous solution.
Proof. We start by noticing that, K being smooth, it satisfies Assumption B.1. Hence, the existence of a weak continuous solution to (B.1) follows from Proposition B.1. It is therefore enough to show the pathwise uniqueness. We may proceed similarly to [25] by considering a 0 = 1, a k−1 > a k for k ≥ 1 with a k−1 a k x −2η dx = k and ϕ k ∈ C 2 (R, R) such that ϕ k (x) = ϕ k (−x), ϕ k (0) = 0 and for x > 0 • ϕ k (x) = 0 for x ≤ a k , ϕ k (x) = 1 for x ≥ a k−1 and ϕ k (x) ∈ [0, 1] for a k < x < a k−1 .
In [1], the proof of [2, Theorem 3.5] is adapted to prove the existence of a non-negative solution for a wide class of admissible input curves g satisfying 6 We therefore define the following set of admissible input curves Remark B.5. Note that any locally square-integrable completely monotone kernel 7 that is not identically zero satisfies Assumption B.2, see [2,Example 3.6]. In particular, this is the case for • the fractional kernel K(t) = t H−1/2 Γ(H+1/2) , with H ∈ (0, 1/2).
• any weighted sum of exponentials K(t) = n i=1 c i e −γ i t such that c i , γ i ≥ 0 for all i ∈ {1, . . . , n} and c i > 0 for some i.
Remark B.6. Theorem B.4 will be used with functions g of the following form where ξ is a non-negative measure of locally bounded variation and c is a non-negative constant. In that case, we may show that (B.3) is satisfied, under Assumption B.2.

C Linear Volterra equation with continuous coefficients
In this section, we consider K ∈ L 2 loc (R + , R) satisfying Assumption B.2 with T = ∞ and recall the definition of G K , that is Theorem C.1. Let K ∈ L 2 loc (R + , R) satisfying Assumption B.2 and g, z, w : R + → R be continuous functions. The linear Volterra equation admits a unique continuous solution χ. Furthermore if g ∈ G K and w is non-negative, then χ is non-negative and Proof. The existence and uniqueness of such solution in χ ∈ L 1 loc (R + , R) is obtained from [2, Lemma C.1]. Because χ is solution of (C.1), it is enough to show the local boundedness of χ to get its continuity. This follows from Grönwall's Lemma A.4 applied on the following inequality for any t ∈ [0, T ] and for a fixed T > 0.
We assume now that g ∈ G K and w is non-negative. The fact that g t 0 ∈ G K , for t 0 ≥ 0, is proved by adapting the computations of the proof of [1, Theorem 3.1] with ν = 0 provided that χ is non-negative. In order to establish the non-negativity of χ, we introduce, for each ε > 0, χ ε as the unique continuous solution of It is enough to prove that χ ε is non-negative, for every ε > 0, and that (χ ε ) ε>0 converges uniformly on every compact to χ as ε goes to zero.
We now provide a version of Theorem C.1 for complex valued solutions. Proof. The existence and uniqueness of a continuous solution is obtained in the same way as in the proof of Theorem C.1. Consider now, for each ε > 0, ψ ε the unique continuous solution of ψ ε = |h 0 | + K * ( (z)ψ + |w| + ε).
As done in the proof of Theorem C.1, ψ ε converges uniformly on every compact to ψ as ε goes to zero. Thus, it is enough to show that, for every ε > 0 and t ≥ 0, |h(t)| ≤ ψ ε (t).
We start by showing the inequality in a neighborhood of zero. Because z, h, w and ψ ε are continuous, we get, taking h 0 = 0, |h(t)| = |w(0)| As |h 0 | is now positive, we conclude that |h| ≤ ψ ε on a neighborhood of zero by the Cauchy-Schwarz inequality.
The following result is a direct consequence of Theorems C.1 and C.3.
Corollary C.4. Let h 0 ∈ C and z, w : R + → C be continuous functions such that (z) ≤ λ for some λ ∈ R. We define h : R + → C as the unique continuous solution of h = h 0 + K * (zh + w).
Then, for any t ∈ [0, T ], where E λ is the canonical resolvent of K with parameter λ.