A RATE-OPTIMAL TRIGONOMETRIC SERIES EXPANSION OF THE FRACTIONAL BROWNIAN MOTION

: Let B ( H ) ( t ) ; t 2 [ ¡ 1 ; 1], be the fractional Brownian motion with Hurst parameter H 2 ¡ 12 ; 1 ¢ . In this paper we present the series representation where a j ; j 2 N [ f 0 g , are constants given explicitly, and » j ; j 2 N [ f 0 g , e » j ; j 2 N , are independent standard Gaussian random variables. We show that the series converges almost surely in C [ ¡ 1 ; 1], and in mean-square (in L 2 (›)), uniformly in t 2 [ ¡ 1 ; 1]. Moreover we prove that the series expansion has an optimal rate of convergence.


Introduction
Let B (H) (t), t ∈ R, denote the H-fractional Brownian motion (H-FBM, FBM). We recall that the H-FBM is the only Gaussian process with H-self-similar and stationary increments, and almost surely (a.s.) continuous sample paths. H-self-similarity means that for all c > 0 i.e. the processes on the two sides of the sign d = have the same finite dimensional distributions. The parameter H, called Hurst parameter, can be in the interval (0, 1). For H ∈ 1 2 , 1 the stationary increment processes of B (H) (t) are long-memory, while for H ∈ 0, 1 2 they are not. Thus the former case is more important, therefore we consider that case, i.e. for the rest of the paper H ∈ 1 2 , 1 is assumed. It is clear that self-similarity makes it possible to derive the FBM on any finite interval [−c, c] from its version on [−1, 1]. By this reason we regard the statements for B (H) (t) defined on [−1, 1]. Since process B (H) (t), t ∈ [−1, 1], has a.s. square integrable sample paths, i.e. B (H) (·) ∈ L 2 ([−1, 1]), it seems to be useful to deal with its Fourier series expansions. For example, one can evidently use the complete orthonormal system and expand B (H) (t) into the form where coefficients ζ j , j ∈ N∪{0}, ζ j , j ∈ N, are obviously random variables with centered Gaussian distribution. However, they are not uncorrelated, meaning that neither are independent. Thus, the obvious expansion (3) is of little use, e.g. when the aim is to simulate FBM. Of course one could use an other system of deterministic functions of some sort, but it is an important requirement that the coefficient random variables be independent. (E.g. the Karhunen-Loève representation would be suitable, if it were known.) Nevertheless, it is far from trivial to find such an expansion for the FBM. In this paper we present one, which is of the form where Γ, B, and γ are the gamma, beta, and complementary (lower) incomplete gamma functions, respectively. Moreover, ξ j , j ∈ N ∪ {0}, ξ j , j ∈ N, are independent standard Gaussian random variables. Note that if one found the form of (4) in which function t is also expanded (with respect to system (2)), i.e. the form of pure sine and cosine expansion (3), one would have no luck since the coefficient random variables would not be independent. It must be stressed that recently Dzhaparidze and van Zanten have given a series expansion for the FBM, which is similar to (4) without the term a 0 tξ 0 , see [1]. It is based on sine and cosine functions also, but the frequencies are the roots of Bessel functions, so it is somewhat more complicated than (4).
The main result of this paper is Theorem 1, in which representation (4) is stated exactly. The proof is based on the fact that the FBM can be approximated by scaled integrated gamma-mixed Ornstein-Uhlenbeck (ΓMOU) processes; see [3] for the ΓMOU process and this approximation. On the other hand the ΓMOU process can be simply expanded into a pure sine and cosine function series of the form where c j , j ∈ N ∪ {0}, are constants, and ξ j , j ∈ N ∪ {0}, ξ j , j ∈ N, are independent Gaussian random variables. Thus, the scaled integrated ΓMOU process can be represented in the form is used frequently hereafter, which means that the quotient of its two sides converges to 1 when N → ∞. We remark that the above mentioned expansion in [1] is rate-optimal in the sense introduced in [4]; the proof can be found in [2]. Rate-optimality means the following: with independent Gaussian random variables ξ j , j ∈ N, and continuous deterministic functions f j (t), j ∈ N. Such a series expansion is called rate-optimal if it holds for its convergence speed that In [4] it is proven that this convergence speed is really optimal, i.e. there can be no series expansion of the form (5) with a convergence speed faster than (6), and convergence speed (6) can be achieved. Now expansion (4) is also rate-optimal. This statement as well as the one about the mean square convergence rate are formulated in Theorem 2. The method for proving the rate-optimality part is taken from [2]. The outline of the rest of the paper is the following. In Section 2 the main theorem, Theorem 1 is stated and proved. Section 3 includes Theorem 2 about the rates of convergences. In Section 4 an application, the simulation of FBM, is presented.

The series expansion
Let us first fix a few notations and conditions. The stochastic processes in this paper are considered on the time interval [−1, 1]. The beta function is denoted by B(p, q), are the incomplete gamma function, and the complementary incomplete gamma function, respectively, with complex argument z. In both cases the principal branches are considered. Do not confuse the former gamma function with the gamma distribution Γ(p, λ) with density function which is also used in the following. The space of continuous functions on [−1, 1] is denoted by C[−1, 1]. Every random element is defined on the same probability space Ω. L 2 (Ω) denotes the space of square integrable random variables. See [3] for the ΓMOU process to be discussed later. The notation ΓMOU H − 1 2 , λ is also used when it is important to indicate the long-memory parameter H − 1 2 ∈ 0, 1 2 2 and scale parameter λ > 0. The ΓMOU H − 1 2 , λ process is denoted by Y λ (t). Not all the constants are denoted individually; in many cases they act simply as 'const.', always understood as positive constants. The main statement of the paper is the following theorem. where Proof. The proof is based on the fact that the H-FBM can be approximated by scaled integrated ΓMOU H − 1 2 , λ processes, see [3]. On the other hand the ΓMOU H − 1 2 , λ process Y λ (t) can be expanded simply into a pure sine and cosine function series of the form where c j (λ), j ∈ N ∪ {0}, are constants, and ξ j , j ∈ N ∪ {0}, ξ j , j ∈ N, are independent standard Gaussian random variables. Hence, the scaled integrated ΓMOU H − 1 2 , λ process can be represented in the form These steps will be detailed in the following.
First the ΓMOU process Y λ (t) is expanded into a sine and cosine function series. For this purpose let X α (t) be a stationary Ornstein-Uhlenbeck (OU) process with parameter α < 0, that is, where B(t) is a Brownian motion. The covariance function of process X α (t), being an even function in L 2 [−1, 1], has the following Fourier cosine series expansion: Now we use the fact that the covariance function of the ΓMOU H − 1 2 , λ process is the mixture of OU covariance functions, where the mixing distribution has the density function see [3], Section 2.1. The mixing parameter is −α. Thus, using (11), for the covariance function of the ΓMOU H − 1 2 , λ process Y λ (t)we have where The exponent 2 used in the previous notations can be clarified by observing the nonnegativity of the right hand sides of (14) and (15). It is also easy to see that in (13) it was legitimate to change the order of the integral and the infinite sum by Fubini's theorem. Due to (13) R Y λ (t − s) can be factorized as which results in the sine and cosine series representation with independent standard Gaussian random variables ξ j , j ∈ N ∪ {0}, ξ j , j ∈ N. Here the function series also converges a.s. in L 2 [−1, 1], on the one hand by virtue of Parseval's relation, on the other hand because owing to Kolmogorov's two series theorem. Now let us integrate (17) to obtain When obtaining (18) from (17), it was legitimate to change the order of the integral and the infinite sum, because-as it was mentioned above-the function series in (17) converges a.s. in L 2 [−1, 1], thus its integral from zero to t, considered as an L 2 [−1, 1] inner product converges a.s.. In the next step the convergence (see [3] Section 3, Property 4) is applied. Together with (18) it yields Now the question is whether the following limits exist, and if they do, whether they have any closed form: lim Using (14) and (12), in the case j = 0 we obtain by Lebesgue's dominated convergence theorem. After a change of variables the inner integral can be put into closed form: Similarly for j ∈ N we have also by Lebesgue's dominated convergence theorem. Using the Fourier transform and a bit of integral calculus one can show that ∞ 0 e −tx 1 The substitutions b = jπ, t = 0, and b = jπ, t = 1 in (25) lead to respectively. Thus, (24) can be rewritten as To sum up, the answers to the question raised in (20) are the following: Returning to (19) we obtain the representation of B (H) (t), t ∈ [−1, 1]: which coincides with (8-9-10). The limit and the infinite sum in (27) could be changed for the following reason. On the one hand for each j ∈ N the convergence in (26) is monotone increasing (see (23)), and on the other hand for all t ∈ [−1, 1] the series in (28) is a.s. absolutely convergent. The latter a.s. absolute convergence follows from the fact that and from Kolmogorov's two series theorem. The series in (28) converges in L 2 (Ω), uniformly in t ∈ [−1, 1]. To prove this, take into account that the terms of the series in (28) are orthogonal random variables, and use (29) to obtain and from Kolmogorov's two series theorem.

Rate of convergence
Let us consider the series expansion of Theorem 1, and introduce a notation for the truncated series: The following theorem is about the rate of mean square (L 2 (Ω)) convergence and C[−1, 1] convergence of B 2) The rate of convergence in C[−1, 1] can be characterized by the asymptotic relation where the inequality follows from (30).
Clearly we have S n (t) = S n,2 n−1 (t), and where n ∈ N is such that 2 n−1 ≤ N + 1 < 2 n . First we prove that The two terms on the right hand side are estimated separately. The first term can be bounded by using the inequality of Lemma 2.2.2 in [5] (with ψ(x) = e x 2 − 1 there) and the equivalence of the norms of Gaussian random variables: Similarly as in (34) we have for the first term of (38). To estimate the second term of (38) we first mention that by (29) Thus we have proved (37).
It seems now that we can shortly complete the proof of (33), since by (36) we have Now applying (37) which was to be proved.
Apart from the constant factors the convergence rates in Theorem 2 are the same as those in Dzhaparidze and van Zanten's expansion, see [1] and [2]. In the former paper the authors mention, and in the latter they prove, that their series expansion is rate-optimal in the sense given in [4]; see Definition 1. By the second statement of Theorem 2, the expansion of this paper is also rate-optimal.
4 Application: simulation of FBM Series representations are especially usable for simulation. The truncated series representation (31) is of the form of the linear combination of easily computable deterministic functions. Only the numerical computation of coefficients a j may involve difficulties because of the complementary incomplete gamma function. We have used the 'gammainc' function of MATLAB to compute the complementary incomplete gamma function. Originally it accepted real arguments only, while in (31) an imaginary argument is needed. However, after rewriting it works well also for the complex input argument z (see (7)). This is because it uses the continued fraction expansion of the complementary incomplete gamma function, and that expansion is convergent for all complex arguments except for negative reals. The simulation algorithm based on the series expansion of this paper is fast and easy to program. We have written the simulation program, and here we demonstrate the fitting of its output to the FBM. The simulated sample path on the time interval [0, 1] based on 1000 equidistant points was converted by the transformation (1) into a sample path on the time interval [0, 1000], and the difference series was taken.
In Figure 1 one of the two graphs is the estimated spectrum (averaged periodograms of 10 independent realizations) of the difference series. The other curve is the theoretical spectrum of the fractional Gaussian noise (FGN, difference series of the FBM). The value of the Hurst parameter H = 0.9, the accuracy ε = 10 −4 (note that the standard deviation of B (H) (1) is ≈1.39), and the value N = 7308 for the truncation limit was computed from the equation ε = the right hand side of (32).