FracSim: An R Package to Simulate Multifractional Lévy Motions

In this article a procedure is proposed to simulate fractional fields, which are non Gaussian counterpart of the fractional Brownian motion. These fields, called real harmonizable (multi)fractional Levy motions, allow fixing the Holder exponent at each point. FracSim is an R package developed in R and C language. Parallel computers have been also used.


Introduction
Irregular phenomena rise up in various fields of scientific research: fluid mechanics, image processing and financial mathematics for example. From a mathematical point of view, one can measure irregularity by the Hölder exponent, which is a real number 0 < H < 1. The lower H, the more irregular the function, whereas H = 1 corresponds heuristically to differentiable functions, which are considered smooth in this field. Of course, there exist many examples of Hölder functions in mathematical literature. But, non specific Hölder functions are often given by the sample path of a stochastic process.
One of the simplest models is the fractional Brownian motion B H , in short FBM, introduced in Kolmogorov (1940) and further developed in Mandelbrot and Van Ness (1968). Then, the multifractional Brownian motion, in short MBM, has been introduced independently in Peltier and Lévy Véhel (1996) and Benassi, Jaffard, and Roux (1997) as a generalization of the FBM. The FBM and the MBM provide very powerful models in applied mathematics. Whereas the pointwise Hölder exponent of the FBM is almost surely equal to a constant, the MBM one is allowed to vary along the trajectory.
On one hand, simulation of fractional Brownian motion is now both theoretically and practically well understood (see Coeurjolly 2000;Bardet, Lang, Oppenheim, Philippe, and Taqqu 2003, for surveys on this problem). On the other hand, methods of generating sample paths of the multifractional Brownian motion are given in Peltier and Lévy Véhel (1996) and in Chan and Wood (1998).
However, it is a well known fact that in some fields of applications data do not fit Gaussian models, see for instance Mallat (1989), Simoncelli (1999), and Vidakovic (1999) for image modelling. More recently other processes that are neither Gaussian nor stable have been proposed to model Internet traffic (Cohen and Taqqu 2004;Kaj and Taqqu 2004).
Real harmonizable fractional Lévy motions, in short RHFLMs, have been defined in Benassi, Cohen, and Istas (2002) in order to obtain non-Gaussian fields that share many properties with the FBM. In particular, their sample paths are locally Hölder and their pointwise Hölder exponent is almost surely equal to a constant. Then in order to allow the pointwise Hölder exponent to vary along the trajectory, the class of real harmonizable multifractional Lévy motions, in short RHMLMs, has been introduced in Lacaux (2002). This class makes up a class of locally asymptotically self-similar fields that includes RHFLMs and the MBM.
The main aim of this paper is to propose an R package, FracSim, composed of R (R Development Core Team 2005) and C (Kernighan and Ritchie 1988) functions to simulate RHFLMs and RHMLMs. Such motions are described from theory and practical aspects in a first part. Then, in a second part, we give some illustrations of simulated processes in one or two dimension using FracSim. The third part is devoted to implementation of parallel computing that are needed for two-dimensional simulations.

(Multi)fractional Lévy motions 2.1. Theory
In this part we recall some theoretical facts explained in Lacaux (2004). A RHMLM X h with multifractional function h is defined as the stochastic integral: with M (dξ) a random Lévy measure without Brownian component. See Benassi et al. (2002) for a precise definition of random Lévy measures.
In Lacaux (2004), RHMLMs with a Brownian part are considered. Since the Brownian part leads to an independent multifractional Brownian motion for which simulations are known we restrict ourselves to the simulation of the non-Brownian part X h . The field X h is an infinitely divisible field and infinitely divisible laws can be represented as a generalized shot noise series. An overview of these representations is given in Rosiński (2001). The control measure, ν(dz) of the random Lévy measure M (dξ) is in general a sigma finite measure on the set of complex numbers. If it has finite mass denoted by ν(C), the generalized shot noise series can be used for simulation purpose. Let us first introduce definitions of random variables used in these series. Let (Z n ) n≥1 , (U n ) n≥1 and (R n ) n≥1 be independent sequences of random variables.
• ν is supposed to be a finite measure, and (Z n ) n≥1 is a sequence of i.i.d. random variables such that L(Z n ) = ν(dz) ν(C) .
• (U n ) n≥1 is a sequence of i.i.d. random variables such that the law of U 1 is the uniform probability on the unit sphere in R d .
• R n is the nth arrival time of a Poisson process with intensity 1.
Note that when d = 1, U n is a symmetric Bernoulli random variable.
Let us now recall the theoretical result that allows us to simulate X h . Next proposition by Lacaux (2004) used in the following is given without proof.
Proposition 2.1. Let K ⊂ R d be a compact set. Then almost surely, the series where denotes the real part of a complex number, converges uniformly on K and where d = denotes equality in distribution.
The RHMLM is then approximated by (2) One can find in Lacaux (2004) theoretical rates of convergence for Y h,N . Theoretical results are given on any compact set K but the algorithms are written on the interval [0, 1] or on the cube [0, 1] 2 for sake of simplicity.

Simulation procedure
Let us first introduce the formula used for the simulation of one-dimensional RHMLM (d = 1). The parameters of the simulation are: • n: the number of terms in the sum; • k: the number of discretization points; • h(t) (t = 1, . . . , k): the value of the multifractional function at each discretization point; • m = ν(C): the mass value, usually fixed at 1. Actually we choose in the simulation to take ν the uniform distribution on the unit circle in C.
If we consider the random variables that are standard to simulate • θ n : n-vector generated according to a uniform distribution in [0, 1], Z n = e i2πθn ; • U n : n-vector generated according to a Bernoulli distribution with equiprobability for one-dimensional simulation; • R n : n-vector of cumulated sums of samples from exponential distribution with parameter λ = 1, and apply (2) When d = 2, the formula (3) becomes where k is the number of discretization points for each axis kx and ky, c 2 = π, (U 1 n , U 2 n ) are uniformly distributed on the unit sphere in R 2 .
In practice, the difficulty of the simulation lies on both the number of terms in the sum n and in the calculation of scalar products between k-vectors. In our cases, n may vary between 10 3 and 10 8 , and k between 100 and 1000, which may imply until 10 14 loops for the computations in the extreme case.

Using FracSim
FracSim is available from the Comprehensive R Archive Network (CRAN, http://CRAN. R-project.org/ or one of its mirrors). Instructions for package installation are given by typing help(INSTALL) or help(install.packages) in R.

Implementation of FracSim
FracSim is a set of R and C functions that allows performing RHFLMs and RHMLMs simulations. Two classes of functions are implemented: • R functions: fracsim.1d() and fracsim.2d() perform initialisations thanks to parameters set by the user (n, k and h(t) t=1,...,k ), call C functions for calculation then get back results for graphical representations.
• C functions: core_1D.c and core_2D.c perform the tasks in the computation that are time consuming because of the number of loops required.
The R environment is the only user interface. The R procedure calls a C subroutine, whose results are 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 command used to produce graphical output in Appendix A. In the following, R commands are preceded by the prompt symbol R>.

One-dimensional fractional motions
To simulate RHFLM, one needs to specify one real parameter 0 < h < 1 or, equivalently, the multifractional function is set to a constant h(t) = h, ∀t ∈ [0, 1].
The result of the function fracsim.1d is an R object of class list. It contains the following elements: • the vector simul.h of the simulated regularity value; • the vector t of discretization points; • the k-vector process whose elements are the RHFLM value at each discretization point.
Thus to produce the graphs on the Figure 1, one has to call the function plot as R> plot(x=X05$t,y=X05$process,type="l") As previously mentioned, we can see on the Figure 1, the greater the regularity, the smoother the sample paths.

One-dimensional multifractional motions
To simulate RHMLMs, one has to give as input for h, either a vector of length k, or a function that will calculate the regularity function at each discretization point. The two following commands give examples of: • an increasing regularity function from 0.1 to 0.9 (Hinc); • a sinusoidal regularity function oscillating according to a linear combination of sinus (Hsin); R> Hinc = seq(from=0.1,to=0.9,length=1000) R> Hsin = 0.25+0.25*sin(seq(0,1,length=1000)*(6*pi)) We can then call the function fracsim.1d() with the vectors previously defined given as parameter to the argument h: R> Xinc = fracsim.1d(h=Hinc,k=1000,n=5000) R> Xsin = fracsim.1d(h=Hsin,k=1000,n=5000) An equivalent way of simulating a process according to a sinusoidal regularity function consists in using a function as input for h. Let us define the function hsin as below To highlight the behaviour of the RHMLM according to the multifractional function, we represent both regularity functions and the corresponding sample paths on the Figure 2.
As previously mentioned for RHFLMs, the greater the regularity, the smoother the trajectory and this can be observed on one graph: for example when we look at the sinusoidal regularity, the corresponding process is much more perturbed when the regularity is close to zero.

Two-dimensional fractional motions
RHFLMs can be simulated by calling the function fracsim.2d() with one value given for h. For instance, R> X2d05 = fracsim.2d(h=0.5,kx=100,ky=100,n=100000) sets h constant and equal to 0.5 on a 100×100 regular grid with 10 5 terms in the sum. The user can simulate on an irregular grid by specifying kx and ky as vectors.
The elements of the object X2d05 provided by the function fracsim.2d() are the same as those given by fracsim.1d(), except that, in this case, the motion values are now located on a grid and are stored in a matrix object.
To represent two-dimensional motions, we now use the R function persp to draw a surface: R> persp(X2d05$process,shade=0.5,phi=30) Arguments shade and phi are optional; they can be set to give a nice representation: shade modify the shade at a surface facet and phi is the colatitude angle defining the viewing direction.
Other functions can be used to represent two-dimensional motions, we propose a comparison of various representations at the end of this section.

Two-dimensional multifractional motions
To simulate two-dimensional RHMLMs, one has to give for h, either a matrix which dimensions are kx × ky, or, as in the one-dimensional case, a function. This function takes two parameters as input and returns a matrix whose elements are calculated on the grid defined by kx and ky vectors. For example, we can define two matrices corresponding to the cases presented in one dimension with the following command: R> Hinc = matrix(rep(seq(0,1,l=100),100),nc=100) R> Hsin = matrix(rep(0.25+0.25*sin(seq(0,1,l=100)*(6*pi)),100),nc=100) One can obtain the graph of the regularity function ( Figure 4) as well by calling the function persp with the same optional arguments used for the representation of fractional twodimensional motions.
Once the matrix containing the regularity function at each point of the grid is built, we can perform the simulation of the motion by calling, as in the fractional case, the function fracsim.2d().
Then fracsim.2d can be called as below:    Figure 5), the surface looks like a landscape where we go from the mountain with several peaks (low regularity) to the plain (high regularity). With a sinusoidal regularity, the 'landscape' looks like an alternation of moutains and valleys.

Graphical representations
The way we represent a two-dimensional motion has an influence on the visual impression provided. If we only use the persp function to draw a surface, it could be interesting to modify point of view by changing phi (vertical rotation) and theta (horizontal rotation) parameters.
However, we can use other kinds of representations such as images or contour line plots as it is done on the Figure 6. Amongst these various representations, it is not easy to decide which is the best; it can depend on the practice of the user. However, let's note that with a high discretization and a regularity close to zero, the visual impression given by the image is smoother than the one given by the sheet (Figure 8, obtained thanks to improvements detailed in Section 4). This can be explained by the inability to distinguish varying intensity in the same color on the image whereas the peaks appear clearly on the sheet. On the other hand, the image of a two-dimensional motion highlights more precisely the global trends of the process, which is more difficult to observe on the sheet. For instance, in the case of a motion with an increasing regularity function along the diagonal, it would be difficult to set the point of view to give a good sheet representation whereas the behaviour of the motion appears clearly on an image ( are described at http://www.cict.fr/sys/Serveurs/ondine.html) for the configuration we performed, and we did not have any trouble with too long calculations.
For two-dimensional motions, the calculation durations order of magnitude can be linearly extrapolated from the basic configuration {kx = ky = 100, n = 10000}. Such a configuration, with 10 8 loops, lasts about 2 minutes on the machine previously mentioned. We can estimate the potential duration for the extreme configuration we would like to perform {kx = ky = 500, n = 10 8 } at about one year. In Section 4, we described out first steps in parallel computing to try to reach this extreme configuration in a 'reasonable' duration.

Parallel computing
In order to go further in our experimental design, we use a supercomputer via the scientific group Calmip (Calcul en Midi-Pyrénées http://www.calmip.cict.fr/). This machine is exclusively dedicated to intensive computation. Its is composed of 68 processors (1.5 GHz) with 136 Go of RAM memory in shared configuration. Technical characteristics are detailed in CALMIP (2005).
To make better use of this machine, we adapted our programs to exploit multi-threadings possibilities via parallel computing. We used an implicit parallelization (OpenMP 2005) that seemed to be a good compromise between automatic and explicit parallelization. The former is easier to implement (just add an argument to the compilation command) but less efficient. The latter, using for example Message Passing Interface, (MPI 2005), is more efficient but more demanding: the programmer has to give each processor a specific task and to ensure inter-processors communications. With OpenMP, one just has to add few lines to precise the piece of code that has to be shared between the processors and then performed in a multi-threading way.
To do so, we convert our function for two-dimensional simulations into one main C program. It is composed of three parts: • initializations and random number generation; • computations; • output to a file.
The first and the third part were previously dealt in R; the second one is equivalent to the C part of the package. This task, which contains the most demanding part of the computation, requires three nested loops: the rows, the columns and the terms of the sum to be caculated in each cell.
Our problem is quite easy to parallelize because of the independence between the calculation performed on each cell of the grid. Thus we can divide the grid in several blocks and give each processor a block to perform computation on. Technically, the loops intructions are nested in an OpenMP environment written in C as follows: #pragma omp parallel private(tmp,tmp1,tmp2,tmp3,l) { #pragma omp for for (i=0;i<kx;i++) /* loop over the rows */ #pragma omp parallel shared(i) { #pragma omp for for (j=0;j<ky;j++) /* loop over the columns */ { for (l=0;l<n;l++) /* loop over the terms of the sum */ { /* * Computations */ } /* End of loop l */ } /* End of loop j */ } /* End of loop i */ } /* End of parallel area */ OpenMP instructions must be preceded by the # character. Details on how to run such a program are given in Appendix B.
First simulations performed in this way are promising. We can divide the computation duration by at least the number of processors involved in the computations. We use from 8 to 16 processors. Furthermore, one processor of the supercomputer is faster than those of the current machine previously used.
So a simulation with kx = ky = 512 and n=10 7 last about 25 hours using (only) 8 processors instead of an estimated time around 2 months (1500 hours) with our previous configuration (R and C functions).

Conclusion
In this article we have proposed an efficient way to produce non Gaussian fields and we hope it can help to model various irregular phenomena. Of special interest is the use of parallel computers that shows that the algorithm is flexible enough. On the other hand, other non-Gaussian models exist and could be simulated in a similar way (Cohen, Lacaux, and Ledoux 2005) and we plan to generalize our techniques to those fields. Moreover statistical estimations of the Hölder exponent (Coeurjolly 2000) are possible for fractional Brownian motions and we plan to apply this to non-Gaussian fields.
We maintain a web page about RHFLMs and RHMLMs. It will always contain the latest release of the FracSim package: http://www.lsp.ups-tlse.fr/FracSim/.
In these calls, the argument ky is optional; if missing, its value is taken to be equal to kx.
• Graphical representations of 2D-RHFLMs (Figure 3) R> par(mfrow=c(1,3)) R> persp(X2d01$process,shade=0.5,phi=30) R> persp(X2d05$process,shade=0.5,phi=30) R> persp(X2d09$process,shade=0.5,phi=30) The function persp produces the plot on a grid determined by the size of the matrix given as the first parameter (X2d0*). We used shade=0.5 to define the shade at a surface facet in order to provide an approximation to daylight illumination as explained in the R help for the persp function. The option phi=30 implies a vertical rotation that sets the point of view above the surface.

B. Run parallel program
The compilation is done using the option -openmp to deal with the OpenMP command preceded by '#'.