Beyond multifractional Brownian motion: new stochastic models for geophysical modelling

Abstract. Multifractional Brownian motion (mBm) has proved to be a useful tool in various areas of geophysical modelling. Although a versatile model, mBm is of course not always an adequate one. We present in this work several other stochastic processes which could potentially be useful in geophysics. The first alternative type is that of self-regulating processes: these are models where the local regularity is a function of the amplitude, in contrast to mBm where it is tuned exogenously. We demonstrate the relevance of such models for digital elevation maps and for temperature records. We also briefly describe two other types of alternative processes, which are the counterparts of mBm and of self-regulating processes when the intensity of local jumps is considered in lieu of local regularity: multistable processes allow one to prescribe the local intensity of jumps in space/time, while this intensity is governed by the amplitude for self-stabilizing processes.


Introduction
Multifractional Brownian motion (mBm) has proved to be a useful model in geophysics, as witnessed by other articles in this special issue. Other references of interest include Gaci and Zaourar (2011), Keylock (2010), Wanliss (2005) and Wei et al. (2004). A more complete, but sill partial, list of works is compiled at the web page http://regularity.saclay.inria.fr/ theory/stochasticmodels/bibliombm, which also includes applications in other fields. While a versatile model, mBm is of course not always the most appropriate one in various situations: schematically, mBm is a good choice when one needs to model signals with varying local regularity (see Sect. 2 for a precise definition of local regularity), and when one has reasons to believe that local regularity has important implications in terms of, for instance, classification or detection. The use of mBm puts at least two restrictions on the studied signals: they must evolve in a continuous fashion, and their local regularity must be independent from the amplitude of the signal.
These assumptions are not verified for some geophysical signals. As we show below, it turns out that amplitude and local regularity are linked by a functional relation in certain cases. Such a feature requires the development of new models, which we call self-regulating processes. Also, it is clear that many geophysical records are not continuous. This calls for the use of alternatives to mBm that allow the presence of jumps and the possibility to tune their intensity in time, either in an exogenous way (similar to what mBm allows for local regularity), or through a functional relation with amplitude (similar to what self-regulating processes achieve again for local regularity).
The aim of this work is to present some stochastic processes obtained as generalizations of mBm which display either of the properties mentioned above, and could be of potential use in geophysics. The remainder of this article is organized as follows: in Sect. 2, we briefly recall some basic facts about fractional and multifractional Brownian motion. Section 3 explains how to build self-regulating processes and describes their application to natural terrain modelling and temperature records analysis. Multistable processes, i.e. discontinuous processes where the local intensity of jumps can be tuned exogenously, are presented in Sect. 4, and their endogenous version are the topics of Sect. 5.
We end this introduction with a word of warning: multifractal processes and multifractal analysis are popular and important topics in geophysics (Lovejoy and Schertzer, 2013). Although they sound the same, "multifractal" and "multifractional" refer to two very distinct concepts and approaches. To give but one example, multifractional Brownian motion, as presented below, does not typically have nontrivial multifractal properties.

Background on fBm and mBm
Let us start by recalling some basic facts about fractional Brownian motion (fBm) (Kolmogorov, 1940;Mandelbrot and Van Ness, 1968). This is a centred Gaussian process that essentially depends on a single parameter, usually denoted by H and called the Hurst exponent, which belongs to (0,1). Its covariance function R H reads where γ H is a positive constant. When H = 1 2 , fBm reduces to standard Brownian motion. Fractional Brownian motion exhibits features that make it a useful model in various fields such as geophysics, financial and teletraffic modelling, image analysis and synthesis, and more. These features include self-similarity, long-range dependence (when H > 1/2), and the ability to match any prescribed constant local regularity. Local regularity is a major theme in this work. We will measure it with the help of the pointwise Hölder exponent, which is defined as follows for a continuous nondifferentiable stochastic process X: one says that X belongs to C α (x 0 ), α ∈ (0, 1), if there exists ǫ > 0 and C ∈ R, such that The pointwise Hölder exponent of X at x 0 is with the convention that sup(∅) = 0. When there is no risk of confusion, we shall write β(x 0 ) in place of β X (x 0 ). This exponent can be roughly understood as follows: the largest absolute increments |X(x) − X(x 0 )] in the neighbourhood of x 0 are of the order of |x − x 0 | β(x 0 ) . The Hölder exponent of fBm is almost surely equal to H at all times. Thus, an fBm with small H will have paths looking rougher than the ones of an fBm with large H .
While fBm is a useful model, the fact that most of its properties are governed by the single number H restricts its application in some situations. For instance, long-range dependent fBm, which require H > 1 2 , must have smoother paths than Brownian motion, which corresponds to H = 1/2. Also, sample paths have almost surely everywhere the same regularity, a feature not typically observed in geophysical signals.
Multifractional Brownian motion was introduced in Peltier and Lévy Véhel (1995) to overcome these limitations. The basic idea is to replace the real H by a function t → h(t) ranging in (0,1). The construction of mBm is best understood through the introduction of a fractional Brownian field. Fix a positive real T . We define the fractional Brownian field B(t, H ) on [0, T ] × (0, 1) as the Gaussian field given by where W denotes a Gaussian white noise measure. For every H in (0,1), the process (B H t ) t∈[0,T ] := (B(t, H )) t∈[0,T ] is a fractional Brownian motion with Hurst parameter H , but when one deals with B(t, H ), H is just seen as an independent variable.
For a deterministic continuous function h : [0, T ] → (0, 1), we call multifractional Brownian motion with functional parameter h the Gaussian process B h (t) := B(t, h(t)). We say that h is the regularity function of the mBm.
Its covariance reads The main properties of mBm are as follows: the pointwise Hölder exponent at any point t of B (h) is almost surely equal to the minimum of h(t) and β h (t), where β h (t) is the pointwise Hölder exponent of h at t. For a smooth h, one thus can control the local regularity of the paths by the value of h. In addition, the increments of mBm display long-range dependence for all non-constant h(t). Finally, when h is C 1 , mBm is tangent to fBm with exponent h(u) in the neighbourhood of any u in the following sense (Falconer, 2003): In other words, mBm looks locally like fBm around each point. These properties show that mBm is a more versatile model than fBm: in particular, it is able to mimic in a more faithful way local properties of geophysical records, financial logs (Bianchi, 2005), biomedical signals (Betrouni et al., 2012) and Internet traffic (Li et al., 2011) by matching their local regularity. This is important, e.g. for purposes of classification, detection or real-time control. The price to pay is of course that one has to deal with the added complexity introduced by having a functional parameter instead of a single number. Various tools exist to deal with the simulation and calibration of mBm. Several of these are available in the Fra-cLab software toolbox (FracLab, 2012), which allows one to synthesize mBms and estimate the regularity function from sampled data. Two useful tools for estimating H (t) are the so-called generalized quadratic variations and incrementsratio statistics. Let us briefly recall the definition of the former. Assume sample points {X( p N ); p ∈ {0, . . . , N − 1}} of a discretized trajectory of the process {X(t)} t∈[0,1] have been observed. Set d 0 = 1, d 1 = −2, d 2 = 1, and fix t in [0, 1].
Then the generalized quadratic variation of X around t is the random quantity: where for a fixed γ ∈ (0, 1). Then a suitable renormalized and centred version of V (2) N (t) tends almost surely to H (t) when N tends to infinity. We refer the reader to Bardet and Surgailis (2013) for more details on estimation issues and a presentation of the increments-ratio statistics. Figure 1 displays an mBm with linear h function increasing from 0.1 to 0.8. The reader who wishes to see mBm "in action" should consult http://regularity.saclay.inria.fr/theory/ stochasticmodels/multifractional-brownian-motion for an animation showing the effect of a varying regularity function on the paths of the process.

Motivation
Modelling with mBm amounts to fixing a target regularity function h, and building a process whose pointwise Hölder exponent at each point t will be equal to h(t). The local regularity is thus set in an exogenous way, in the sense that h is prescribed in an independent manner.
In Echelard et al. (2010); , experimental findings were reported indicating that, for certain natural phenomena such as electrocardiograms or natural terrains, there seems to exist a link between the amplitude of the measurements and their pointwise regularity. This intriguing fact calls for the development of new models where the regularity would be obtained in an endogenous way: in other words, the Hölder exponent at each point would be a function of the value of the process at this point. With such models, one could, for instance, synthesize numerical terrains which would automatically be more irregular at high altitudes and smooth in valleys.
Stochastic processes with this property are generically termed self-regulating. Some instances have been studied in Barrière et al. (2012). They satisfy a functional relation of the form β X (t) = g(X(t)) for all t, where g is a smooth deterministic function defined on [0, 1] and ranging in a subset [a, b] of (0, 1). In the sequel, the interval [a, b] is fixed, and we shall denote k g = sup t∈[0,1] |g ′ (t)|.

Self-regulating multifractional process
One construction of a self-regulating process (srp) starts from a fractional field The intuitive idea is simple: start from a "line" on B, i.e. a plain fBm, that we will denote Z 0 . Since g(Z 0 ) is (almost surely) a continuous function ranging in [a, b], we can use it as a regularity function of an "mBm". In other words, we consider the process defined as Z 1 = B(t, g(Z 0 (t))). Similarly, we construct Z 2 = B(t, g(Z 1 (t))). Continuing this procedure, we obtain a sequence of processes Z n , for n ∈ N. With additional technical tricks, and under some conditions, it is possible to prove that this sequence converges, and that its limit Z is indeed self-regulating. This makes intuitive sense, since, at the limit, we should have Z = B(t, g(Z(t))), which exactly means that the pointwise Hölder exponent of Z at time t is equal to g(Z(t)).
Without going into full details, let us briefly state the main steps in the actual construction of Z. We will make use of a crucial fact about B, namely that it is smooth in the second variable H : more precisely, the function H → B (H, t) is almost surely continuously differentiable over [a, b]. Moreover, the random variable is almost surely finite. We will need the following notation: for X a continuous non-constant field defined as K, and α ′ , β ′ two real numbers, we set In words, X β ′ α ′ is just an affine rescaling of X so that it ranges between α ′ and β ′ .
Fix now two real numbers α < β and choose α ′ (ω), β ′ (ω) two random variables such that α ≤ α ′ (ω) < β ′ (ω) ≤ β and We define a stochastic operator α ′ ,β ′ almost surely on the set of continuous processes from [0, 1] to [α, β] as follows: One can show that α ′ ,β ′ possesses a unique fixed point, which we shall denote Z g . By definition, Z g verifies . (3) In addition, Z g is indeed self-regulating, i.e. for all t, the pointwise Hölder exponent of Z g is almost surely equal to g(Z g (t)). Note that, by construction, Z g ranges in [α, β]. By choosing adequately the interval of definition of g, one can control the values taken by the process. With some additional work, one can also force Z to be close to a prescribed smooth shape. The reader who wishes to see the effect of prescribing shapes on srps will find an animation at http://regularity.saclay.inria. fr/theory/stochasticmodels/self-regulating-processes.
For simplicity, we have presented the construction of the self-regulating multifractional process in one dimension. However, the theory applies in higher dimensions without essential modifications. Figure 2 shows a bi-dimensional srp that mimics a natural terrain with g(Z) = (1 − Z) 2 . Notice how, due to the shape of g, the synthetic terrain automatically adjusts its regularity to be smoother at lower altitudes and rougher at higher ones.

Midpoint-displacement self-regulating process
A quite different construction of an srp is based on a waveletlike decomposition. Recall the definition of the "triangle" function: Fig. 3).
Let Z j k be i.i.d. random variables following an N (0, 1) law. It is well known that is a representation of Brownian bridge. We remark that the factor 1/2 in the expression 2 −j/2 ϕ j k Z j k corresponds to the constant pointwise Hölder exponent of the process. Heuristically, the term 2 −j/2 ϕ j k entails a variation with amplitude 2 −j/2 and duration 2 −j . In other words, variations of time length h = 2 −j are of the order of h 1/2 . It is easy to prove that the modified process j,k 2 −jβ ϕ j k (t)Z j k has almost surely everywhere Hölder exponent β for β ∈ (0, 1). One can take advantage of this fact to build, in an iterative way, a selfregulating process. As above, g is a C 1 function ranging in [a, b] ⊂ (0, 1).
As an initialization, set X −1 ≡ 0 and define the sequence of processes (X j ) j ∈N by One can show that, almost surely, the sequence (X j ) j ∈N converges uniformly to a continuous process X, which we shall call self-regulating random midpoint displacement process (srmdp). This process verifies the following property: almost surely, for all t, β X (t) = g(X(t)). In other words, X it is indeed self-regulating. An example of an srmdp with g(x) = 1 1+5x 2 is displayed in Fig. 4. One verifies that, due to the form of g, the process is smoother when it is close to 0 and gets rougher as its amplitude becomes large in absolute value.

Statistical estimation
Checking whether empirical data are well modelled by an srp amounts to verifying whether there is a functional relation between the amplitude and the regularity at all times/locations. In practice, this can be done by visually inspecting the phase plane (Z, β Z ). If the points in this plane organize themselves along well-defined curves, this is a sign of self-regulation. If they appear randomly scattered, then a self-regulating model is not relevant. We show in Fig. 5 such a phase plane plot for the self-regulating process Z displayed in Fig. 4. One can see how the points roughly follow the curve g(Z) = 1 1+5z 2 . In Fig. 6, we plot β Z (t) as a function of Z(σ (t)), where σ is a random permutation. Obviously, the self-regulating property is lost through shuffling, and this is a typical graph for a non-self-regulating process.   Fig. 4 and theoretical self-regulating function. The blue circles are the points (Z(t), β(t)) and the solid black line is the theoretical g function.
In many practical situations, however, the plain phase plane is not sufficient to detect self-regulation because a few outliers can easily blur the relation between amplitude and exponents. A more quantitative display is needed. A clearer picture is obtained by regularly binning (or "histogramming") the phase plane and displaying in each bin the number of couples (X(t), β X (t)) that fall within it. Figure 11, which shows  . 6. Phase plane β Z (t) as a function of a random shuffling of Z(t) for the process shown in Fig. 4: each blue circle is a point (Z(σ (t)), β(t)) where σ is a random permutation. both the plain and the binned phase plane for temperature data, illustrates this.
When visual inspection confirms the presence of selfregulation, the self-regulating function g can be estimated in two ways: a non-parametric one, which essentially refines the image obtained in the binned phase plane, and a www.nonlin-processes-geophys.net/20/643/2013/ Nonlin. Processes Geophys., 20, 643-655, 2013 parametric one in the case of a midpoint displacement selfregulating process. Let us briefly explain this second method: let N = 2 j be the number of samples and choose a sequence (ε N ) N such that as N → ∞. Fix x in the observed range of X. Let s 1 , . . . , s n j denote the real numbers of the form (k + 1/2)2 −j such that X j −1 (s i ) ∈ [x − ε N , x + ε N ] and k 1 , . . . , k n j denote those integers k such that s i = (k i + 1/2)2 −j . Set finally T j,n j = n j i=1 X j (s i ) − X j −1 (s i ) 2 and definê g j (x) := log 2 (n j ) 2j − log 2 (T j,n j ) 2j .
Thenĝ j (x) converges almost surely to g(x) when j tends to infinity. We refer the reader to  for more details on these statistical aspects. Both of these estimation procedures, as well as synthesis methods for self-regulating midpoint displacement processes and self-regulating multifractional processes, are available in FracLab. An example of a non-parametric estimation using the process shown in Fig. 4 is displayed in Fig. 7.

Application to the characterization of natural terrains
As a first application of srps in geophysics, we briefly recall in this section the findings of Echelard et al. (2010). We have studied four mountainous zones: two younger ones (Himalaya and Rocky Mountains), and two older ones (Massif Central and Tibesti). The data were obtained from the internet site of the United States Geological Survey. The resolution is 3 arc seconds (approx. 90 m). They have the following dimensions: Tibesti: 6097 × 6393; Massif Central: 3807×3835; Rocky Mountains: 9274×6072; Himalaya: 5620 × 6767. In order to check for a possible self-regulating relation, we represent, as explained above, the "altitude-exponent" plane. As one can see in Fig. 8 in the case of the Rocky Mountains, the results are not convincing. Indeed, scattered exponents with the same altitude are observed. The points do not fit the graph of a function. Rather, they form a roughly elliptic blob. The same kind of results were obtained with the other data. This tends to indicate that, at this scale, terrains are not homogeneous enough to be represented by self-regulating processes.
The situation is, however, different when one studies smaller regions. We have made experiments with windows of size 512 2 pixels. To obtain statistically significant results, we have moved such windows by steps of 50 pixels for the four mountainous zones. Figure 9 displays a typical example of the results we have obtained: we have found that, consistently for most windows, a clear relationship between altitude and exponent is indeed present. For the Himalaya and Rocky Mountains, the estimated g function decreases, while it increases for Massif Central and Tibesti, as can be seen in Fig. 9. A possible explanation of this difference is that it might be due to erosion: high-altitude points in old mountains are smoothed out by erosion, and thus high altitudes translates into more smoothness in this case. In contrast, high altitudes zones in young mountains are typically more irregular than low altitude ones. Further investigations are needed to confirm the explanation for these findings.

Temperature records
We now analyse daily temperature records, and show that they also typically satisfy a self-regulating property. Our data come from the CDIAC facility (see Menne et al. (2011) for details). We studied minimum, mean, and maximum temperatures at more than 100 sites scattered over the USA. Figure 10 shows an example of a temperature log at the Corning (Arkansas, USA) station along with the estimated regularities. It is apparent that both graphs vary approximately in phase. To quantify this impression, we display in Fig. 11 a comparison between a plain phase plane, a binned one, and a binned one with randomly shuffled data for the minimum temperatures at Greensboro, Alabama. One clearly sees in the middle graph that the time evolution of temperatures at this site exhibit a self-regulating property, with the variations of temperature being more irregular at low temperatures and vice-versa. Figure 12 shows the binned phase planes for the maximum temperatures at Lees Ferry, Arizona, and mean temperatures at Corning, Arkansas, and Santa Cruz, California. These illustrate the three kinds of shapes that we have observed in all our experiments. By far the most frequent is the top one: roughly 75 % of the sites we have analysed exhibit this shape, which is the sign that, for such records, the regularity of the time evolution of temperatures increases when temperature increases. As can be seen in Fig. 13, which displays a sample of 42 binned phase planes among the ones we have computed, this is particularly apparent for the sites at Ardmore, Oklahoma, Conway, South Carolina, and Dickson, Tennessee. Around 15 % of the records exhibit a shape reminiscent of a crescent in the phase plane. This is the case at sites like Lees Ferry, Arizona (Fig. 12), and Tooele, Utah (Fig. 13). This means that, at these sites, temperatures vary in a more regular way when they are extreme, and less smoothly when they are close to their mean. For the remaining 10 % sites or so, the self-regulating property is not present. This is often due to the fact that, except for a small proportion of days, the temperatures remain within a rather small interval: this is, for instance, the case for Santa Cruz, California (Fig. 12), and Newport, Oregon (Fig. 13). As a result, the relation amplitude-regularity does not manifest itself since there is not enough variability in the records. The full characteristics of the studied sites are given in the Appendix. Geophysical explanations of our findings remain to be found; we hope to have raised the interest of researchers in this field,

Massif central
Rocky mountains Tibesti Himalaya who are able understand where the self-regulating property of temperature records comes from, and how it can be used for further analysis of such data.

Multistable processes
Multifractional Brownian motion and self-regulating processes allow one to model a wide spectrum of irregular natural phenomena, provided they do not display jumps. Indeed, almost surely, the paths of the processes studied in the previous sections are continuous. In view of analysing discontinuous phenomena that appear in various areas of geophysics, other classes of stochastic processes must be considered. The simplest processes of this kind, which are the counterpart of Brownian motion, are called stable motions. The name comes from the fact that, similarly to Brownian motion, they enjoy some invariance property in law under linear combination. The interested reader should consult Samorodnitsky and Taqqu (1994) for a thorough description of these processes. Stable motions are characterized by various parameters, the most important of which is the stability exponent, usually denoted α, which ranges in (0, 2). This exponent measures the intensity of jumps. Rather than explaining the precise  mathematical meaning of this statement, we display in Fig. 14 four stable motions with increasing value of α. We believe these graphs make rather intuitive what is meant by "intensity of jumps", with the first graph displaying large jumps while the last showing very small ones. Very much in the same way that, for fBm, the pointwise regularity is constant, the intensity of jumps remains the same all along a trajectory of a stable motion. Using the same approach that generalizes fBm to mBm by letting H evolve in time, we can extend stable motions into multistable ones by letting the stability exponent α vary in time. Again, we will not go into mathematical details, for which the interested reader should consult Le Guével and Lévy Véhel (2012). Rather, we display in Fig. 15  synthesize and estimate multistable processes are available in FracLab.

Self-stabilizing processes
In the Gaussian frame, we have moved from fBm to mBm to self-regulating processes. Similarly, in the stable case, the next step after having generalized from stable to multistable motions is to define self-stabilizing processes, i.e. processes X verifying α X (t) = g(X(t)) at all points, for some smooth function g ranging in (0, 2), and where α X (t) is the local jump intensity exponent of X. We do not go into mathematical details, but mention that self-stabilizing processes are obtained in the same way as self-regulating ones: the first method is a fixed point approach, which starts from a field of stable motions and iterates an operator similar to α ′ ,β ′ defined in Sect. 3.2.1. The second one is in the spirit of the iterative construction of Sect. 3.2.2, replacing the normal random variables Z j k with stable ones having adequately chosen α index. Figure 16 displays two self-stabilizing processes with respectively linear and sine prescribed shapes obtained with FracLab. Top: a self-stabilizing process with prescribed linear shape, which has larger intensity of jumps when it is small. Bottom: a selfstabilizing process with prescribed sine shape, which has larger intensity of jumps when it is close to 0.