Eduardo Abi Jaber
email: eduardo.abi-jaber@univ-paris1.fr.

The characteristic function of Gaussian stochastic volatility models: an analytic expression

Keywords: Gaussian processes, Volterra processes, non-Markovian Stein-Stein/Schöbel-Zhu models, rough volatility

Stochastic volatility models based on Gaussian processes, like fractional Brownian motion, are able to reproduce important stylized facts of financial markets such as rich autocorrelation structures, persistence and roughness of sample paths. This is made possible by virtue of the flexibility introduced in the choice of the covariance function of the Gaussian process. The price to pay is that, in general, such models are no longer Markovian nor semimartingales, which limits their practical use. We derive, in two different ways, an explicit analytic expression for the joint characteristic function of the log-price and its integrated variance in general Gaussian stochastic volatility models. Such analytic expression can be approximated by closed form matrix expressions. This opens the door to fast approximation of the joint density and pricing of derivatives on both the stock and its realized variance using Fourier inversion techniques. In the context of rough volatility modeling, our results apply to the (rough) fractional Stein-Stein model and provide the first analytic formulae for option pricing known to date, generalizing that of Stein-Stein, Schöbel-Zhu and a special case of Heston.

Introduction

In the realm of risk management in mathematical finance, academics and practitioners have been always striving for explicit solutions to option prices and hedging strategies in their models. Undoubtedly, finding explicit expressions to a theoretical problem can be highly satisfying in itself; it also has many practical advantages such as: reducing computational time (compared to brute force Monte-Carlo simulations for instance); achieving a higher precision for option prices and hedging strategies; providing a better understanding of the 1 role of the parameters of the model and the sensitivities of the prices and strategies with respect to them. As one would expect, explicit expressions usually come at the expense of sacrificing the flexibility and the accuracy of the model. In a nutshell, the aim of the present paper is to show that analytic expressions for option prices can be found in a highly flexible class of non-Markovian stochastic volatility models.

From Black-Scholes to rough volatility

In their seminal paper, [START_REF] Black | The pricing of options and corporate liabilities[END_REF] derived closed form solutions for the prices of European call and put options in the geometric Brownian motion model where the dynamics of the stock price S are given by: dS t = S t σdB t , S 0 > 0,

(1.1) with B a standard Brownian motion and σ the constant instantaneous volatility parameter.

Although revolutionary, the model remains very simple: it drifts away from the reality of financial markets characterized by non-Gaussian returns, fat tails of stock prices and their volatilities, asymmetric option prices (i.e. the implied volatility smile and skew). . . see [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF]. Since then a large and growing literature has been developed to refine the [START_REF] Black | The pricing of options and corporate liabilities[END_REF] model. One notable direction is stochastic volatility modeling where the constant volatility σ in (1.1) is replaced by a Markovian stochastic process (σ t ) t≥0 . In their celebrated paper, [START_REF] Stein | Stock price distributions with stochastic volatility: an analytic approach[END_REF] modeled (σ t ) t≥0 by a mean-reverting Brownian motion of the form

dσ t = κ(θ -σ t )dt + νdW t , (1.2)
where W is a standard Brownian motion independent of B. Remarkably, they obtained closed-form expressions for the characteristic function of the log-price, which allowed them to recover the density as well as option prices by Fourier inversion of the characteristic function. Later on the model has been extended by [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF] to account for the leverage effect, i.e. an arbitrary correlation between W and B. Similar formulas for the characteristic function of the log-price to those of Stein-Stein are derived for the non-zero correlation case.

Prior to the extension by [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF], [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] took a slightly different approach to include the leverage effect by introducing a model deeply rooted in the Stein-Stein model. Heston observed that the instantenous variance process V t = σ 2 t in the Stein-Stein model with θ = 0 follows a CIR process thanks to Itô's formula,1 so that the Stein-Stein model can be recast in the following form

dS t = S t V t dB t , dV t = (ν 2 -2κV t )dt + 2ν V t dW t , (1.3)
where B = ρW + 1 -ρ2 W ⊥ with ρ ∈ [-1, 1] and W ⊥ a Brownian motion independent of W . Such model remains tractable as it was shown earlier in the context of bond pricing with uncertain inflation by Cox et al. (1985, Equations ( 51)-( 52)). 2 Heston (1993) carried on by deriving closed form expressions for the characteristic function of the log-price, which made his model one of the most, if not the most, popular model among practitioners. As one would expect, the expressions of [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] and [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF] share a lot of similarities and they perfectly agree when θ = 0 in (1.2), see Lord and Kahl (2006, equation (44)). Such analytical tractability motivated the development of the theory of finite-dimensional Markovian affine processes, see [START_REF] Duffie | Affine processes and applications in finance[END_REF].

Unfortunately, Markovian stochastic volatility models, such as the Heston and the Stein-Stein models, are not flexible enough: they generate an auto-correlation structure which is too simplistic compared to empirical observations. Indeed, several empirical studies have documented the persistence in the volatility time series, see [START_REF] Andersen | Intraday periodicity and volatility persistence in financial markets[END_REF]; [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF]. More recently, [START_REF] Gatheral | Volatility is rough[END_REF] and [START_REF] Bennedsen | Decoupling the short-and long-term behavior of stochastic volatility[END_REF] show that the sample paths of the realized volatility are rougher than standard Brownian motion at any realistic time scale as illustrated on Figure 1-(a). From a pricing perspective, continuous semi-martingale models driven by a standard Brownian motion fail to reproduce the power-law decay of the at-the-money skew of option prices as shown on Figure 1-(b), see [START_REF] Carr | The finite moment log stable process and option pricing[END_REF]; [START_REF] Fouque | Multiscale stochastic volatility asymptotics[END_REF]; [START_REF] Lee | Implied volatility: Statics, dynamics, and probabilistic interpretation[END_REF]; [START_REF] Alòs | On the short-time behavior of the implied volatility for jump-diffusion models with stochastic volatility[END_REF]; [START_REF] Bayer | Pricing under rough volatility[END_REF]; [START_REF] Fukasawa | Asymptotic analysis for stochastic volatility: martingale expansion[END_REF][START_REF] Fukasawa | Volatility has to be rough[END_REF].

Date

(a) S&P realized volatility q q q q q q q q q q 0.5 for the S&P index on June 20, 2018 (red dots) and a power-law fit t → 0.35 × t -0.41 . Here k := ln(K/S 0 ) stands for the log-moneyness and T for the time to maturity.

These studies have motivated the need to enhance conventional stochastic volatility models with richer auto-correlation structures. This has been initiated in [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF] by replacing the driving Brownian motion of the volatility process by a fractional Brownian motion W H :

W H t = 1 Γ(H+1/2) t 0 (t -s) H-1/2 dW s + 1 Γ(H+1/2) 0 -∞ ((t -s) H-1/2 -(-s) H-1/2 )dW s
where H ∈ (0, 1) is the Hurst exponent: H > 1/2 corresponds to positively correlated returns, H < 1/2 to negatively correlated increments and H = 1/2 reduces to the case of standard Brownian motion. Sample paths of W H are locally Hölder continuous of any order strictly less than H, thereby less regular than standard Brownian motion. Initially [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF] considered the case H > 1/2. However, a smaller Hurst index H ≈ 0.1 allows to match exactly the regularity of the volatility time series and the exponent in the power-law decay of the at-the-money skew measured on the market (Figure 1). Consequently models involving the fractional kernel t → t H-1/2 with H < 1/2 have been dubbed "rough volatility models" by [START_REF] Gatheral | Volatility is rough[END_REF].

The price to pay is that, in general, such models are no longer Markovian nor semimartingales, which limits their practical use and make their mathematical analysis quite challenging. This has initiated a thriving branch of research. 3 The need for fast pricing in such non-Markovian models is therefore, more than ever, crucial. One breakthrough in that direction was achieved by El Euch and Rosenbaum (2019) who came up with a rough version of the [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] model after convolving the dynamics (1.3) with a fractional kernel to get

V t = V 0 + 1 Γ(H + 1/2) t 0 (t -s) H-1/2 (θ -κV s )ds + ν V s dW s , (1.4)
for H ∈ (0, 1/2). Remarkably, they show that an analogous formula for the characteristic function of the log price to that of [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] continue to hold modulo a fractional deterministic Riccati equation. From a theoretical perspective, the rough Heston model falls into the broader class of non-Markovian affine Volterra processes developed in Abi Jaber et al. ( 2019); Abi Jaber (2021), and can be recovered as a projection of infinite dimensional Markovian affine processes as illustrated in Abi Jaber and El Euch (2019a); [START_REF] Cuchiero | Generalized Feller processes and markovian lifts of stochastic Volterra processes: the affine case[END_REF]; [START_REF] Gatheral | Affine forward variance models[END_REF].

Although the rough Heston model can be efficiently implemented (Abi Jaber, 2019b; Abi Jaber and El Euch, 2019b; [START_REF] Callegaro | Rough but not so tough: fast hybrid schemes for fractional Riccati equations[END_REF][START_REF] Gatheral | Rational approximation of the rough Heston solution[END_REF], no closed-form solution for the fractional deterministic Riccati equation and whence for the characteristic function is known to date, which has to be contrasted with the conventional [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] model. One possible explanation could be that, unlike the Markovian case, squares of fractional Brownian motion have different dynamics than (1.4), so that the marginals of the process (1.4) are not chi-square distributed, except for the case H = 1/2.

The main objective of the paper is to rely on squares of general Gaussian processes with arbitrary covariance structures by considering the non-Markovian extension of the [START_REF] Stein | Stock price distributions with stochastic volatility: an analytic approach[END_REF] and the [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF] models. We will show that the underlying Gaussianity makes the problem highly tractable and allows to recover analytic expressions for the joint Fourier-Laplace transform of the log price and the integrated variance in general, which would agree with that of Stein-Stein, Schöbel-Zhu and Heston under the Markovian setting. Such models have been already considered several times in the context of non-Markovian and rough volatility literature [START_REF] Cuchiero | Markovian lifts of positive semidefinite affine Volterra type processes[END_REF][START_REF] Gulisashvili | Extreme-strike asymptotics for general Gaussian stochastic volatility models[END_REF][START_REF] Harms | Affine representations of fractional processes with applications in mathematical finance[END_REF][START_REF] Horvath | Asymptotic behaviour of randomised fractional volatility models[END_REF] but there has been no derivation of the analytic form of the characteristic function. Our methodology takes a step further the recent derivation in Abi Jaber (2019a) for the Laplace transform of the integrated variance and that of Abi Jaber et al. (2021) where the Laplace transform of the forward co-variance curve enters in the context of portfolio optimization.

The Gaussian Stein-Stein model and main results

For T > 0, we will consider the following generalized version of the Stein-Stein model:

dS t = S t X t dB t , S 0 > 0,
(1.5)

X t = g 0 (t) + T 0 K(t, s)κX s ds + T 0 K(t, s)νdW s , (1.6)
with B = ρW + 1 -ρ 2 W ⊥ , ρ ∈ [-1, 1], κ, ν ∈ R, g 0 a suitable deterministic input curve, K : [0, T ] 2 → R a measurable kernel and (W, W ⊥ ) a two-dimensional Brownian motion. Under mild assumptions on its covariance function, every Gaussian process can be written in the form (1.6) with κ = 0, see [START_REF] Sottinen | Stochastic analysis of Gaussian processes via Fredholm representation[END_REF]. Such representation is known as the Fredholm representation. We will be chiefly interested in two classes of kernels K:

• Symmetric kernels, i.e. K(t, s) = K(s, t) for all s, t ≤ T , for which the integration in (1.6) goes up to time T , meaning that X is not necessarily adapted to the filtration generated by W . In this case, the stochastic integral

• 0 XdB cannot be defined in a dynamical way as an Itô integral whenever ρ = 0. We will make sense of (1.5)-(1.6) in a static sense in Section 2.

• Volterra kernels, i.e. K(t, s) = 0 whenever s ≥ t, for which integration in (1.6) goes up to time t, which is more in line with standard stochastic volatility modeling and for which the stochastic integral

• 0 XdB can be defined in the usual Itô sense, see Section 3. For instance, the conventional mean reverting Stein-Stein model (1.2) can be recovered by setting g 0 (t) = X 0 -κθt, κ ≤ 0 and by considering the Volterra kernel K(t, s) = 1 s<t . The fractional Brownian motion with a Hurst index H ∈ (0, 1) can be represented using the Volterra kernel

K(t, s) = 1 s<t (t -s) H-1/2 Γ(H + 1 2 ) 2 F 1 H - 1 2 ; 1 2 -H; H + 1 2 ; 1 - t s ,
where 2 F 1 is the Gauss hypergeometric function; and the Riemman-Liouville fractional Brownian motion corresponds to the case K(t, s) = 1 s<t (t -s) H-1/2 /Γ(H + 1/2).

For suitable u, w ∈ C, we provide the following analytical expression for the conditional joint Fourier-Laplace transform of the log-price and the integrated variance: [START_REF] Fredholm | Sur une classe d'équations fonctionnelles[END_REF] determinant (see Appendix A.1), g t the adjusted conditional mean given by

E exp u log S T S t + w T t X 2 s ds F t = exp ( g t , Ψ t g t L 2 ) det (Φ t ) 1/2 , (1.7) with f, h L 2 = T 0 f (s)h(s)ds, det the
g t (s) = 1 t≤s E X s - T t K(s, r)κX r dr F t , s, t ≤ T ;
(1.8) and Ψ t a linear operator acting on L 2 ([0, T ], R) defined by

Ψ t = (id -bK * ) -1 a id -2a Σt -1 (id -bK) -1 , t ≤ T, (1.9)
where K denotes the integral operator induced by K, K * the adjoint operator,4 id denotes the identity operator, i.e. (idf ) = f for all f ∈ L 2 ([0, T ], C),

a = w + 1 2 (u 2 -u), b = κ + ρνu, (1.10)
and Σt the adjusted covariance integral operator defined by

Σt = (id -bK) -1 Σ t (id -bK * ) -1 , (1.11)
with Σ t defined as the integral operator associated with the covariance kernel

Σ t (s, u) = ν 2 T t K(s, z)K(u, z)dz, t ≤ s, u ≤ T, (1.12)
and finally Φ is defined by

Φ t = (id -bK)(id -2a Σt )(id -bK) if K is a symmetric kernel id -2a Σt if K is a Volterra kernel .
At first glance, the expressions for Φ seem to depend on the class of the kernel, but they actually agree. Indeed, for Volterra kernels, i.e. K(t, s) = 0 for s ≥ t, det(id -bK) = det(id -bK * ) = 1 so that using the relation (Simon, 1977, Theorem 3.8)

det((id + F )(id + G)) = det(id + F ) det(id + G): det((id -bK)(id -2a Σt )(id -bK * )) = det(id -2a Σt ).
As already mentioned, we prove (1.7) for two classes of kernels:

• Symmetric nonnegative kernels: we provide an elementary static derivation of (1.7) for t = 0 and κ = 0, based on the spectral decomposition of K which leads to the decomposition of the characteristic function as an infinite product of independent Wishart distributions. The operator Ψ 0 appears naturally after a rearrangement of the terms. The main result is collected in Theorem 2.3.

• Volterra kernels: under some L 2 -continuity and boundedness condition, we adopt a dynamical approach to derive the conditional characteristic function (1.7) via Itô's formula on the adjusted conditional mean process (g t (s)) t≤s . The main result is stated in Theorem 3.3. This is the class of kernels which is more suited for financial applications.

From the numerical perspective, we will show in Section 4.1 that the expression (1.7) lends itself to approximation by closed form solutions using finite dimensional matrices after a straightforward discretization of the operators in the form

E exp u log S T S 0 + w T 0 X 2 s ds ≈ exp T n g n Ψ n 0 g n det(Φ n 0 ) 1/2
where g n ∈ R n and Φ n 0 , Ψ n 0 ∈ R n×n are entirely determined by (g 0 , K, ν, κ, u, w) and det is the standard determinant of a matrix, we refer to Section 4.1. We illustrate the applicability of these formulas on an option pricing and calibration example by Fourier inversion techniques in a (rough) fractional Stein-Stein model in Section 4.2.

Notations

Fix T > 0. We let K denote R or C. We denote by •, • L 2 the following product

f, g L 2 = T 0 f (s) g(s)ds, f, g ∈ L 2 ([0, T ], K) . We note that •, • L 2 is an inner product on L 2 ([0, T ], R), but not on L 2 ([0, T ], C). We define L 2 [0, T ] 2 , K to be the space of measurable kernels K : [0, T ] 2 → K such that T 0 T 0 |K(t, s)| 2 dtds < ∞. For any K, L ∈ L 2 [0, T ] 2 , K we define the -product by (K L)(s, u) = T 0 K(s, z)L(z, u)dz, (s, u) ∈ [0, T ] 2 , (1.13)
which is well-defined in L 2 [0, T ] 2 , K due to the Cauchy-Schwarz inequality. For any kernel K ∈ L 2 [0, T ] 2 , K , we denote by K the integral operator induced by the kernel K that is

(Kg)(s) = T 0 K(s, u)g(u)du, g ∈ L 2 ([0, T ], K) .
K is a linear bounded operator from L 2 ([0, T ], K) into itself. If K and L are two integral operators induced by the kernels K and L in L 2 [0, T ] 2 , K , then KL is also an integral operator induced by the kernel K L.

When K = R, we denote by K * the adjoint kernel of K for •, • L 2 , that is K * (s, u) = K(u, s), (s, u) ∈ [0, T ] 2 ,
and by K * the corresponding adjoint integral operator.

The square-root of a complex number √ z is defined through its main branch, i.e.

√ z = |z|e iarg(z)/2 with z = |z|e iarg(z) such that arg(z) ∈ (-π, π].

Symmetric kernels: an elementary static approach

We provide an elementary static derivation of the joint Fourier-Laplace transform in the special case of symmetric kernels with κ = 0. We stress that, although the case of symmetric kernels is not of interest for practical applications, it naturally leads through direct computations to the analytic expression (1.7) in terms of the operator Ψ given in (1.9). Later on, in Section 3, such expressions are shown to hold in the more practical case of Volterra kernels using a dynamical approach.

Definition 2.1. A linear operator K from L 2 ([0, T ], R) into itself is symmetric nonnega- tive if K = K * and f, Kf L 2 ≥ 0, for all f ∈ L 2 ([0, T ], R). Whenever K is an integral operator induced by some kernel K ∈ L 2 [0, T ] 2 ,
R , we will say that K is symmetric nonnegative. In this case, it follows that K = K * a.e. and

T 0 T 0 f (s) K(s, u)f (u)duds ≥ 0, ∀f ∈ L 2 ([0, T ], R) . K is said to be symmetric nonpositive, if (-K) is symmetric nonnegative.
Throughout this section, we fix T > 0 and we consider the case of symmetric kernels having the following spectral decomposition

K(t, s) = n≥1 λ n e n (t)e n (s), t, s ≤ T, (2.1)
where (e n ) n≥1 is an orthonormal basis of L 2 ([0, T ], R) for the inner product f, g

L 2 = T 0 f (s)g(s)ds and λ 1 ≥ λ 2 ≥ . . . ≥ 0 with λ n → 0, as n → ∞, such that n≥1 λ n < ∞. (2.2)
Such decomposition is possible whenever the operator K is the (nonnegative symmetric) square-root of a covariance operator C which is generated by a continuous kernel. This is known as Mercer's theorem, see Shorack and Wellner (2009, Theorem 1, p.208) and leads to the so-called Kac-Siegert/Karhunen-Loève representation of the process X, see [START_REF] Kac | On the theory of noise in radio receivers with square law detectors[END_REF]; [START_REF] Karhunen | Zur spektraltheorie stochastischer prozesse[END_REF]; [START_REF] Loeve | Probability theory: foundations, random sequences[END_REF]. In this case, one can show that any squareintegrable Gaussian process X with mean g 0 and covariance C admits the representation (1.6) with κ = 0 on some filtered probability space supporting a Brownian motion W , see Sottinen and Viitasaari (2016, Theorem 12). We start by making precise how one should understand (1.5)-(1.6) in the case of symmetric kernels and κ = 0. We rewrite (1.5) in the equivalent form

log S t = log S 0 - 1 2 t 0 X 2 s ds + ρ t 0 X s dW s + 1 -ρ 2 t 0 X s dW ⊥ s .
(2.3)

We fix T > 0, g 0 ∈ L 2 ([0, T ], R) and a complete probability space (Ω, F, Q) supporting a two dimensional Brownian motion (W, W ⊥ ) and, for each t ≤ T , we set

X t = g 0 (t) + T 0 K(t, s)νdW s .
We note that (2.1)-(2.2) imply that K ∈ L 2 ([0, T ] 2 , R) so that the stochastic integral T 0 K(t, s)νdW s is well-defined as an Itô integral for almost every t ≤ T and X has sample paths in

L 2 ([0, T ], R) almost surely. Setting F t = F X t ∨ F W ⊥ t
where (F Y t ) t≥0 stands for the filtration generated by the process Y , we have that W ⊥ is still a Brownian motion with respect to (F t ) t≥0 and, up to a modification, X is progressively measurable5 with respect to (F X t ) t≥0 (and whence w.r.t. the enlarged filtration (F t ) t≥0 ) so that

• 0 X s dW ⊥ s
is well defined as an Itô integral with respect to (F t ) t≥0 . If ρ = 0, (2.3) is therefore welldefined in the classical way. However, for ρ = 0, since X is not necessarily adapted to the filtration generated by W (and vice versa), W is no longer necessarily a Brownian motion with respect to the extended filtration F X ∨ F W , and one cannot make sense of the stochastic integral

• 0 XdW in the usual dynamical sense. We provide a static interpretation of (2.3) valid only at the terminal time T . To this end, since g 0 ∈ L 2 ([0, T ], R), we can write g 0 = n≥1 g 0 , e n e n . Making use of (2.1), we first observe that, an application of Fubini's theorem (Veraar, 2012, Theorem 2.2), justified by the fact that

T 0 n≥1 E T 0 | λ n e n (t)e n (s)| 2 ds dt = n≥1 λ n T 0 e n (t) 2 dt ≤ n≥1 λ n < ∞, yields that X t = g 0 (t) + T 0 K(t, s)νdW s = n≥1 g 0 , e n + λ n νξ n e n (t), dt ⊗ Q -a.e. (2.4)
where ξ n = T 0 e n (s)dW s , for each n ≥ 1. Since (e n ) n≥1 is an orthonormal family in L 2 , (ξ n ) n≥1 is a sequence of independent standard Gaussian random variables that are F W T measurable. We set

N T = n≥1 g 0 , e n + λ n νξ n ξ n .
(2.5)

Remark 2.2. We note that N T plays the role of T 0 X s dW s , since a formal interchange leads to

N T = n≥1 g 0 , e n + λ n νξ n T 0 e n (s)dW s (=) T 0 n≥1 g 0 , e n + λ n νξ n e n (s)dW s (=) T 0 X s dW s .
Obviously, since ξ n are not adapted the integral • 0 ξ n e n (s)dW s cannot be defined in the non-anticipative sense.

Finally, we take as definition for the log-price at the terminal time T :

log S T = log S 0 - 1 2 T 0 X 2 s ds + ρN T + 1 -ρ 2 T 0 X s dW ⊥ s , S 0 > 0, (2.6)
which is an F W T ∨ F T -measurable random variable. We state our main result of the section on the representation of the characteristic function for symmetric kernels.

Theorem 2.3. Let K be as in (2.1), g 0 ∈ L 2 ([0, T ], R) and set κ = 0. Fix u, w ∈ C such that (u) = 0 and (w) ≤ 0. Then,

E exp u log S T S 0 + w T 0 X 2 s ds = exp ( g 0 , Ψ 0 g 0 L 2 ) det (Φ 0 ) 1/2 , (2.7)
with Ψ 0 and Σ0 respectively given by (1.9) and (1.11), for (a, b) as in (1.10) (with κ = 0), that is

a = w + 1 2 (u 2 -u), b = ρνu,
and

Φ 0 = (id -bK)(id -2a Σ0 )(id -bK).
The rest of the section is dedicated to the proof of Theorem 2.3. The key idea is to rely on the spectral decomposition (2.1) to decompose the characteristic function as an infinite product of independent Wishart distributions. The operators Σ0 and Ψ 0 will then appear naturally after a rearrangement of the terms.

In the sequel, to ease notations, we drop the subscript L 2 in the product •, • L 2 . We will start by computing the joint Fourier-Laplace transform of T 0 X 2 s ds, N T . Furthermore, the representation (2.4) readily leads to

T 0 X 2 s ds = n≥1 g 0 , e n + λ n νξ n 2 .
(2.8)

Lemma 2.4. Let K be as in (2.1), g 0 ∈ L 2 ([0, T ], R), set κ = 0 and fix α, β ∈ C such that (α) ≤ 0, (β) = 0.

(2.9)

Then,

E exp α T 0 X 2 s ds + βN T = exp α + β 2 2 n≥1 g 0 ,en 2 1-2βν √ λn-2αν 2 λn n≥1 1 -2βν √ λ n -2αν 2 λ n . (2.10) Proof. Define U T = α T 0 X 2 s ds + βN T . We first observe that (2.9) yields that |exp (U T )| = exp( (U T )) ≤ 1, so that E [exp (U T )
] is finite. By virtue of the representations (2.5) and (2.8), we have

U T = n≥1 α ξ2 n + β ξn ξ n ,
where ξn = g 0 , e n + ν √ λ n ξ n , for each n ≥ 1. Setting Y n = ( ξn , ξ n ) , it follows that (Y n ) n≥1 are independent such that each Y n is a two dimensional Gaussian vector with mean µ n and covariance matrix Σ n given by

µ n = g 0 , e n 0 and Σ n = ν 2 λ n ν √ λ n ν √ λ n 1 .
Furthermore, we have

U T = n≥1 Y n w n Y n , with w n = α β 2 β 2 0 .
By successively using the independence of Y n and the well-known expression for the characteristic function of the Wishart distribution, see for instance Abi Jaber (2019a, Proposition A.1), we get

E [exp(U T )] = E   exp   n≥1 Y n w n Y n     = n≥1 E exp Y n w n Y n = n≥1 exp tr w n (I 2 -2Σ n w n ) -1 µ n µ n det (I 2 -2Σ n w n ) 1/2
.

We now compute the right hand side. We have

(I 2 -2Σ n w n ) = 1 -2αν 2 λ n -βν √ λ n -βν 2 λ n -2αν √ λ n -β 1 -βν √ λ n so that det(I 2 -2Σ n w n ) = 1 -2βν λ n -2αν 2 λ n and (I 2 -2Σ n w n ) -1 = 1 1 -2βν √ λ n -2αν 2 λ n 1 -βν √ λ n βν 2 λ n 2αν √ λ n + β 1 -2αν 2 λ n -βν √ λ n .
Straightforward computations lead to the claimed expression (2.10).

Relying on the spectral decomposition (2.1), we re-express the quantities entering in (2.10) in terms of suitable operators.

Lemma 2.5. Let K be as in (2.1), set κ = 0 and fix α, β ∈ C as in (2.9). Then, the following operator defined by (1.9) with a = α + β 2 2 and b = νβ:

Ψ α,β 0 = (id -bK * ) -1 a id -2 Σ0 a -1 (id -bK) -1 , t ≤ T,
admits the following decomposition

Ψ α,β 0 = n≥1 α + β 2 2 1 -2βν √ λ n -2αν 2 λ n e n , • e n (2.11)
and

det 1 α + β 2 2 Ψ α,β 0 = n≥1 1 1 -2βν √ λ n -2αν 2 λ n , (2.12)
with the convention that 0/0 = 1. In particular,

E exp α T 0 X 2 s ds + βN T = det 1 α + β 2 2 Ψ α,β 0 1/2 exp g 0 , Ψ α,β 0 g 0 . (2.13)
Proof. Throughout the proof, we will make use of the following rule for computing the decomposition of a product of operators in terms of the orthonormal basis (e n ) n≥1 : for K and L in the form K = It follows from (2.1) that

(id -bK) = n≥1 1 -b λ n e n , • e n . Since (β) = 0, (1 -b √ λ n ) = 1 = 0 for each n ≥ 1, so that (id -bK) is invertible with an inverse given by (id -bK) -1 = n≥1 1 1 -b √ λ n e n , • e n .
(2.14)

Similarly, recalling (1.12), (2.1) leads to the representation of Σ 0 = ν 2 KK * :

Σ 0 = n≥1 ν 2 λ n e n , • e n ,
so that Σ0 given by (1.11) reads

Σ0 = n≥1 ν 2 λ n 1 -b √ λ n 2 e n , • e n .
Whence,

id -2a Σ0 = n≥1 1 -b √ λ n 2 -2aν 2 λ n 1 -b √ λ n 2 e n , • e n .
Recalling that a = α

+ β 2 2 and b = νβ, 1 -b √ λ n 2 -2aν 2 λ n = 1 -2νβ √ λ n -2αν 2 λ n .
Since (α) ≤ 0 and (β) = 0, we have that (1

-2νβ √ λ n -2αν 2 λ n ) > 0 so that id -2a Σ0 is invertible with an inverse given by id -2a Σ0 -1 = n≥1 1 -νβ √ λ n 2 1 -2νβ √ λ n -2αν 2 λ n e n , • e n .
The representations (2.11)-(2.12) readily follows after composing by (id -bK * ) -1 a from the left, by (id -bK) -1 from the right and recalling (2.14). Finally, combining these expressions with (2.10), we obtain (2.13). This ends the proof.

We can now complete the proof of Theorem 2.3.

Proof of Theorem 2.3. It suffices to prove that

E exp u log S T S 0 + w T 0 X 2 s ds = E exp α T 0 X 2 s ds + βN T , (2.15)
where

α = w + 1 2 (u 2 -u) - ρ 2 u 2 2 and β = ρu.
Indeed, if this the case, then

(α) = (w) + 1 2 (ρ 2 -1) (u) 2 ≤ 0,
so that an application of Lemma 2.5 yields the expression (2.7). It remains to prove (2.15) by means of a projection argument. Conditional on (F X t ∨ F W t ) t≤T , by independence of X and W ⊥ , the random variable T 0 X s dW ⊥ s is centered gaussian with variance T 0 X 2 s ds so that

M T : = E exp u 1 -ρ 2 T 0 X s dW ⊥ s (F X t ∨ F W t ) t≤T = exp u 2 (1 -ρ 2 ) 2 T 0 X 2 s ds .
(2.16)

A successive application of the tower property of the conditional expectation on the expression (2.6) yields that

E exp u log S T S 0 + w T 0 X 2 s ds = E E exp u log S T S 0 + w T 0 X 2 s ds (F X t ∨ F W t ) t≤T = E exp w - u 2 T 0 X 2 s ds + ρu T 0 X s dW s M T
leading to (2.15) due to (2.16). This ends the proof.

3 Volterra kernels: a dynamical approach

In this section, we treat the class of Volterra kernels which are practically relevant in mathematical finance. We will consider the class of Volterra kernels of continuous and bounded type in L 2 in the terminology of Gripenberg et al. (1990, Definitions 9.2.1, 9.5.1 and 9.5.2).

Definition 3.1. A kernel K : [0, T ] 2 → R is a Volterra kernel of continuous and bounded type in L 2 if K(t, s) = 0 whenever s ≥ t and

sup t∈[0,T ] T 0 |K(t, s)| 2 ds < ∞, lim h→0 T 0 |K(u + h, s) -K(u, s)| 2 ds = 0, u ≤ T. (3.1)
The following kernels are of continuous and bounded type in L 2 .

Example 3.2. (i) Any convolution kernel of the form

K(t, s) = k(t -s)1 s<t with k ∈ L 2 ([0, T ], R). Indeed, sup t≤T T 0 |K(t, s)| 2 ds = sup t≤T t 0 |k(t -s)| 2 ds ≤ T 0 |k(s)| 2 ds < ∞
yielding the first part of (3.1). The second part follows from the L 2 -continuity of k, see (Brezis, 2010, Lemma 4.3).

(ii) For H ∈ (0, 1),

K(t, s) = 1 s<t (t -s) H-1/2 Γ(H + 1 2 ) 2 F 1 H - 1 2 ; 1 2 -H; H + 1 2 ; 1 - t s ,
where 2 F 1 is the Gauss hypergeometric function. Such kernel enters in the Volterra representation (1.6) of the fractional Brownian motion whose covariance function is [START_REF] Decreusefond | Stochastic analysis of the fractional Brownian motion[END_REF]. In this case,

Σ 0 (s, u) = (K K * )(s, u) = 1 2 (s 2H + u 2H -|s -u| 2H ), see
sup t≤T T 0 |K(t, s)| 2 ds = sup t≤T Σ 0 (t, t) ≤ T 2H
and by developing the square

T 0 |K(u + h, s) -K(u, s)| 2 ds = Σ 0 (u + h, u + h) -2Σ 0 (u + h, u) + Σ 0 (u, u)
which goes to 0 as h → 0.

(iii) Continuous kernels K on [0, T ] 2 . This is the case for instance for the Brownian Bridge W T 1 conditioned to be equal to W T 1 0 at a time T 1 : for all T < T 1 , W T 1 admits the Volterra representation (1.6) on [0, T ] with the continuous kernel K(t, s) = 1 s<t (T 1 -t)/(T 1 -s), for all s, t ≤ T .

(iv) If K 1 an K 2 satisfy (3.1) then so does K 1 K 2 by an application of Cauchy-Schwarz inequality.

Throughout this section, we fix a probability space (Ω, F, (F t ) t≤T , Q) supporting a two dimensional Brownian motion (W, W ⊥ ) and we set B = ρW + 1 -ρ 2 W ⊥ . For any Volterra kernel K of continuous and bounded type in L 2 , and any g 0 ∈ L 2 ([0, T ], R), there exists a progressively measurable R × R + -valued strong solution (X, S) to (1.5)-(1.6) such that

sup t≤T E [|X t | p ] < ∞, p ≥ 1, (3.2)
we refer to Theorem A.3 below for the proof. It follows in particular from (3.2) that T 0 X 2 s ds < ∞ almost surely, so that X has sample paths in L 2 ([0, T ], R). We now state our main result on the representation of the Fourier-Laplace transform for Volterra kernels under the following additional assumption on the kernel:

sup t≤T T 0 |K(s, t)| 2 ds < ∞. (3.3) Theorem 3.3. Let g 0 ∈ L 2 ([0, T ], R)
and K be a Volterra kernel as in Definition 3.1 satisfying (3.3). Fix u, w ∈ C, such that 0 ≤ (u) ≤ 1 and (w) ≤ 0. Then,

E exp u log S T S t + w T t X 2 s ds F t = exp (φ t + g t , Ψ t g t L 2 ), (3.4)
for all t ≤ T , with Ψ t given by (1.9) for (a, b) as in (1.10) and

φ t = - T t Tr(Ψ t Σt )dt, (3.5)
where Σt is the strong derivative6 of t → Σ t induced by the kernel

Σt (s, u) = -ν 2 K(s, t)K(u, t), a.e.
and Tr is the trace operator, see Appendix A.1.

Proof. We refer to Appendix B.

The following remark establishes the link between φ and the Fredholm determinant.

Remark 3.4. Assume u, w are real. We recall the definition

Φ t = id -2 Σt a, t ≤ T,
and that Σt is an integral operator of trace class with continuous kernel by virtue of Lemma A.5 below so that the determinant det(Φ t ) is well defined and non-zero by the invertibility of (id -2 Σt a), see Lemma A.6 and Simon (1977, Theorem 3.9). We set

φ t = log(det(Φ t ) -1/2 ) = - 1 2 log(det(Φ t )). (3.6)
Differentiation using the logarithmic derivative of the Fredholm's determinant (see (Gohberg and Krein, 1978, Chap IV, p.158 (1.3))) and (1.11) yields

φt = Tr a id -2 Σt a -1 Σt = Tr a id -2 Σt a -1 (id -bK) -1 Σt (id -bK * ) -1 .
Finally, using (1.9) and the identity Tr(F G) = Tr(GF ), we obtain φt = Tr(Ψ t Σt ).

(3.7)

When u, w are complex numbers, the definition of (3.6) requires the use of several branches of the complex logarithm. For numerical implementation, to prevent complex discontinuities, one should either use (3.6) with multiple branches or stick with the discretization of expression (3.5). We refer to section 4.1 for the numerical implementation.

Finally, for K(t, s) = 1 s<t and an input curve of the form

g 0 (t) = X 0 + θt, t ≥ 0, (3.8)
for some X 0 , θ ∈ R, one recovers from Theorem 3.3 the well-known closed form expressions of [START_REF] Stein | Stock price distributions with stochastic volatility: an analytic approach[END_REF] and [START_REF] Schöbel | Stochastic volatility with an Ornstein-Uhlenbeck process: an extension[END_REF], and that of [START_REF] Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF] when θ = 0.

Corollary 3.5. Assume that K(t, s) = 1 s<t and that g 0 is of the form (3.8), then, the expression (3.4) reduces to

E exp u log S T S t + w T t X 2 s ds F t = exp A(t) + B(t)X t + C(t)X 2 t (3.9)
where A, B, C solve the following system of (Backward) Riccati equations

Ȧ = -θB - 1 2 ν 2 B 2 -ν 2 C, A(T ) = 0, Ḃ = -2θC -(κ + ρνu + 2ν 2 C)B, B(T ) = 0, Ċ = -2ν 2 C 2 -2(κ + ρνu)C -w - 1 2 (u 2 -u), C(T ) = 0.
In particular, (A, B, C) can be computed in closed form as in Lord and Kahl (2006, Equations (43)-( 44)-( 45)).

Sketch of proof. The characteristic function is given by (3.4). Assume that K(t, s) = 1 s<t and g 0 is as in (3.8). Then,

X s = X t + (s -t)θ + s t κX u du + s t νdW u , s ≥ t,
so that taking conditional expectation yields

g t (s) = 1 t≤s (X t + (s -t)θ) .
It follows that

g t , Ψ t g t L 2 = Ã(t) + B(t)X t + C(t)X 2 t with Ã(t) = θ 2 1 t≤• (•-t), Ψ t 1 t≤• (•-t) L 2 , B(t) = 2θ 1 t≤• (•-t), Ψ t 1 t≤• L 2 , C(t) = 1 t≤• , Ψ t 1 t≤• .
Combined with (3.6) and (B.7) below, we obtain (3.9) with A such that A T = 0 and

Ȧ(t) = Ȧ(t) + Tr(Ψ t Σt )
with Tr the trace of an operator (see Appendix A.1 below) and Σt (s, u) = -ν 2 1 t≤s∧u .

Using the operator

Numerical illustration

In this section, we make use of the analytic expression for the characteristic function in (1.7) to price options. We first present an approximation of the formula (1.7) using closed form expressions obtained from a natural discretization of the operators. Throughout this section, we consider the case t = 0 and we fix a Volterra kernel K, i.e. K(t, s) = 0 if s ≥ t, as in Section 3.

A straightforward approximation by closed form expressions

The expression (1.7) lends itself to approximation by closed form solutions by a simple discretization of the operator Ψ 0 given by (1.9) à la [START_REF] Fredholm | Sur une classe d'équations fonctionnelles[END_REF]. Fix n ∈ N and let t i = iT /n, i = 0, 1, . . . , n be a partition of [0, T ]. Discretizing the -product given in (1.13) yields the following approximation for Ψ 0 by the n × n matrix:

Ψ n 0 = a I n -b(K n ) -1 I n -2 aT n Σn -1 (I n -bK n ) -1 ,
where I n is the n × n identity matrix, K n is the lower triangular matrix with components

K n ij = 1 j≤i-1 t j t j-1 K(t i-1 , s)ds, 1 ≤ i, j ≤ n, (4.1)
and Σn = (

I n -bK n ) -1 Σ n I n -b(K n ) -1
with Σ n the n × n discretized covariance matrix, recall (1.12), given by

Σ n ij = ν 2 T 0 K(t i-1 , s)K(t j-1 , s)ds, 1 ≤ i, j ≤ n. (4.2)
Defining the n-dimensional vector g n = (g 0 (t 0 ), . . . , g 0 (t n-1 )) , the discretization of the inner product •, • L 2 leads to the approximation

E exp u log S T + w T 0 X 2 s ds ≈ exp u log S 0 + T n g n Ψ n 0 g n det(Φ n 0 ) 1/2 (4.3) with Φ n 0 = I n -2a T n
Σn .

Remark 4.1. Recalling Remark 3.4, one needs to be careful with the numerical implementation of the square root of the determinant that appears in equation (4.3) to avoid complex discontinuities, either by switching the sign of the determinant each time it crosses the axis of negative real numbers or by discretizing (3.5) which would require the computation of Ψ t for several values of t but has the advantage of being analytic on the whole domain. We refer to [START_REF] Mayerhofer | Reforming the wishart characteristic function[END_REF] for more details for finite-dimensional Wishart distributions.

Remark 4.2. Depending on the smoothness of the kernel, other quadrature rules might be more efficient for the choice of the discretization of the operator and the approximation of the Fredholm determinant based on the so-called Nyström method, see for instance [START_REF] Bornemann | On the numerical evaluation of distributions in random matrix theory: a review[END_REF][START_REF] Bornemann | On the numerical evaluation of Fredholm determinants[END_REF]; [START_REF] Corlay | The Nyström method for functional quantization with an application to the fractional Brownian motion[END_REF]; [START_REF] Kang | Nyström-Clenshaw-Curtis quadrature for integral equations with discontinuous kernels[END_REF].

Remark 4.3. For the case u = 0 and κ = 0, the previous approximation formulas agree with the ones derived in (Abi Jaber, 2019a, Section 2.3) where a numerical illustration for the integrated squared fractional Brownian motion is provided.

Option pricing in the fractional Stein-Stein model

In this section, we illustrate the applicability of our results on the following fractional Stein-Stein model based on the Riemann-Liouville fractional Brownian motion with the Volterra convolution kernel K(t, s) = 1 s<t (t -s) H-1/2 /Γ(H + 1/2):

dS t = S t X t dB t , S 0 > 0, X t = g 0 (t) + κ Γ(H + 1/2) t 0 (t -s) H-1/2 X s ds + ν Γ(H + 1/2) t 0 (t -s) H-1/2 dW s , with B = ρW + 1 -ρ 2 W ⊥ , for ρ ∈ [-1, 1], κ, ν ∈ R
and a Hurst index H ∈ (0, 1). For illustration purposes we will consider that the input curve g 0 , which can be used in general to fit at-the-money curves observed in the market, has the following parametric form7

g 0 (t) = X 0 + 1 Γ(H + 1/2) t 0 (t -s) H-1/2 θds = X 0 + θ t H+1/2 Γ(H + 1/2)(H + 1/2)
. (4.4)

Remark 4.4. It would have also been possible to take instead of the fractional Riemman-Liouville Brownian motion the true fractional Brownian motion by considering

X t = g 0 (t) + ν Γ(H + 1/2) t 0 (t -s) H-1/2 2 F 1 H -1/2, 1/2 -H; H + 1/2, 1 - t s dW s ,
where 2 F 1 is the Gaussian hypergeometric function.

Taking H < 1/2 allows one to reproduce the stylized facts observed in the market as in Figure 1. Indeed, the simulated sample paths of the instantaneous variance process X 2 with H = 0.1 in Figure 2 has the same regularity as the realized variance of the S&P in Figure 1-(a). In the case H < 1/2, we refer to the model as the rough Stein-Stein model. We now move to pricing. The expression (1.7) for the joint characteristic function allows one to recover the joint density p T (x, y) of log S T , T 0 X 2 s ds by Fourier inversion:

p T (x, y) = 1 2π R 2 e -i(z 1 x+z 2 y) E exp iz 1 log S T + iz 2 T 0 X 2 s ds dz 1 dz 2 ,
but also to price derivatives on the stock price and the integrated variance by Fourier inversion techniques, see [START_REF] Carr | Option valuation using the fast Fourier transform[END_REF]; [START_REF] Fang | A novel pricing method for european options based on Fourier-cosine series expansions[END_REF]; [START_REF] Lewis | A simple option formula for general jump-diffusion and other exponential lévy processes[END_REF] among many others. In the sequel we will make use of the cosine method of [START_REF] Fang | A novel pricing method for european options based on Fourier-cosine series expansions[END_REF] to price European call options on the stock S combined with our approximation formulae of Sections 4.1. We start by observing that the kernel Σ 0 is given in the following closed form

Σ 0 (s, u) = ν 2 Γ(H + 1/2) 2 s∧u 0 (s -z) H-1/2 (u -z) H-1/2 dz = ν 2 Γ(α)Γ(1 + α) s α u 1-α 2 F 1 1, 1 -α; 1 + α; s u , s ≤ u,
where α = H + 1/2 and 2 F 1 is the Gaussian hypergeometric function, see for instance Malyarenko (2012, page 71).8 Fix n ∈ N and a given partition 0 = t 0 < t 1 < . . . < t n = T .

It follows that the n × n matrices (4.1)-(4.2) can be computed in closed form:

K n ij = 1 j≤i-1 1 Γ(1 + α) [(t i-1 -t j-1 ) α -(t i-1 -t j ) α ] , 1 ≤ i, j ≤ n, Σ n ij = ν 2 Γ(α)Γ(1 + α) t i-1 α t 1-α j-1 2 F 1 1, 1 -α; 1 + α; t i-1 t j-1 , Σ n ji = Σ n ij , 1 ≤ i ≤ j ≤ n,
with the convention that 0/0 = 0. We note that K n is lower triangular with zeros on the diagonal and that the symmetric matrix Σ n has zeros on its first row and first column.

The final ingredient to compute (4.3) is the vector g n whose elements are given by:

g i n = g 0 (t i-1 ) = X 0 + θ t i-1 α Γ(1 + α) , 1 ≤ i ≤ n.
As a sanity check, we visualize on Figure 3 the convergence of the approximation methods on the implied volatility for H = 0.2 and H = 0.5 with the uniform partition t i = iT /n. The benchmark is computed for H = 0.5 via the cosine method with the closed form expressions for the characteristic function of the conventional Stein-Stein model, see [START_REF] Lord | Why the rotation count algorithm works[END_REF]; and for H = 0.2 using Monte Carlo simulation. The smaller the maturity the faster the convergence. Other discretization rules might turn out to be more efficient and would require less points to achieve the same accuracy, which makes the implementation even faster, recall Remark 4.2. The main challenge for applying such methods is the singularity of the kernel at s = t when H < 1/2 and is left for future research. (ii) the implied volatility surface of the S&P accross several maturities for in Figure 5.

Both calibration lead to Ĥ < 0.5 indicating that the rough regime of the fractional Stein-Stein model is coherent with the observations on the market. q q q q q q q q q q q q q q q q q q q q 0.5 A Trace, determinants and resolvents

A.1 Trace and determinants

In this section we recall classical results on operator theory in Hilbert spaces regarding mainly their trace and their determinant. For further details we refer to [START_REF] Gohberg | Introduction to the theory of linear nonselfadjoint operators[END_REF]; [START_REF] Gohberg | Traces and determinants of linear operators[END_REF]; [START_REF] Simon | Notes on infinite determinants of Hilbert space operators[END_REF][START_REF] Simon | Trace ideals and their applications[END_REF]; [START_REF] Smithies | Integral equations[END_REF], and also Bornemann (2009, Section 2 and 3). Let A be a linear compact operator acting on L 2 ([0, T ], C). Then, the operator A has a countable spectrum9 denoted by sp(A) = (λ n (A)) n≤N (A) , where N (A) is either a finite integer or infinity. Whenever A is a linear operator induced by a kernel

A ∈ L 2 ([0, T ] 2 , C), A is a Hilbert-Schmidt operator on L 2 ([0, T ], C
) into itself and is in particular compact.

The trace and the determinant are two important functionals on the space of compact operators. Such quantities are defined for operators of trace class. A compact operator A is said to be of trace class if the quantity

Tr A = n≥1 Av n , v n (A.1)
is finite for a given orthonormal basis (v n ) n≥1 . It can be shown that the quantity on the right hand side of (A.1) is independent of the choice of the orthonormal basis and will be called the trace of the operator A. Furthermore, Lidskii's theorem [START_REF] Simon | Trace ideals and their applications[END_REF] see Brislawn (1988, Proposition 3).

Furthermore, the equivalence The determinant (A.3) is named after [START_REF] Fredholm | Sur une classe d'équations fonctionnelles[END_REF] who defined it for the first time for integral operators with continuous kernels.

n≥1 (1 + |λ n |) < ∞ ⇔ n≥1 |λ n | < ∞,

A.2 Resolvents

For a kernel K ∈ L 2 ([0, T ] 2 , C), we define its resolvent R T ∈ L 2 ([0, T ] 2 , C) by the unique solution to

R T = K + K R T , K R T = R T K. (A.4)
In terms of integral operators, this translates into

R T = K + KR T , KR T = R T K. (A.5)
In particular, if K admits a resolvent, (id -K) is invertible and

(id -K) -1 = id + R T . (A.6)
Lemma A.2. Any K as in Definition 3.1 admits a resolvent kernel R T which is again a Volterra kernel and satisfies (3.1).

Proof. It follows from (3.1) that K is a Volterra kernel of continuous and bounded type in the terminology of Gripenberg et al. (1990, Definitions 9.5.1 and 9.5.2). But since, we are considering kernels on the compact set [0, T ], then every kernel of bounded and continuous type is of bounded and uniformly continuous type, see Gripenberg et al. (1990, p.243, paragraph 1). An application of Gripenberg et al. (1990, Theorem 9.5.5-(ii)), yields that K admits a resolvent kernel R T which is again a Volterra kernel of bounded and continuous type. In particular,

|R T | L 1 := sup t≤T T 0 |R T (t, s)|ds < ∞.
It remains to prove that R T inherits condition (3.1) from K using the resolvent equation (A.4). We first show that

T 0 T 0 |R T (t, s)| 2 dtds < ∞. (A.7)
An application of Jensen's inequality on the normalized measure

|R T (t, z)|dz/ T 0 |R T (t, z )|dz yields T 0 T 0 |(R T K)(t, s)| 2 dtds ≤ |R T | L 1 sup r≤T T 0 |K(r, s)| 2 ds T 0 T 0 |R T (t, z)|dtdz < ∞.
Combined with the resolvent equation (A.4) and the first part of (3.1), we obtain (A.7). Using (A.7) and the Cauchy-Schwarz inequality we now get

T 0 |(K R T )(t, u)| 2 du ≤ sup t ≤T T 0 |K(t , z)| 2 dz T 0 T 0 |R T (u, z)| 2 dudz < ∞, t ≤ T,
which combined with (3.1) and (A.7) gives

sup t≤T T 0 |R T (t, s)| 2 ds < ∞,
which shows that R T satisfies the first condition in (3.1). Finally, another application of the Cauchy-Schwarz inequality, for all t, h ≥ 0, shows that

T 0 |(K R T )(t + h, s) -(K R T )(t, s)| 2 ds ≤ T 0 T 0 |R T (u , s)| 2 du ds T 0 |K(t + h, u) -K(t, u)| 2 du
where the left hand side goes to 0 as h → 0 from the second part of (3.1). Combined with the resolvent equation (A.4), we can deduce that lim h→0 T 0 |R T (t + h, u) -R T (t, u)| 2 du = 0, which yields the second condition in (3.1) for R T .

Using the resolvent we can provide the explicit solution to the system (1.5)-(1.6).

Theorem A.3. Fix T > 0, g 0 ∈ L 2 ([0, T ], R) and a kernel K as in Definition 3.1. Then, there exists a unique progressively measurable strong solution (X, S) to (1.5)-(1.6) on [0, T ] given by

X t = g 0 (t) + t 0 R κ T (t, s)g 0 (s)ds + 1 κ t 0 R κ T (t, s)νdW s , (A.8) S t = S 0 exp - 1 2 t 0 X 2 s ds + t 0 X s dB s , (A.9)
Lemma A.5. Fix b ∈ C and a kernel K as in Definition 3.1. Then, (id-bK) is invertible. Furthermore, for all t ≤ T , Σt given by (1.11) is an integral operator of trace class with continuous kernel and can be re-written in the form

Σt = (id -bK t ) -1 Σ t (id -bK * t ) -1 (A.10)
where K t is the integral operator induced by the kernel K t (s, u) = K(s, u)1 u≥t , for s, u ≤ T .

Proof. Lemma A.2 yields the existence of the resolvent R b T of bK which is again a Volterra kernel of continuous and bounded type. Whence, (A.6) yields that (id -bK) is invertible with an inverse given by (id + R b T ). To prove (A.10), we fix t ≤ T and we observe that since Σ t (s, u) = 0 whenever s ∧ u ≤ t, we have

(R b T Σ t )(s, u) = T t R b T (s, z)Σ t (z, u)dz = (R b t,T Σ t )(s, u),
where we defined the kernel

R b t,T (s, u) = R b T (s, u)1 u≥t . Similarly, Σ t (R b T ) * = Σ t (R b t,T ) * . Using the resolvent equation (A.4) of R b T , it readily follows that R b t,T is the resolvent of bK t so that (id -bK t ) -1 = (id + R b t,T
). Combining all of the above leads to

Σt = (id -bK) -1 Σ t (id -bK * ) -1 = (id + R b T )Σ t (id + R b T ) * = Σ t + R b T Σ t + Σ t (R b T ) * + R b T Σ t (R b T ) * = Σ t + R b t,T Σ t + Σ t (R b t,T ) * + R b t,T Σ t (R b t,T ) * (A.11) = (id + R b t,T )Σ t (id + R b t,T ) * = (id -bK t ) -1 Σ t (id + bK * t ) -1 ,
which proves (A.10). Furthermore, it can be readily deduced from (A.11) that Σt is an integral operator of trace class with continuous kernel: the trace class property follows from the fact that the product of two Hilbert-Schmidt operators is of trace class; the continuity of the kernel follows from the fact that both K and R b T satisfy (3.1), recall Lemma A.2.

Lemma A.6. Fix a, b ∈ C such that (a) ≤ -(b) 2 2ν 2 .
Let t ≤ T and K be a kernel as in Definition 3.1. Then, (id -2 Σt a) is invertible and Ψ t given by (1.9) is well-defined. Furthermore, if (a) = (b) = 0 then, Ψ t is a symmetric nonpositive operator in the sense of Definition 2.1.

Proof. • Using Lemma A.5, we write

(id -2a Σt ) = (id -bK t ) -1 A t (id -bK * t ) -1 with A t = (id -bK t ) (id -bK * t ) -2aΣ t = id -bK t -bK * t + b 2 K t K * t -2aΣ t .
It suffices to prove that A t is invertible, that is 0 / ∈ sp(A t ). Taking real parts and observing that Σ

t = ν 2 K t K * t yields (A t ) = id -(b)K t -(b)K * t + (b) 2 K t K * t -(b) 2 K t K * t -2 (a)Σ t = (id -(b)K t ) (id -(b)K t ) * -2 (a) + (b) 2 ν 2 Σ t = I + II
The operator I is symmetric nonnegative and invertible so that sp(I) ⊂ (0, ∞). Furthermore, since 2 (a) + (b) 2 ν 2 ≤ 0 by assumption and Σ t is symmetric nonnegative we have sp(II) ∈ [0, ∞). It follows that sp( (A t )) ∈ (0, ∞), showing that 0 / ∈ sp(A t ) and that A t is invertible. Whence, (id -2a Σt ) is invertible. Combined with Lemma A.5, we obtain that Ψ t is well-defined.

• Assume that (a) = (b) = 0. Σt defined as in (1.11) is clearly a symmetric nonnegative operator with a continuous kernel on [0, T ] 2 , recall Lemma A.4, an application of Mercer's theorem (Shorack and Wellner, 2009, Theorem 1, p.208) yields the existence of an orthonormal basis (e n ) n≥1 of L 2 ([0, T ], R) and nonnegative eigenvalues (λ n ) n≥1 such that

Σt = n≥1 λ n e n , • L 2 e n . Whence, id -2a Σt = n≥1 (1 -2aλ n ) e n , • L 2 e n .
Since a ≤ 0, (1 -2aλ n ) ≥ 1 > 0, for each n ≥ 1, so that the inverse of (id -2a Σt ) is a symmetric nonnegative operator given by id -2a Σt

-1 = n≥1 1 1 -2aλ n e n , • L 2 e n .
Finally, Ψ t is clearly symmetric and for any

f ∈ L 2 ([0, T ], R) f, Ψ t f L 2 = a f , id -2a Σt -1 f L 2 ≤ 0,
with f = (id -bK) -1 f . This shows that Ψ t is nonpositive.

B Proof of Theorem 3.3

This section is dedicated to the proof of Theorem 3.3. We fix T > 0, a Volterra kernel K as in Definition 3.1 satisfying (3.3) and u, w ∈ C, such that 0 ≤ (u) ≤ 1 and (w) ≤ 0. It follows that a, b defined by (1.10) satisfy

(a) + (b) 2 2ν 2 = (w) + 1 2 ( (u) 2 -(u)) + 1 2 (ρ 2 -1) (u) 2 ≤ 0,
so that an application of Lemma A.6 yields that Ψ t is well-defined. We now collect from Abi Jaber et al. (2021, Lemma 5.8) further properties of t → Ψ t . In particular, its link with an operator Riccati equation. We recall that t → Ψ t is said to be strongly differentiable at time t ≥ 0, if there exists a bounded linear operator Ψt from

L 2 ([0, T ], C) into itself such that lim h→0 1 h Ψ t+h -Ψ t -h Ψt op = 0, where G op = sup f ∈L 2 ([0,T ],C) Gf L 2 f L 2 .
Lemma B.1. Fix a kernel K as in Definition 3.1 satisfying (3.3). Then, for each t ≤ T , Ψ t given by (1.9) is a bounded linear operator from L 2 ([0, T ], R) into itself. Furthermore, (i) Ψt := (-aid+Ψ t ) is an integral operator induced by a symmetric kernel ψt (s, u) such that

sup t≤T [0,T ] 2 | ψt (s, u)| 2 dsdu < ∞. (ii) For any f ∈ L 2 ([0, T ], R), (Ψ t f 1 t )(t) =(aid + bK * Ψ t )(f 1 t )(t),
where 1 t : s → 1 t≤s .

(iii) t → Ψ t is strongly differentiable and satisfies the operator Riccati equation

Ψt = 2Ψ t Σt Ψ t , t ∈ [0, T ] (B.1) Ψ T = a (id -bK * ) -1 (id -bK) -1
where Σt is the strong derivative of t → Σ t induced by the kernel Σt (s, u) = -ν 2 K(s, t)K(u, t), a.e.

Proof. The proof follows from a straighforward adaptation of the proof of Abi Jaber et al.

(2021, Lemma 5.6).

Using the previous lemma and observing that the adjusted conditional mean given in (1.8) has the following dynamics g t (s) = 1 t≤s g 0 (s) + • We now integrate in s to obtain the right hand side in (B.4). We let N = {(t, ω) : ∃s ∈ [0, T ] such that (B.5) does not hold}. Then, N is a null set and we fix (t, ω) ∈ [0, T ]×Ω\N .

In the sequel, all the equalities are written for this particular ω. First, using that Σt (s, s) = -ν 2 K(s, t) 2 and recalling that with N t = S u 0 exp -u 2 2 t 0 X 2 s ds + u t 0 X s dB s which can be shown to be a true martingale by a similar argument to that used in Abi Jaber et al. (2019, Lemma 7.3). Finally, we have showed that the local martingale M is bounded by a martingale, which gives that M is also a true martingale. The proof is complete.

Ψ =

Figure 1 :

 1 Figure 1: (a) Realized volatility of the S&P downloaded from https://realized. oxford-man.ox.ac.uk/ with an estimated Hurst index of Ĥ = 0.11. (b) Term structure of the at-the-money skew of the implied volatility ∂σ implicit (k,T ) ∂k k=0 for the S&P index on June 20, 2018 (red dots) and a power-law fit t → 0.35 × t -0.41 . Here k := ln(K/S 0 )

n≥1a

  n e n , • e n , L = n≥1 b n e n , • e n , the composition is given by KL = n≥1 a n e n , m≥1 b m e m e n , • e n = n≥1 a n b n e n , • e n .

Figure 2 :

 2 Figure 2: One simulated sample path of the stock price S and the instantaneous variance process X 2 in the rough Stein-Stein model with parameters: X 0 = 0, 1, κ = 0, θ = 0.01, ν = 0.02, ρ = -0.7 and H = 0.1.

Figure 3 :

 3 Figure 3: Convergence of the implied volatility slices for short (T = 0.05 year) and long maturities (T = 1 year) of the operator discretization of Section 4.1 towards: (i) the explicit solution of the conventional Stein-Stein model (H = 0.5 upper graphs); (ii) the 95% Monte-Carlo confidence intervals (H = 0.2 lower graphs). The parameters are X 0 = θ = 0.1, κ = 0, ν = 0.25 and ρ = -0.7.

Figure 4 :

 4 Figure 4: Term structure of the at-the-money skew for the S&P index on June 20, 2018 (red dots) and for the rough Stein-Stein model with calibrated parameters (4.5) (blue circles with dashed line).

Figure 5 :

 5 Figure 5: The implied volatility surface of the S&P index (red) and the calibrated fractional

  The product of two Hilbert-Schmidt operators K and L is of trace class. If in addition, both K and L are integral operators on L 2 ([0, T ]) induced by K and L, then

  allows one to define a determinant functional for a trace class operator A bydet(id + zA) = N (A) n=1 (1 + zλ n (A)),for all z ∈ C. If in addition A is an integral operator induced by a continuous kernel A, then one can show that det(id + zA) s i , s j )) 1≤i,j≤n ] ds 1 . . . ds n .(A.3)

  we derive in the next lemma the dynamics of t → g t , Ψ t g t L 2 .Lemma B.2. The dynamics of t → g t , Ψ t g t L 2 are given byd g t , Ψ t g t L 2 = g t , Ψt g t L 2 -aX 2 t -2uρνX t (K * Ψ t ) (g t )(t) -Tr Ψ t Σt dt + 2ν ((K * Ψ t ) g t ) (t)dW t , dt × Q -a.e (B.2)An application of Itô's lemma on the square yieldsdḡ t (s) 2 = ν 2 K(s, t) 2 + 2ḡ t (s)K(s, t)κX t dt + 2ḡ t (s)K(s, t)νdW t , dt × Q -a.e. application of the Leibniz rule combined with the fact that ḡt (t) = X t for almost every (t, ω) and Lemma B.1-(iii) yields that t → ( Ψt g t )(s) is a semimartingale on [0, s) with the following dynamicsd( Ψt g t )(s) = -ψt (s, t)X t + = -X t ψt (s, t) + ( Ψt g t )(s) + X t ( Ψt K(•, t)κ)(s) dt + ( Ψt K(•, t)ν)(s)dW t , dt × Q -a.e.where we used that Ψt = Ψt and that K(u, t) = 0 for all u ≤ t. Moreover, the quadratic covariation between t → ḡt (s) and t → ( Ψt g t )(s) is given byd ḡ• (s), ( Ψ• g • )(s) t = ν 2

	so that an T	0	t ψt (s, u)ḡ t (u)du + K(s, u)κX u du +	0 T	t	K(s, u)νdW u ψt (s, u)K(u, t)κX u du dt
				t			t
	+	T	ψt (s, u)K(u, t)νdudW t
	t
	Proof. We first set
	ḡt (s) = g 0 (s) +	0	t	K(s, u)κX u du + 0 T ψt (s, u)K(u, t)K(s, t)dudt t 0 K(s, u)νdW u ,	(B.3)
	so that using Lemma B.1-(i), we can write = -
	g Furthermore, we write
			( Ψt g t )(s) =	T	ψt (s, u)ḡ t (u)du
							t

t , Ψ t g t L 2 = T t aḡ t (s) 2 ds + ḡt (s)( Ψt g t )(s) ds.

The Leibniz rule yields

d g t , Ψ t g t L 2 = -aḡ t (t) 2 -ḡt (t)( Ψt g t )(t) dt + T t d aḡ t (s) 2 + ḡt (s)( Ψt g t )(s) ds, dt × Q a.e.

(B.4)

• We first compute the dynamics of t → aḡ t (s)

2 ds + ḡt (s)( Ψt g t )(s). We fix s ∈ [0, T ]. It follows from (B.3), that dḡ t (s) = K(s, t)κX t dt + K(s, t)νdW t , dt × Q -a.e. T 0 ψt (s, u) Σt (u, s)dudt = -Ψt Σt (•, s) (s)dt.

Whence, combining the previous three identities, we get the dynamics of U t (s) := aḡ t (s) 2 + ḡt (s)( Ψt g t )(s) :

dU t (s) = adḡ t (s) 2 + dḡ t (s)( Ψt g t )(s) + ḡt (s)d( Ψt g t )(s) + d ḡ• (s), ( Ψ• g • )(s) t = aν 2 K(s, t) 2 dt + 2aḡ t (s)K(s, t)κX t dt + X t κK(s, t)( Ψt g t )(s)dt + ḡt (s)( Ψt g t )(s)dt -ḡt (s) ψt (s, t)X t dt + ḡt (s)X t ( Ψt K(•, t)κ)(s)dt -Ψt Σt (•, s) (s)dt + 2aḡ t (s)K(s, t)ν + νK(s, t)( Ψt g t )(s) + ḡt (s)( Ψt K(•, t)ν)

(

s) dW t = I(s) + II(s) + III(s) + IV(s) + V(s) + VI(s) + VII(s) dt + (VIII(s) + IX(s) + X(s)) dW t , dt × Q -a.e. (B.5)

  aid + Ψ (B.6) (3.4) remains valid on {(u, w) ∈ C 2 : 0 ≤ (u) ≤ 1 and (v) ≤ 0}. Fix u ∈ [0, 1], w ∈ R -. SetU t = u log S t + w ds + φ t + g t , Ψ t g t L 2 , (B.8)and M t = exp(U t ). It suffices to prove that M is a martingale. Indeed, if this is the case, then observing that the terminal value of M isM T = u log S T + wand writing the martingale propertyE[M T |F t ] = M t , for t ≤ T , yields (3.4).Step 1. We prove that M is a local martingale by expliciting its dynamics. We first observe thatdM t =M t dU t +Combined with the dynamics (B.2) and the fact that a = w + 1 2 (u 2 -u), we get thatdU t = g t , Ψt g t L 2 -Ψ t ) (g t )(t)) 2 = -4 g t , Ψ t Σt Ψ t g t L 2 ,we get that the drift part in (B.9) is given byM t g t , Ψt -2Ψ t Σt Ψ t g t L 2 + φt -Tr Ψ t Σt = 0,by virtue of the Riccati equations (B.1) and (B.7). This shows that M is a local martingale.Step 2. It remains to argue that the local martingale M is a true martingale. To this end, we fix t ≤ T . An application of the second part of Lemma A.6 yields that Ψ t is a symmetric nonpositive operator so that, recall (B.7), g t , Ψ t g t L 2 ≤ 0 and φ t = -Σs )ds ≤ 0.Whence, since w ≤ 0 and 0 ≤ u ≤ 1, it follows from (B.8) thatU t ≤ u log S t |M t | = exp(U t ) ≤ exp(u log S t ) ≤ N t

	t t -2uρνX Observing that 0 X 2 s T 0 1 2 d U t . X 2 s ds Using (1.5), we have d log S t = -1 2 X 2 u 2 2 X 2 4ν 2 ((K T t t u 2 0 X 2 s ds + u 0 X s dB s ≤ u log S 0 -u 2 2 t 0 X 2 s ds + u t 0 X s dB s Tr(Ψ s = u log S 0 -Therefore,	(B.9)
	t

t dt + ρX t dW t + 1 -ρ 2 X t dW ⊥ t . t (K * Ψ t ) (g t )(t) + φt -Tr Ψ t Σt dt + (ρuX t + 2ν (K * Ψ t ) (g t )(t)) dW t + u 1 -ρ 2 X t dW ⊥ t , so that d U t = u 2 X 2 t + 4ρuνX t (K * Ψ t ) (g t )(t) + 4ν 2 ((K * Ψ t ) (g t )(t)) 2 dt. *

Squares of Brownian motion constitute the building blocks of squared Bessel processes, seeRevuz and Yor (1999, Chapter XI).

The long-term level of the variance ν 2 in (1.3) can be replaced by a more general coefficient θ ≥ 0.

Refer to https://sites.google.com/site/roughvol/home for references.

cf. below for detailed notations.

Every jointly measurable and adapted process admits a progressively measurable modification, see[START_REF] Ondreját | On existence of progressively measurable modifications[END_REF].

See Lemma B.1 below.

In conventional Markovian stochastic volatility models, the input curve g0 is usually in the parametric form (4.4). However, if one is interested in a practical implementation, then more general forms of g0 (non-parametric) would allow more flexibility (by making θ time dependent for instance). The advantage is that g0 can be estimated from the market to match certain term structures today (e.g. term structure of forward variance, etc. . . ). For illustration purposes here, and since a comparison with the standard Stein-Stein model is given, we restrict to such parametric forms of g0.

Note that in the case of Remark 4.4, the expression for the covariance function simplifies to Σ0(s, u) =ν 2 2 (s 2H + u 2H -|s -u| 2H ).

We recall that the spectrum sp(A) is defined as the set of points λ ∈ C for which there does not exist a bounded inverse operator (λid -A) -1 .

where R κ T is the resolvent kernel of κK with the convention that R κ T /κ = K when κ = 0. In particular, (3.2) holds.

Proof. If κ = 0, the existence is trivial. Fix κ = 0. An application of Lemma A.2 on the kernel κK yields the existence of a resolvent R κ T satisfying (3.1). We define X as in (A.8) and we write it in compact form:

where we used the notation R κ T (νdW

We first observe that X admits a progressively measurable modification. Indeed, the stochastic integral N is adapted as an Itô integral and it is mean-square continuous, i.e. E[|N t -N s | 2 ] → 0 as s → t by virtue of Itô's isometry and the fact that R κ T satisfies (3.1) (see Lemma A.2). Therefore, N admits a progressively measurable modification. We now show that X solves (1.6). Using (A.6), composing both sides by (id + R κ T ) -1 = (id -κK) and invoking stochastic Fubini's theorem yield

where we used the resolvent equation (A.5) for the last equality. This shows that

yielding that X is a strong solution of (1.6). Furthermore, (3.2) follows from the fact that sup s≤T T 0 |R κ T (s, u)| 2 du < ∞ combined with the Burkholder-Davis-Gundy inequality. One can therefore define S as in (A.9) and it is immediate that S solves (1.5) by an application of Itô's formula. The uniqueness statement follows by reiterating the same argument backwards: by showing that any solution X to (1.6) is of the form (A.8) using the resolvent equation.

We now justify in the three following lemmas that the quantities (id -bK) and id -2a Σt appearing in the definition of t → Ψ t in (1.9) are invertible so that Ψ t is well-defined for any kernel K as in Definition 3.1.

Proof. An application of the Cauchy-Schwarz inequality yields the first part. The second part follows along the same lines as in the proof of Abi Jaber (2019a, Lemma 3.2).

we obtain that 10

On the other hand, we have

Therefore, summing the above, plugging in (B.4), using Lemma B.1-(ii) and recalling (B.6) and that b = κ + uρν and ḡt (t) = X t yield

leading to the claimed dynamics (B.2).

We can now complete the proof of Theorem 3.3. We recall that φ given in (3.5) solves φt = Tr(Ψ t Σt ).

(B.7)

Proof of Theorem 3.3. It suffices to prove that (3.4) holds for all 0 ≤ u ≤ 1 and w ≤ 0 to obtain the claimed expression by analytic continuation. Indeed, the left hand side in (3.4) is analytic in (u, w) in an open region ( (u), (w)) ∈ (u -, u + )×(w -, w + ) by general results on the analycity of characteristic functions, see Widder (2015, Theorem II.5a). The right hand side is also analytic in (u, w) since resolvents are analytic: they are given by power series. Therefore, if (3.4) holds for all 0 ≤ u ≤ 1 and w ≤ 0, then by analytic continuation 10 The operator Ψt Σt = a Σt + Ψt Σt is of trace class: (i) Σt is of trace class since it can be written as product of two Hilbert-Schmidt integral operators Σt = Kt Kt * with Kt(s, z) = K(s, t)/ √ T , so that (A.2) yields Tr( Σt) = T 0 Σt(s, s)ds; (ii) Ψt Σt is of trace class as product of two Hilbert-Schmidt integral operators so that (A.2) yields Tr( Ψt Σt) = T 0 T 0 ψt(s, z) Σt(z, s)dzds.