Hybrid scheme for Brownian semistationary processes

We introduce a simulation scheme for Brownian semistationary processes, which is based on discretizing the stochastic integral representation of the process in the time domain. We assume that the kernel function of the process is regularly varying at zero. The novel feature of the scheme is to approximate the kernel function by a power function near zero and by a step function elsewhere. The resulting approximation of the process is a combination of Wiener integrals of the power function and a Riemann sum, which is why we call this method a hybrid scheme. Our main theoretical result describes the asymptotics of the mean square error of the hybrid scheme and we observe that the scheme leads to a substantial improvement of accuracy compared to the ordinary forward Riemann-sum scheme, while having the same computational complexity. We exemplify the use of the hybrid scheme by two numerical experiments, where we examine the finite-sample properties of an estimator of the roughness parameter of a Brownian semistationary process and study Monte Carlo option pricing in the rough Bergomi model of Bayer et al. [Quant. Finance 16(6), 887-904, 2016], respectively.


Introduction
We study simulation methods for Brownian semistationary (BSS) processes, first introduced by Barndorff-Nielsen and Schmiegel [8,9], which form a flexible class of stochastic processes that are able to capture some common features of empirical time series, such as stochastic volatility (intermittency), roughness, stationarity and strong dependence. By now these processes have been applied in various contexts, most notably in the study of turbulence in physics [7,16] and in finance as models of energy prices [4,11]. A BSS process X is defined via the integral representation g(t − s)σ(s)dW (s), (1.1) where W is a two-sided Brownian motion providing the fundamental noise innovations, the amplitude of which is modulated by a stochastic volatility (intermittency) process σ that may depend on W . This driving noise is then convolved with a deterministic kernel function g that specifies the dependence structure of X. The process X can also be viewed as a moving average of volatilitymodulated Brownian noise and setting σ(s) = 1, we see that stationary Brownian moving averages are nested in this class of processes.
In the applications mentioned above, the case where X is not a semimartingale is particularly relevant. This situation arises when the kernel function g behaves like a power-law near zero; more specifically, when for some α ∈ (− 1 2 , 1 2 ) \ {0}, Here we write "∝" to indicate proportionality in an informal sense, anticipating a rigorous formulation of this relationship given in Section 2.2 using the theory of regular variation [15], which plays a significant role in our subsequent arguments. The case α = − 1 6 in (1.2) is important in statistical modeling of turbulence [16] as it gives rise to processes that are compatible with Kolmogorov's scaling law for ideal turbulence. Moreover, processes of similar type with α ≈ −0.4 have been recently used in the context of option pricing as models of rough volatility [1,10,18,20], see Sections 2.5 and 3.3 below. The case α = 0 would (roughly speaking) lead to a process that is a semimartingale, which is thus excluded.
Under (1.2), the trajectories of X behave locally like the trajectories of a fractional Brownian motion with Hurst index H = α + 1 2 ∈ (0, 1) \ { 1 2 }. While the local behavior and roughness, measured in terms of Hölder regularity, of X are determined by the parameter α, the global behavior of X (e.g., whether the process has long or short memory) depends on the behavior of g(x) as x → ∞, which can be specified independently of α. This should be contrasted with fractional Brownian motion and related self-similar models, which necessarily must conform to a restrictive affine relationship between their Hölder regularity (local behavior and roughness) and Hurst index (global behavior), as elucidated by Gneiting and Schlather [21]. Indeed, in the realm of BSS processes, local and global behavior are conveniently decoupled, which underlines the flexibility of these processes as a modeling framework.
In connection with practical applications, it is important to be able to simulate the process X. If the volatility process σ is deterministic and constant in time, then X will be strictly stationary and Gaussian. This makes X amenable to exact simulation using the Cholesky factorization or circulant embeddings, see, e.g., [2,Chapter XI]. However, it seems difficult, if not impossible, to develop an exact method that is applicable with a stochastic σ, as the process X is then neither Markovian nor Gaussian. Thus, in the general case one must resort to approximative methods. To this end, Benth et al. [13] have recently proposed a Fourier-based method of simulating BSS processes, and more general Lévy semistationary (LSS) processes, which relies on approximating the kernel function g in the frequency domain.
In this paper, we introduce a new discretization scheme for BSS processes based on approximating the kernel function g in the time domain. Our starting point is the Riemann-sum discretization of (1.1). The Riemann-sum scheme builds on an approximation of g using step functions, which has the pitfall of failing to capture appropriately the steepness of g near zero. In particular, this becomes a serious defect under (1.2) when α ∈ (− 1 2 , 0). In our new scheme, we mitigate this problem by approximating g using an appropriate power function near zero and a step function elsewhere. The resulting discretization scheme can be realized as a linear combination of Wiener integrals with respect to the driving Brownian motion W and a Riemann sum, which is why we call it a hybrid scheme. The hybrid scheme is only slightly more demanding to implement than the Riemann-sum scheme and the schemes have the same computational complexity as the number of discretization cells tends to infinity.
Our main theoretical result describes the exact asymptotic behavior of the mean square error (MSE) of the hybrid scheme and, as a special case, that of the Riemann-sum scheme. We observe that switching from the Riemann-sum scheme to the hybrid scheme reduces the asymptotic root mean square error (RMSE) substantially. Using merely the simplest variant the of hybrid scheme, where a power function is used in a single discretization cell, the reduction is at least 50% for all α ∈ (0, 1 2 ) and at least 80% for all α ∈ (− 1 2 , 0). The reduction in RMSE is close to 100% as α approches − 1 2 , which indicates that the hybrid scheme indeed resolves the problem of poor precision that affects the Riemann-sum scheme.
To assess the accuracy of the hybrid scheme in practice, we perform two numerical experiments. Firstly, we examine the finite-sample performance of an estimator of the roughness index α, introduced by Barndorff-Nielsen et al. [6] and Corcuera et al. [16]. This experiment enables us to assess how faithfully the hybrid scheme approximates the fine properties of the BSS process X. Secondly, we study Monte Carlo option pricing in the rough Bergomi stochastic volatility model of Bayer et al. [10]. We use the hybrid scheme to simulate the volatility process in this model and we find that the resulting implied volatility smiles are indistinguishable from those simulated using a method that involves exact simulation of the volatility process. Thus we are able propose a solution to the problem of finding an efficient simulation scheme for the rough Bergomi model, left open in the paper [10].
The rest of this paper is organized as follows. In Section 2 we recall the rigorous definition of a BSS process and introduce our assumptions. We also introduce the hybrid scheme, state our main theoretical result concerning the asymptotics of the mean square error and discuss an extension of the scheme to a class of truncated BSS processes. Section 3 briefly discusses the implementation of the discretization scheme and presents the numerical experiments mentioned above. Finally, Section 4 contains the proofs of the theoretical and technical results given in the paper.

Brownian semistationary process
Let (Ω, F, {F t } t∈R , P) be a filtered probability space, satisfying the usual conditions, supporting a (two-sided) standard Brownian motion W = {W (t)} t∈R . We consider a Brownian semistationary process where σ = {σ(t)} t∈R is an {F t } t∈R -predictable process with locally bounded trajectories, which captures the stochastic volatility (intermittency) of X, and g : (0, ∞) → [0, ∞) is a Borel measurable kernel function.
To ensure that the integral (2.1) is well-defined, we assume that the kernel function g is square integrable, that is, ∞ 0 g(x) 2 dx < ∞. In fact, we will shortly introduce some more specific assumptions on g that will imply its square integrability. Throughout the paper, we also assume that the process σ has finite second moments, E[σ(t) 2 ] < ∞ for all t ∈ R, and that the process is covariance stationary, namely, These assumptions imply that also X is covariance stationary, that is, However, the process X need not be strictly stationary as the dependence between the volatility process σ and the driving Brownian motion W may be time-varying.

Kernel function
As mentioned above, we consider a kernel function that satisfies g(x) ∝ x α for some α ∈ (− 1 2 , 1 2 )\{0} when x > 0 is near zero. To make this idea rigorous and to allow for additional flexibility, we formulate our assumptions on g using the theory of regular variation [15] and, more specifically, slowly varying functions.
To this end, recall that a measurable function L : Moreover, a function f (x) = x β L(x), x ∈ (0, 1], where β ∈ R and L is slowly varying at 0, is said to be regularly varying at 0, with β being the index of regular variation. A key feature of slowly varying functions, which will be very important in the sequel, is that they can be sandwiched between polynomial functions as follows. If δ > 0 and L is slowly varying at 0 and bounded away from 0 and ∞ on any interval (u, 1], u ∈ (0, 1), then there exist constants C δ ≥ C δ > 0 such that (2. 2) The inequalities above are an immediate consequence of the so-called Potter bounds for slowly varying functions, see [15, Theorem 1.5.6(ii)] and (4.1) below. Making δ very small therein, we see that slowly varying functions are asymptotically negligible in comparison with polynomially growing/decaying functions. Thus, by multiplying power functions and slowly varying functions, regular variation provides a flexible framework to construct functions that behave asymptotically like power functions. Our assumptions concerning the kernel function g are as follows: where L g : (0, 1] → [0, ∞) is continuously differentiable, slowly varying at 0 and bounded away from 0. Moreover, there exists a constant C > 0 such that the derivative L g of L g satisfies (A2) The function g is continuously differentiable on (0, ∞), so that its derivative g is ultimately monotonic and satisfies (Here, and in the sequel, we use h(x) < ∞. Additionally, analogous notation is later used for sequences and computational complexity.) In view of the bound (2.2), these assumptions ensure that g is square integrable. It is worth pointing out that (A1) accommodates functions L g with lim x→0 L g (x) = ∞, e.g., L g (x) = 1 − log x.
The assumption (A1) influences the short-term behavior and roughness of the process X. A simple way to assess the roughness of X is to study the behavior of its variogram (also called the second-order structure function in turbulence literature) as h → 0. Note that, by covariance stationarity, Under our assumptions, we have the following characterization of the behavior of V X near zero, which generalizes a result of Barndorff-Nielsen [3, p. 9] and implies that X has a locally Hölder continuous modification. Therein, and in what follows, we write a(x) ∼ b(x), x → y, to indicate that lim x→y a(x) b(x) = 1. The proof of this result is carried out in Section 4.1.
which implies that V X is regularly varying at zero with index 2α + 1.
Motivated by Proposition 2.2, we call α the roughness index of the process X. Ignoring the slowly varying factor L g (h) 2 in (2.2), we see that the variogram V (h) behaves like h 2α+1 for small values of h, which is reminiscent of the scaling property of the increments of a fractional Brownian motion (fBm) with Hurst index H = α + 1 2 . Thus, the process X behaves locally like such an fBm, at least when it comes to second order structure and roughness. (Moreover, the factor 1 2α+1 + ∞ 0 ((y + 1) α − y α ) 2 dy coincides with the normalization coefficient that appears in the Mandelbrot-Van Ness representation [24, Theorem 1.3.1] of an fBm with H = α + 1 2 .) Let us now look at two examples of a kernel function g that satisfies our assumptions.
with parameters α ∈ (− 1 2 , 1 2 ) \ {0} and λ > 0, has been used extensively in the literature on BSS processes. It is particularly important in connection with statistical modeling of turbulence, see Corcuera et al. [16], but it also provides a way to construct generalizations of Ornstein-Uhlenbeck (OU) processes with roughness that differs from the usual semimartingale case α = 0, while mimicking the long-term behavior of an OU process. Moreover, BSS and LSS processes defined using the gamma kernel have interesting probabilistic properties, see [25]. An in-depth study of the gamma kernel can be found in [3]. Setting L g (x) := e −λx , which is slowly varying at 0 since lim x→0 L g (x) = 1, it is evident that (A1) holds. Since g(x) decays exponentially fast to 0 as x → ∞, it is clear that also (A3) holds. To verify (A2), note that g satisfies where lim x→∞ (( α x − λ) 2 − α x 2 ) = λ 2 > 0, so g is ultimately increasing with Thus, Example 2.4 (Power-law kernel). Consider the kernel function with parameters α ∈ (− 1 2 , 1 2 ) \ {0} and β ∈ (−∞, − 1 2 ). The behavior of this kernel function near zero is similar to that of the gamma kernel, but g(x) decays to zero polynomially as x → ∞, so it can be used to model long memory. In fact, it can be shown that if β ∈ (−1, − 1 2 ), then the autocorrelation function of X is not integrable. Clearly, (A1) holds with L g (x) := (1 + x) β−α , which is slowly varying at 0 since lim x→0 L g (x) = 1. Moreover, note that we can write where K g (x) := (1 + x −1 ) β−α satisfies lim x→∞ K g (x) = 1. Thus, also (A3) holds. We can check (A2) by computing , so g is ultimately increasing. Additionally, we note that

Hybrid scheme
Let t ∈ R and consider discretizing X(t) based on its integral representation (2.1) on the grid G n (t) := {t, t − 1 n , t − 2 n , . . .} for n ∈ N. To derive our discretization scheme, let us first note that if the volatility process σ does not vary too much, then it is reasonable to use the approximation that is, we keep σ constant in each discretization cell. (Here, and in the sequel, "≈" stands for an informal approximation used for purely heuristic purposes.) If k is "small", then due to (A1) we may approximate as the slowly varying function L g varies "less" than the power function y → y α near zero, cf. (2.2). If k is "large", or at least k ≥ 2, then choosing b k ∈ [k − 1, k] provides an adequate approximation For completeness, we also allow for κ = 0, in which case we require that b 1 ∈ (0, 1] and interpret the first sum on the right-hand side of (2.6) as zero. To make numerical implementation feasible, we truncate the second sum on the right-hand side of (2.6) so that both sums have N n ≥ κ + 1 terms in total. Thus, we arrive at a discretization scheme for X(t), which we call a hybrid scheme, given by is a sequence of real numbers, evaluation points, that must satisfy b k ∈ [k − 1, k] \ {0} for each k ≥ κ + 1, but otherwise can be chosen freely.
As it stands, the discretization grid G n (t) depends on the time t, which may seem cumbersome with regard to sampling X n (t) simultaneously for different times t. However, note that whenever times t and t are separated by a multiple of 1 n , the corresponding grids G n (t) and G n (t ) will intersect. In fact the hybrid scheme defined by (2.7) and (2.8) can be implemented efficiently, as we shall see in Section 3.1, below. Since the degenerate case κ = 0 with b k = k for all k ≥ 1 corresponds to the usual Riemann-sum discretization scheme of X(t) with (Itō type) forward sums from (2.8). Henceforth, we denote the associated sequence {k} ∞ k=κ+1 by b FWD , where the subscript "FWD" alludes to forward sums. However, including terms involving Wiener integrals of a power function given by (2.7), that is having κ ≥ 1, improves the accuracy of the discretization considerably, as we shall see. Having the leeway to select b k within the interval [k − 1, k] \ {0}, so that the function g(t − ·) is evaluated at a point that does not necessarily belong to G n (t), leads additionally to a moderate improvement.
The trunction in the sum (2.8) entails that the stochastic integral (2.1) defining X is truncated at t − Nn n . In practice, the value of the parameter N n should be large enough to mitigate the effect of truncation. To ensure that the truncation point t − Nn n tends to −∞ as n → ∞ in our asymptotic results, we introduce the following assumption: (A4) For some γ > 0, N n ∼ n γ+1 , n → ∞.

Asymptotic behavior of mean square error
We are now ready to state our main theoretical result, which gives a sharp description of the asymptotic behavior of the mean square error (MSE) of the hybrid scheme as n → ∞. We defer the proof of this result to Section 4.2.
Corollary 2.7 (Covariance structure). Suppose that the assumptions of Theorem 2.5 hold. Then for any s, t ∈ R and ε > 0, Proof. Let s, t ∈ R. Applying the Cauchy-Schwarz inequality, we get In Theorem 2.5, the asymptotics of the MSE (2.11) are determined by the behavior of the kernel function g near zero, as specified in (A1). The condition (2.9) ensures that error from approximating g near zero is asymptotically larger than the error induced by the truncation of the stochastic integral (2.1) at t − Nn n . In fact, different kind of asymptotics of the MSE, where truncation error becomes dominant, could be derived when (2.9) does not hold, under some additional assumptions, but we do not pursue this direction in the present paper.
While the rate of convergence in (2.11) is fully determined by the roughness index α, which may seem discouraging at first, it turns out that the quantity J(α, κ, b), which we shall call the asymptotic MSE, can vary a lot, depending on how we choose κ and b, and can have a substantial impact on the precision of the approximation of X. It is immediate from (2.12) that increasing κ will decrease J(α, κ, b). Moreover, for given α and κ, it is straightforward to choose b so that J(α, κ, b) is minimized, as shown in the following result.
, and consequently the asymptotic MSE induced by the discretization, is minimized by the sequence b * given by To understand how much increasing κ and using the optimal sequence b * from Proposition 2.8 improves the approximation, we study numerically the asymptotic root mean square error (RMSE) J(α, κ, b). In particular, we assess how much the asymptotic RMSE decreases relative to RMSE of the forward Riemann-sum scheme (κ = 0 and b = b FWD ) using the quantity The results are presented in Figure 1. We find that employing the hybrid scheme with κ ≥ 1 leads to a substantial reduction in the asymptotic RMSE relative to the forward Riemann-sum scheme when α ∈ (− 1 2 , 0). Indeed, when κ ≥ 1, the asymptotic RMSE, as a function of α, does not blow up as α → − 1 2 , while with κ = 0 it does. This explains why the reduction in the asymptotic RMSE approaches 100% as as α → − 1 2 . When α ∈ (0, 1 2 ), the improvement achieved using the hybrid scheme is more modest, but still considerable. Figure 1 also highlights the importance of using the optimal sequence b * , instead of b FWD , as evaluation points in the scheme, in particular when α ∈ (0, 1 2 ). Finally, we observe that increasing κ beyond 2 does not appear to lead to a significant further reduction. Indeed, in our numerical experiments, reported in Section 3.2 and 3.3 below, we observe that using κ = 1, 2 already leads to good results. Remark 2.9. It is non-trivial to evaluate the quantity J(α, κ, b) numerically. Computing the integral in (2.12) explicitly, we can approximate J(α, κ, b) by This approximation is adequate when α ∈ (− 1 2 , 0), but its accuracy deteriorates when α → 1 2 . In particular, the singularity of the function α → J(α, κ, b) at 1 2 is difficult to capture using J N (α, κ, b) with numerically feasible values of N . To overcome this numerical problem, we introduce a correction term in the case α ∈ (0, 1 2 ). The correction term can be derived informally as follows. By the mean value theorem, and since b * k ≈ k − 1 2 for large k, we have is the Hurwitz zeta function, which can be evaluated using accurate numerical algorithms.
Remark 2.10. Unlike the Fourier-based method of Benth et al. [13], the hybrid scheme does not require truncating the singularity of the kernel function g when α ∈ (− 1 2 , 0), which is beneficial to maintaining the accuracy of the scheme when α is near − 1 2 . Let us briefly analyze the effect of truncating the singularity of g on the approximation error, cf. [13, pp. 75-76]. Consider, for any ε > 0, the modified BSS process defined using the truncated kernel function Adapting the proof of Theorem 2.5 in a straightforward manner, it is possible to show that, under (A1) and (A3), for any t ∈ R. While the rate of convergence, as ε ↓ 0, of the MSE that arises from replacing g with g ε is analogous to the rate of convergence of the hybrid scheme, it is important to note that the factorJ(α) blows up as α ↓ − 1 2 . In fact,J(α) is equal to the first term in the series that defines J(α, 0, b FWD ) and which indicates that the effect of truncating the singularity, in terms of MSE, is similar to the effect of using the forward Riemann-sum scheme to discretize the process when α is near − 1 2 . In particular, the truncation threshold ε would then have to be very small in order to keep the truncation error in check.

Extension to truncated Brownian semistationary processes
It is useful to extend the hybrid scheme to a class of non-stationary processes that are closely related to BSS processes. This extension is important in connection with an application to the so-called rough Bergomi model, which we discuss in Section 3.3, below. More precisely, we consider processes of the form where the kernel function g, volatility process σ and driving Brownian motion W are as before.
We call Y a truncated Brownian semistationary (T BSS) process, as Y is obtained from the BSS process X by truncating the stochastic integral in (2.1) at 0. Of the preceding assumptions, only (A1) and (A2) are needed to ensure that the stochastic integral in (2.14) exists -in fact, of (A2), only the requirement that g is differentiable on (0, ∞) comes into play. The T BSS process Y does not have covariance stationary increments, so we define its (timedependent) variogram as Extending Proposition 2.2, we can describe the behavior of h → V Y (h, t) near zero as follows. The existence of a locally Hölder continuous modification is then a straightforward consequence. We omit the proof of this result, as it would be straightforward adaptation of the proof of Proposition 2.2. (i) The variogram of Y satisfies for any t ≥ 0, is regularly varying at zero with index 2α + 1.
Note that while the increments of Y are not covariance stationary, the asymptotic behavior of V Y (h, t) is the same as that of V X (h) as h → 0 (cf. Proposition 2.2) for any t > 0. Thus, the increments of Y (apart from increments starting at time 0) are locally like the increments of X.
We define the hybrid scheme to discretize Y (t), for any t ≥ 0, as whereY In effect, we simply drop the summands in (2.7) and (2.8) that correspond to integrals and increments on the negative real line. We make remarks on the implementation of this scheme in Section 3.1, below. The MSE of hybrid scheme for the T BSS process Y has the following asymptotic behavior as n → ∞, which is, in fact, identical to the asymptotic behavior of the MSE of the hybrid scheme for BSS processes. We omit the proof of this result, which would be a simple modification of the proof of Theorem 2.5.
Theorem 2.12 (Asymptotics of mean square error). Suppose that (A1) and (A2) hold, and that for some δ > 0, Then for all t > 0, where J(α, κ, b) is as in Theorem 2.5 Remark 2.13. Under the assumptions of Theorem 2.12, the conclusion of Corollary 2.7 holds mutatis mutandis. In particular, the covariance structure of the discretized T BSS process approaches that of Y when n → ∞.

Practical implementation
Simulating the BSS process X on the equidistant grid {0, 1 n , 2 n , . . . , nT n } for some T > 0 using the hybrid scheme entails generating .

(3.4)
In order to simulate (3.2) and (3.3), it is instrumental to note that the κ + 1-dimensional random vectors W n i := W n i , W n i,1 , . . . , W n i,κ , i = −N n , −N n + 1, . . . , nT − 1, are i.i.d. according to a multivariate Gaussian distribution with mean zero and covariance matrix Σ given by for j = 2, . . . , κ + 1, and for j, k = 2, . . . , κ + 1 such that j < k, where 2 F 1 stands for the Gauss hypergeometric function, see, e.g., [17, p. 56] for the definition. (When k < j, set Σ j,k = Σ k,j .) For the convenience of the reader, we provide a proof of (3.5) in Section 4.3. Thus, {W n i } nT −1 i=−Nn can be generated by taking independent draws from the multivariate Gaussian distribution N κ+1 (0, Σ). If the volatility process σ is independent of W , then {σ n i } nT −1 i=−Nn can be generated separately, possibly using exact methods. (Exact methods are available, e.g., for Gaussian processes, as mentioned in the introduction, and diffusions, see [14].) In the case where σ depends on W , simulating {W n i } nT −1 i=−Nn and {σ n i } nT −1 i=−Nn is less straightforward. That said, if σ is driven by a standard Brownian motion Z, correlated with W , say, one could rely on a factor decomposition i=−Nn thereafter. This approach has, however, the caveat that it induces an additional approximation error, not quantified in Theorem 2.5.
Remark 3.1. In the case of the T BSS process Y , introduced in Section 2.5, the observations Y n ( i n ), i = 0, 1, . . . , nT , given by the hybrid scheme (2.15) can be computed via .
In the hybrid scheme, it typically suffices to take κ to be at most 3. Thus, in (3.4), the first sumX n ( i n ) requires only a negligible computational effort. By contrast, the number of terms in the second sumX n ( i n ) increases as n → ∞. It is then useful to note that  is O(n log n). Figure 2 presents examples of trajectories of the BSS process X using the hybrid scheme with κ = 1, 2 and b = b * . We choose the kernel function g to be the gamma kernel (Example 2.3) with λ = 1. We also discretize X using the Riemann-sum scheme, κ = 0 with b ∈ {b FWD , b * } (that is, the forward Riemann-sum scheme and its counterpart with optimally chosen evaluation points). We can make two observations: Firstly, we see how the roughness parameter α controls the regularity properties of the trajectories of X -as we decrease α, the trajectories of X become increasingly rough. Secondly, and more importantly, we see how the simulated trajectories coming from the Riemann-sum and hybrid schemes can be rather different, even though we use the same innovations for the driving Brownian motion. In fact, the two variants of the hybrid scheme (κ = 1, 2) yield almost identical trajectories, while the Riemann-sum scheme (κ = 0) produces trajectories that are comparatively smoother, this difference becoming more apparent as α approaches − 1 2 . Indeed, in the extreme case with α = −0.499, both variants of the Riemann-sum scheme break down and yield anomalous trajectories with very little variation, while the hybrid scheme continues to produce accurate results. The fact that the hybrid scheme is able to reproduce the fine properties of rough BSS processes, even for values of α very close to − 1 2 , is backed up by a further experiment reported in the following section.

Estimation of the roughness parameter
Suppose that we have observations X( i m ), i = 0, 1, . . . , m, of the BSS process X, given by (2.1), for some m ∈ N. Barndorff-Nielsen et al. [6] and Corcuera et al. [16] discuss how the roughness index α can be estimated consistently as m → ∞. The method is based on the change-of-frequency (COF) statistics which compare the realized quadratic variations of X, using second-order increments, with two different lag lengths. Corcuera et al. [16] have shown that under some assumptions on the process X, which are similar to (A1), (A2) and (A3) albeit slightly more restrictive, it holds that An in-depth study of the finite sample performance of this COF estimator can be found in [12].
To examine how well the hybrid scheme reproduces the fine properties of the BSS process in terms of regularity/roughness, we apply the COF estimator to discretized trajectories of X, where the kernel function g is again the gamma kernel (Example 2.3) with λ = 1, generated using the hybrid scheme with κ = 1, 2, 3 and b = b * . We consider the case where the volatility process satisfies σ(t) = 1, that is, the process X is Gaussian. This allows us to quantify and control for the intrinsic bias and noisiness, measured in terms of standard deviation, of the estimation method itself, by initially applying the estimator to trajectories that have been simulated using an exact method based on the Cholesky factorization. We then study the behavior of the estimator when applied to a discretized trajectory, while decreasing the step size of the discretization scheme. More precisely, we simulateα(X n , m), where m = 500 and X n is the hybrid scheme for X with n = ms and s ∈ {1, 2, 5}. This means that we computeα(X n , m) using m observations obtained by subsampling every s-th observation in the sequence X n ( i n ), i = 0, 1, . . . , n. As a comparison, we repeat these simulations substituting the hybrid scheme with the Riemann-sum scheme, using The results are presented in Figure 3. We observe that the intrinsic bias of the estimator with m = 500 observations is negligible and hence the bias of the estimates computed from discretized trajectories is then attributable to approximation error arising from the respective discretization scheme, where positive (resp. negative) bias indicates that the simulated trajectories are smoother (resp. rougher) than those of the process X. Concentrating first on the baseline case s = 1, we note that the hybrid scheme produces essentially unbiased results when α ∈ (− 1 2 , 0), while there is moderate bias when α ∈ (0, 1 2 ), which disappears when passing from κ = 1 to κ = 3, even for values of α very close to 1 2 . (The largest value of α considered in our simulations is α = 0.49; one would expect the performance to weaken as α approaches 1 2 , cf. Figure 1, but this range of parameter values seems to be of limited practical interest.) The standard deviations exhibit a similar pattern. The corresponding results for the Riemann-sum scheme are clearly inferior, exhibiting significant bias, while using optimal evaluation points (b = b * ) improves the situation slightly. In particular, the bias in the case α ∈ (− 1 2 , 0) is positive, indicating too smooth discretized trajectories, which is connected with the failure of the Riemann-sum scheme with α near − 1 2 , illustrated in Figure 2. With s = 2 and s = 5, the results improve with both schemes. Notably, in the case s = 5, the performance of the hybrid scheme even with κ = 1 is on a par with the exact method. However, the improvements with the Riemann-sum scheme are more meager, as considerable bias persists when α is near − 1 2 .

Option pricing under rough volatility
As another experiment, we study Monte Carlo option pricing in the rough Bergomi (rBergomi) model of Bayer et al. [10]. In the rBergomi model, the logarithmic spot variance of the price of the underlying is modelled by a rough Gaussian process, which is a special case of (2.14). By virtue of the rough volatility process, the model fits well to observed implied volatility smiles [10, pp. 893-895].
More precisely, the price of the underlying in the rBergomi model with time horizon T > 0 is defined, under an equivalent martingale measure identified with P, as using the spot variance process Above, S(0) > 0, η > 0 and α ∈ (− 1 2 , 0) are deterministic parameters, and Z is a standard Brownian motion given by where ρ ∈ (−1, 1) is the correlation parameter and {W ⊥ (t)} t∈[0,T ] is a standard Brownian motion independent of W . The process {ξ 0 (t)} t∈[0,T ] is the so-called forward variance curve [10, p. 891], which we assume here to be flat, ξ 0 (t) = ξ > 0 for all t ∈ [0, T ]. We aim to compute using Monte Carlo simulation the price of a European call option struck at K > 0 with maturity T , which is given by (3.10) The approach suggested by Bayer et al. [10] involves sampling the Gaussian processes Z and Y on a discrete time grid using exact simulation and then approximating S and v using Euler discretization. We modify this approach by using the hybrid scheme to simulate Y , instead of the computationally more costly exact simulation. As the hybrid scheme involves simulating increments of the Brownian motion W driving Y , we can conveniently simulate the increments of Z, needed for the Euler discretization of S, using the representation (3.9).   Table 1.
We map the option price C(S(0), K, T ) to the corresponding Black-Scholes implied volatility IV(S(0), K, T ), see, e.g., [19]. Reparameterizing the implied volatility using the log-strike k := log(K/S(0)) allows us to drop the dependence on the initial price, so we will abuse notation slightly and write IV(k, T ) for the corresponding implied volatility. Figure 4 displays implied volatility smiles obtained from the rBergomi model using the hybrid and Riemann-sum schemes to simulate Y , as discussed above, and compares these to the smiles obtained using an exact simulation of Y via Cholesky factorization. The parameter values are given in Table 1. They have been adopted from Bayer et al. [10], who demonstrate that they result in realistic volatility smiles. We consider two different maturities: "short", T = 0.041, and "long", T = 1.
We observe that the Riemann-sum scheme (κ = 0, b ∈ {b FWD , b * }) is able capture the shape of the implied volatility smile, but not its level. Alas, the method even breaks down with more extreme log-strikes (the prices are so low that the root-finding algorithm used to compute the implied volatility would return zero). In contrast, the hybrid scheme with κ = 1, 2 and b = b * yields implied volatility smiles that are indistinguishable from the benchmark smiles obtained using exact simulation. Further, there is no discernible difference between the smiles obtained using κ = 1 and κ = 2. As in the previous section, we observe that the hybrid scheme is indeed capable of producing very accurate trajectories of T BSS processes, in particular in the case α ∈ (− 1 2 , 0), even when κ = 1.

Proofs
Throughout the proofs below, we rely on two useful inequalities. The first one is the Potter bound for slow variation at 0, which follows immediately from the corresponding result for slow variation at ∞ [15, Theorem 1.5.6]. Namely, if L : (0, 1] → (0, ∞) is slowly varying at 0 and bounded away from 0 and ∞ on any interval (u, 1], u ∈ (0, 1), then for any δ > 0 there exists a constant C δ > 0 such that The second one is the elementary inequality which can be easily shown using the mean value theorem. Additionally, we use the following variant of Karamata's theorem for regular variation at 0. Its proof is similar to the one of the usual Karamata's theorem for regular variation at ∞ [15, Proposition 1.5.10].

Proof of Proposition 2.2
Proof of Proposition 2.2. (i) By the covariance stationarity of the volatility process σ, we may express the variogram V (h) for any h ≥ 0 as Invoking (A1) and Lemma 4.1, we find that We may clearly assume that h < 1, which allows us to work with the decomposition According to (A2), there exists M > 1 such that x → |g (x)| is non-increasing on [M, ∞). Thus, using the mean value theorem, we deduce that which in turn implies that Making a substitution y = x h , we obtain By the definition of slow variation at 0, lim h→0 G h (y) = (y + 1) α − y α 2 , y ∈ (0, ∞).
(ii) To show existence of the modification, we need a localization procedure that involves an ancillary process We check first that F is locally bounded under (A1) and (A2), which is essential for localization. To this end, let T ∈ (0, ∞), and write for any t ∈ [−T, T ], and M > 1 is such that x → |g (x)| is non-increasing on [M, ∞), as in the proof of (i).
Since g is continuous on (0, ∞) and σ locally bounded, we have for any t ∈ [−T, T ], Further, when t ∈ [−T, T ], Thus, for any t ∈ [−T, T ] almost surely, as we have where we change the order of expectation and integration relying on Tonelli's theorem and where the final equality follows from the covariance stationarity of σ. So we can conclude that F is indeed locally bounded. Let now m ∈ N and, for localization, define a sequence of stopping times that satisfies τ m,n ↑ ∞ almost surely as n → ∞ since both F and σ are locally bounded. (We follow the usual convention that inf ∅ = ∞.) Consider now the modified BSS process that coincides with X on the stochastic interval −m, τ m,n . The process X † m,n satisfies the assumptions of [5, Lemma 1], so for any p > 0 there exists a constantĈ p > 0 such that Using the upper bound in (2.2), we can deduce from (i) that for any δ > 0 there are constants C δ > 0 and h δ > 0 such that (4.8) Applying (4.8) to (4.7), we get We may note that p(α + 1 2 − δ 2 − 1 p ) > 0 for small enough δ and large enough p and, in particular, as δ ↓ 0 and p ↑ ∞. Thus it follows from the Kolmogorov-Chentsov theorem [22,Theorem 3.22] that X † m,n has a modification with locally φ-Hölder continuous trajectories for any φ ∈ (0, α + 1 2 ). Moreover, a modification of X on R, having locally φ-Hölder continuous trajectories for any φ ∈ (0, α + 1 2 ), can then by constructed from these modifications of X † m,n , m ∈ N, n ∈ N, by letting first n → ∞ and then m → ∞.

Proof of Theorem 2.5
As a preparation, we shall first establish an auxiliary result that deals with the asymptotic behavior of certain integrals of regularly varying functions.
Proof. We only prove (i) as (ii) can be shown similarly. By the definition of slow variation at 0, the function In view of the dominated convergence theorem, it suffices to find an integrable dominant for the functions f n , n ∈ N. The construction of the dominant is quite similar to the one seen in the proof of Proposition 2.2, but we provide the details for the convenience of the reader. Using the Potter bound (4.1) and the inequality (u + v) 2 ≤ 2u 2 + 2v 2 , we find that for any where we choose δ ∈ (0, α + 1 2 ). When k ≥ 2, we have x ≥ 1 and b ≥ 1, so is a bounded function of x on [k − 1, k]. When k = 1, we have x ≤ 1 and b ≤ 1, implying that where 2(α − δ) > −1 with our choice of δ, so f is an integrable function on (0, 1]. Proof of Theorem 2.5. Let t ∈ R be fixed. It will be convenient to write X n (t) as Moreover, we introduce an ancillary approximation of X(t), namely, By Minkowski's inequality, we have which together, after taking squares, imply that where Using the Itō isometry, and recalling that σ is covariance stationary, we obtain To analyze the behavior of D n , we substitute y = ns and invoke (A1), yielding where, by Lemma 4.2(ii), we have lim n→∞ k k−1 y 2α L g (k/n) L g (1/n) − L g (y/n) L g (1/n) 2 dy = 0 for any k = 1, . . . , κ. Thus, we find that lim n→∞ D n n −(2α+1) L g (1/n) 2 = 0. (4.12) The asymptotic behavior of the term D n is more delicate to analyze. By (A1), and substituting y = ns, we can write Let us study the asymptotic behavior of the sum n k=κ+1 A n,k as n → ∞. By Lemma 4.2, we have for any k ∈ N, To be able to then deduce, using the dominated convergence theorem, that we seek a sequence {A k } ∞ k=κ+1 ⊂ [0, ∞) such that 0 ≤ A n,k ≤ A k , k = κ + 1, . . . , n, n ∈ N.

Proof of Equation (3.5)
In order to prove (3.5), we rely on the following integral identity for the Gauss hypergeometric function 2 F 1 .
Proof of Equation (3.5). Let α ∈ (− 1 2 , 1 2 ) \ {0} and let j, k = 2, . . . , κ + 1 be such that j < k. By the Itō isometry, we have where we have substituted x = ns. The second integral above can now be expressed as where the second equality follows from Lemma 4.3, which is applicable to both integrals on the second line as 0 < j − 1 < k − 1 and 0 ≤ j − 2 < k − 2.