On Fractional Gaussian Random Fields Simulations

To simulate Gaussian ﬁelds poses serious numerical problems: storage and computing time. The midpoint displacement method is often used for simulating the fractional Brownian ﬁelds because it is fast. We propose an eﬀective and fast method, valid not only for fractional Brownian ﬁelds, but for any Gaussian ﬁelds. First, our method is compared with midpoint for fractional Brownian ﬁelds. Second, the performance of our method is illustrated by simulating several Gaussian ﬁelds. The software FieldSim is an R package developed in R and C and that implements the procedures on which this paper focuses.

The simulation of fractional Gaussian processes is not difficult in one dimension (e.g. the surveys of Bardet, Lang, Oppenheim, Philippe, and Taqqu (2003); Coeurjolly (2000)). Let us recall the numerical complexity of some classical methods : the Cholesky method has a complexity of O(N 3 ), where N is the size of the simulated sample path, Levinson's a complexity of O(N 2 log N ) and Wood and Chan's a complexity of O(N log N ). These complexities do not pose problems in dimension one, but become problematic in higher dimensions. These methods are indeed no longer tractable: these algorithms are time and memory expensive. Since its introduction by Lévy (1965) for Brownian motion, the random midpoint displacement method has been intensively used for generating fractional Brownian field (Fournier, Fussel, and Carpenter 1982;Peitgen and Saupe 1988;Voss 1985). This method is a rough approximation and is very fast.
Our approach is based on exact simulation plus a fast step, that is an improvement of the midpoint method. Let us be more precise. The aim is to simulate a Gaussian field over a fine grid. We first simulate the field in an exact way on a rough grid via the Cholesky method. We then propose a refinement of the midpoint method. The field is simulated, using a set of neighbors, and not only the nearest neighbors. This requires a computation of local weighting coefficients. This is a major difference with the midpoint method. For midpoint, these local coefficients are fixed to 1/4. For our method, these local coefficients are exactly determined from the second-order structure. Contrary to the midpoint method, this allows us to simulate fields with arbitrary covariance as well rough as smooth.
One then needs to compare our methods called fieldsim with the midpoint method. Since the real aim is to simulate rough phenomenon, we estimate the Hölder index H of the simulation. The closer the estimation of H is, the better the simulation is. Our conclusion is the following. First, on fractional Brownian fields, fieldsim is slightly better than midpoint. Since the midpoint method is faster, this means that midpoint simulator is nevertheless relevant for fractional Brownian fields. Note that the Hölder index is only one criterion that can be used to assess the accuracy of the simulation. The procedure fieldsim produces fields with the desired Hölder index. But we do not focus on the difference between the simulated field covariance and the specified one. Second, we illustrate the good performances of our method by simulating several Gaussian fields, for which midpoint is no longer applicable.
The paper proceeds as follows. In the second section, we present our procedure fieldsim. Then we compare with the midpoint method. In the third section, we present our method on many examples. We recall in the appendix some classical results on the estimation of the Hölder index and our R package. Tables and figures are postponed at the end of the paper.

Method
After introducing some notations, we present the both steps of the procedure fieldsim: accurate and refined steps. Then we recall the random midpoint displacement method which is comparable with the procedure proposed here. We underline the drawbacks of this approach compared with ours.

Notations
Let d be a positive integer and X(·) = X(M ), M ∈ R d , be a real valued non stationary field with zero mean and second order moments. In this paper we are only concerned with the second order properties of the field X(·). It is convenient to use a geometrical approach by considering the following Hilbert space M, with the inner product U, V = E {U V } = Cov {U, V }. The elements of M are the linear combinations, with real coefficients, of elements of X(M ), M ∈ R d and their limits for mean square convergence. So the covariance function R(·, ·) is defined by: This function is nonnegative definite (n.n.d.), that is for all n ≥ 1, for all real scalars λ 1 , . . . , λ n , and for all M 1 , . . . , M n ∈ R d , n i,j=1 Conversely, for any n.n.d. function R(·, ·), there exists an unique centered Gaussian field of second order structure given by R(·, ·).

Accurate simulation step
Let us recall that the goal of this paper is to give a procedure that yields discretization of sample path of the Gaussian field associated with any n.n.d. function R(·, ·). In the sequel, we denote by X(·) this sample path. Here we present the accurate simulation part of our procedure. Given a (regular) space discretization {M i , i ∈ I} of size n I , the problem consists in giving a sample of a centered Gaussian vector of size n I : (X(M i )) i∈I of covariance matrix R given by R i,j = R(M i , M j ), i, j ∈ I. There exist many approaches to do that. We chose the procedure given by Dégerine and Lambert-Lacroix (2003) (see Theorem 2). This algorithm is based on Cholesky decomposition of the matrix R in the sequential manner.

Refined simulation step
We need to introduce some additional notations. Let X X I (M ), denote the orthogonal projection of X(M ) on the closed linear subspace X I = sp{X(M i ), i ∈ I}, i.e. the linear predictor of X(M ) given X(M i ), i ∈ I. The partial innovation X(M ) − X X I (M ) is denoted by ε X I (M ). Since ε X I (M ) is uncorrelated with any variables of the space X I , we can obtain "accurate" simulation of X(M ) by X X I (M )+ V ar(ε X I (M ))U where U is a centered and reduced Gaussian variable independent of X(M i ), i ∈ I. Notice that the coefficients weighting the variables X(M i ), i ∈ I in X X I (M ) and the variance of the partial innovation may be determined from the second order structure of the sequence X(M i ), i ∈ I, X(M ) (see Dégerine and Lambert-Lacroix (2003) for details). The drawback of this approach is when the simulated sequence size increases, we have to stock more and more quantities (filters of several partial innovation and associated variances) and to do more and more calculus. Even if that can be done in the case d = 1, it becomes numerically unfeasible when d ≥ 2. To overcome this problem, a natural approach consists in replacing in the previous procedure the indexes set I by a set of indexes of neighbors of M . We denote by N M this set. Notice that X X N M (M ) is the best linear combination of variables of X N M approximating X(M ) in the sense that the variance of X(M ) − X X N M (M ) is minimum. If we have to use only some variables of the set X N M in order to obtain simulation of X(M ), the best way is to use X X N M (M ) + V ar(ε X N M (M ))U.
Let us remark that such a simulated process does not admit anymore R(·, ·) as a covariance function, but a covariance function that is a good approximation of R(·, ·).

Comparison with the random midpoint displacement method
The midpoint displacement method was developed for fractional Brownian field (see Fournier et al. (1982); Peitgen and Saupe (1988); Voss (1985)). For example in the scalar case (d = 1), the covariance of a standard fractional Brownian motion is given by with t and s in [0, 1]. The parameter H, which is usually called Hurst parameter, is a real in (0, 1]. For H = 1/2, we obtain the Brownian motion (with R(t, s) = inf(t, s)). In the case d > 1, the absolute values in the definition of the covariance function below, are replaced by the Euclidean norm over R d .
This procedure is based on use of a recursive subdivision approach. For example in the scalar case, the procedure begins by assigning the null value to X(0) and to X(1). That means that one simulates more or less the bridge associated with X. To fairly compare the both approaches, we will generate X(1) as a standard Gaussian variable instead of X(1) = 0. Then the segment [0, 1] is divided in half to obtain one new grid point. The process at this new point is simulated as average of its two neighbors X(0) and X(1) plus a Gaussian variable with zero mean and appropriate variance. Subdivision continues to a desired level of recursion and the procedure is repeated for each new sub-grid. Precisely, for the iteration n = 1, 2, . . ., the process at the new grid points 2 n−j , j = 2k + 1, 0 ≤ k ≤ n−1 2 , is given by

is a fractional Brownian motion of Hurst parameter H.
When d = 2, one can use for example rectangular ground plane or triangular surfaces. For the first case, the procedure begins by assigning the null values to each of the four corners of the rectangular ground plane. Similarly to the scalar case, we generate some standard Gaussian variables for the three corners. Then the original boundaries of the ground plane are divided in half to obtain five new grid points. Simulations at these new grid positions are obtained as field averages at the two nearest neighbors of these points plus random value as for the scalar case. For the center point, one can use the field average at the four corners. The procedure is repeated for each new sub-grid. For the second case, two equilateral triangles are used, set side by side to form a parallelogram. At each level of recursion, the triangles are subdivided into successively smaller triangular surfaces.
The first difference between this approach and ours concerns the initialization since the procedure proposed here begins with accurate simulation step. Furthermore, the coefficients weighting the field at the two (resp. four) M neighbors in X X I (M ) (where I is constituted of the M neighbors) are generally not equal to 1/2 (resp. 1/4). Precisely, in the case d = 1 and H = 1/2, let us consider t, s and u such that 0 < s < u < t ≤ 1. We obtain X X {s,t} (u) = 1/2(X(s) + X(t)) and ε X {s,t} (u) = ε X I (u) for any I that contains s and t and does not contain any points in ]s, t[. So we need only two neighbors in order to simulate the process at u. On the other hand, if s = 0 then X X {s,t} (u) = uX(t). If we take in our procedure only the point 0 and 1 for the accurate step and use after two neighbors in the refined step, we obtain a procedure to simulate the process in a accurate way. That is not the case for the random midpoint displacement method. Indeed each time the zero point is concerned as neighbor of u, the coefficient weighted the process at the neighbor on the left-hand side is equal to 1/2 instead of u. Otherwise when d = 1 or H = 1/2, the coefficients used in the projection depend on the position and are not equal to 1/2 (resp. 1/4) when two (resp. four) neighbors are concerned. So the random midpoint displacement method does not use all the information contained in the neighbors. However, the both approaches to simulate fractional Brownian fields (see Section 3.1) leads to comparable results; but the procedure fieldsim can be applied for any Gaussian random field, provided that its covariance function is known.

Numerical results
In this section, we illustrate the method proposed here through simulations for several classes of fields.

Fractional Brownian fields
The standard fractional Brownian fields are defined through their covariance function (e.g. Samorodnitsky and Taqqu (1994)): where the Hurst parameter H is real in (0, 1]. The case H = 1 is degenerated and will be omitted in the following. For various size of sample (N + 1) 2 (N = 64, 128, 256 and 512) and various values of H (H = 0.1, 0.3, 0.5, 0.7 and 0.9), we generate 100 paths of fractional Brownian fields, discretized uniformly on [0, 1] 2 . We use the both procedures midpoint and fieldsim. For all the simulations generated by fieldsim, we have chosen in the accurate simulation step, a regular space discretization of size 25 and in the refined simulation step a number of neighbors equal to 8. Figure 1 summarizes typical paths for N = 64, H = 0.1, 0.5 and 0.9. The lower H, the more irregular is the field. For each path, the parameter H is estimated via the generalized quadratic variations introduced in Istas and Lang (1997). We denote byĤ N this estimator. Some Boxplots (see Figure 2) illustrate the results. In order to explore the quality of these both simulation methods, we propose to use the same approach as in Coeurjolly (2000) consisting in testing the value of H. Precisely, we have the following asymptotic normality result: where the convergence is in distribution and γ 2 H is some constant. Some details such as the expression ofĤ N and the one of γ 2 H are given in the Appendix. Table 1 gives the constant γ 2 H for different values of H.   The both methods work well since the level 95% is almost always reached. One can notice that fieldsim is slightly better than midpoint. However, considering the computing times (see Table 3), it seems preferable to use midpoint to simulate fractional Brownian fields.

Two parameters fractional Brownian fields
The standard bi-fractional Brownian fields are defined through their covariance function (see Houdré and Villa (2003)): where the Hurst parameter H is real in (0, 1) and K in (0, 1]. For K = 1 the corresponding process is a standard fractional Brownian field of Hurst parameter H, but, when K < 1, increments of the process start to be non stationary but remain locally self-similar of order HK.
Different paths for size sample N = 64 and parameters H = 0.5, K = 1 and H = 0.9, K = 0.55 are plotted on Figure 5.

Fractional Brownian sheets
The standard fractional Brownian sheets are defined through their covariance function (see Kamont (1996)): Fractional Brownian sheets do not have stationary increments but have stationary increments with respect to each variable. Therefore they are anisotropic fields but they are self-similar of index d i=1 H i . Typical paths of corresponding fields of parameters H 1 = 0.5, H 2 = 0.5 and H 1 = 0.5, H 2 = 0.9 are also plotted on Figure 5.

Space-Time deformed fractional Brownian fields
The space-time fractional Brownian fields are defined through their covariance function (see Bégyn (2006)): The metric space (D d , ρ) is a model of hyperbolic space (e.g. Helgason (1962)). When d = 2, (D 2 , ρ) is the Poincaré's disk. The Hyperbolic fractional Brownian field (in short HFBF) is the centered Gaussian field with covariance function where O is the origin of R d . The HFBM exists if and only if 0 < H ≤ 1/2 (Istas (2005)  Let X(·) be a fractional Brownian field with Hurst parameter H. Let us suppose that one has a sample path of this process over the grid (k, l) T /N, k, l = 0, . . . , N. We consider the generalized variations based on the discrete second derivatives, where a −1 = a 1 = 1 and a 0 = −2. The associated estimator of the Hurst parameter is defined by (see Istas and Lang (1997) Without restriction, one can consider only even N, soĤ N is well defined. It can be shown (see Biermé (2005)) that , where the convergence is in distribution. The constant γ 2 H is given by: where C 1 , C 2 and C 3 are positive constants given by Constants C 1 , C 2 and C 3 are computed in Biermé (2005) in a general setting. Let us give them in our case in a more tractable way: where , that leads to the expression (1) for C 1 .
Concerning the relationship for C 2 , let us recall that for any centered Gaussian variables X and Y , we have, That leads to V ar (V N (X)) = 2 because of 1 i=−1 a i = 0. Changing variables k − k into k" and l − l into l" provides, what leads to the result provided that the sequence of terms (u k l ) 2 converges. To prove that, we can express u k l as follows, and using Taylor's expansion of the function f (x) = (1 + x) H up to order 4, we obtain On the other hand, for large enough values of k or l, | i ,j i,j (k, l)| can be bounded by 1 and then there exists some constant C such that (u k l ) 2 ≤ C(k 2 +l 2 ) 2H−8 . Since H ∈ (0, 1], the series of terms (k 2 + l 2 ) 2H−8 (consequently (u k l ) 2 ) converges. The expression given for C 3 is obtained in a similar way.

Implementation of FieldSim
FieldSim is a set of R functions that allows performing simulations of Gaussian Fractional Fields with known covariance function. Three classes of functions are implemented: • Simulation functions midpoint and fieldsim that perform simulations of the path of the process. Procedure fieldsim uses Levinson-Durbin algorithm, thanks to the covariance function and the number of accurate and refined steps of the algorithm, parameters set by the user. C2D function has been implemented for the covariance function of multifractional Brownian Fields.
• Estimation functions: quadvar yields estimation of the Hurst parameter of a fractional Brownian field using the quadratic variations method Istas and Lang (1997). locquadvar performs estimation of the multifractional function at a given point for a multifractional Brownian field using the procedure described in Lacaux (2004). H.test gives results relating to the test presented at Subsection 3.1.
• C subroutine vf performs the tasks that are consuming because of the number of loops. gamma2 and quadvaraux are other internal functions.
The R environment is the only user interface. Function fieldsim calls the C subroutine vf whose result is returned to R.
In order to make it easier for the reader not used to R language, we detail the call to functions and the commands used to produce graphical outputs.

Fractional Brownian Fields
To simulate Fractional Brownian Fields (see Section 3.1), one needs to specify the Hurst parameter 0 < H < 1 and the covariance function R(M, M simulates the corresponding process, with Elevel accurate simulation steps and Rlevel refined simulation steps with nbneighbor neighbors. That corresponds to a sample path of size (2 Elevel+Rlevel + 1) 2 with a grid of size (2 Elevel + 1) 2 for the accurate simulation steps.
The result of the function fieldsim is a R Object of class list. It contains the following elements: • vectors Zrow, Zcol of x and y coordinates and matrix Z of the simulated path of the process; • real time that gives the CPU time.
Thus to produce the graph such as Figure 3, one calls the function persp as

Two Parameters Fractional Brownian Fields
To simulate Two Parameters Fractional Brownian Fields (see Figure 5), one needs to specify the parameters 0 < H < 1, 0 < K ≤ 1 and the covariance function (see Section 3.3). For instance,

Space-Time deformed Fractional Brownian Fields
To simulate Space-time deformed Fractional Brownian Fields (see Figure 5), one has to specify the Hurst parameter 0 < H < 1, functional parameters τ and σ and the covariance function (see Section 3.5). For instance