Stochastic models for fractal processes

Journal of Statistical Planning and Inference 80 (1999) 123–135

www.elsevier.com/locate/jspi

Stochastic models for fractal processes 1 a Centre

V.V. Anha; ∗ , C.C. Heydeb , Q. Tienga

in Statistical Science and Industrial Mathematics, Queensland University of Technology, GPO Box 2434, Brisbane, Q 4001, Australia b School of Mathematical Sciences, Australian National University, Canberra, ACT 0200, Australia Received 15 December 1997; accepted 22 June 1998

Abstract This paper considers the situation where a stochastic process may display both long-range dependence (LRD) and intermittency. The existence of such a process is established in Anh et al. (1999). Existing works have commonly paid attention either to LRD or intermittency quite separately. This paper oers a convenient framework to study both eects simultaneously. A method is given to estimate and separate the two eects. The wavelet theory plays an essential role in this procedure. Numerical experiments on fractional Brownian motion and multiplicative c 1999 Elsevier Science B.V. All rights cascade processes con rm the power of the method. reserved.

1. Introduction A commonly-used model for the dynamics of turbulent processes is the Itˆo stochastic dierential equation (SDE), which generates a Markov process X (t); t¿0; with values in Rn ; n¿1 (Stroock and Varadhan, 1979). Higher-dimensional SDEs generating Markov random elds (MRF) in the sense of P. Levy have also been investigated (see, e.g., Gikman (1976), Anh et al. (1997)). The conditions for the Markov property of random elds arising as solutions of linear stochastic partial dierential equations (SPDE) are detailed in Rozanov (1979, 1982). Recent studies have indicated that data in a large number of elds display long-range dependence (LRD) (Beran, 1992; Peters, 1994; Anh and Lunney, 1995). A fundamental example is fractional Brownian motion (fBm) BH with Hurst index H; 1=2 ¡ H ¡ 1. In addition to LRD, it has been found that many processes in nance and 2-D turbulence in particular exhibit a high degree of intermittency. A typical example is a multiplicative cascade process (Davis et al., 1994). A need therefore arises for the development of ecient methods and models to deal with both LRD and intermittency in the data. 1

Partially supported by the Australian Research Council grant A69531724. Corresponding author. Tel.: +61-7-3864-2317; fax: +61-7-3864-2310. E-mail address: [email protected] (V.V. Anh)

∗

c 1999 Elsevier Science B.V. All rights reserved. 0378-3758/99/$ - see front matter PII: S 0 3 7 8 - 3 7 5 8 ( 9 8 ) 0 0 2 4 6 - 8

124

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

Traditional models of Markov diusion, such as those exempli ed by SDEs, SPDEs and MRFs mentioned above, are not directed at generating LRD and intermittency. The theory of fractional generalised random elds of Ruiz-Medina et al. (1997) oers a suitable framework to study these phenomena. In fact, Ruiz-Medina et al. (1997) introduced the new concepts of -generalised random elds (-GRF) and -duality based on the theory of potential spaces L2 = I (L2 (Rn ));  ∈ R; where I = (I −
)−=2 is the Bessel potential of fractional order ;
being the Laplace operator (Stein (1970), pp. 134 –138). Random elds with LRD and second-order intermittency can be constructed and studied from -GRFs by exploiting the relationship between the Riesz potential and the Bessel potential. A key result of Ruiz-Medina et al. (1997) is the factorisation of the covariance operators of the -GRF X and its dual eld X˜ ; which leads to the abstract representation ˜ = (I−˜()); X˜ (L)

 ∈ [H −˜(T )]∗

(1.1)

for X˜ ; where T is a C ∞ -bounded domain of Rn ; L˜ is a continuous linear operator, (·) is a generalised white noise on L2 (T ); ˜ is the minimum fractional duality order of X ; and H −˜(T ) is a fractional Sobolev space. The stochastic Eq. (1.1) of fractional diusion extends the SPDE of Markov diusion. The covariance operator of X˜ de ned by Eq. (1.1) is −1 −1 R˜ = [I˜L˜ ]∗ [I˜L˜ ]:

(1.2)

A fundamental subclass of -GRFs is obtained by considering ˜ =  + ; ; ∈ Q, and −1 by selecting the operator L˜ of (1.2) to be the composition I− J ∗ ; where J ∗ is the adjoint of the Riesz potential J = (−
)− =2 ; 0 ¡ ¡ n: That is, R˜ = J I∗ I J ∗ :

(1.3)

The spectral density of the random elds with covariance function de ned by (1.3) is then given by f(!) =

1 1 1 ; (2)n |!|2 (1 + |!|2 )

0 ¡ ¡ n;

¿0; ! ∈ Rn ;

(1.4)

(see Anh et al. (1999)). Since this function involves the Fourier transforms of the Riesz kernel and the Bessel kernel, these random elds are named fractional Riesz–Bessel −1 motion (fRBm). It is shown in Anh et al. (1999) that the selection L˜ = I− J ∗ satis es the -duality condition for 0 ¡ ¡ n; thus establishing the existence of fRBm. This paper will consider stochastic processes (i.e. n = 1) with stationary increments and having spectral density of the form (1.4). It will be seen in Section 2 that the exponent determines the self-similarity (SS)=LRD of the process, while the exponent  indicates its intermittency. It is noted that fBm has spectral density of the form (1.4) with  = 0 (Yaglom (1986), p. 407). The form (1.4) means that fRBm may exhibit both SS=LRD and intermittency ( ¿ 0). In other words, both eects can be studied simultaneously on the basis of model (1.4). It is the purpose of this paper to develop a method to estimate  and of Eq. (1.4). Existing works commonly considered the estimation of the LRD exponent alone based

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

125

on a Geweke–Porter-Hudak-type procedure (Geweke and Porter-Hudak, 1983). We shall estimate  and simultaneously using the wavelet transform of stationary increment processes. In our approach, the Haar wavelet plays a prominent role. The estimation method is detailed in Section 3. A simulation experiment is reported in Section 4 on the performance of this method. 2. Wavelet transform of stochastic processes with stationary increments Let {X (t); t ∈ R} be a zero-mean continuous (in mean square) process. We shall assume that X (t) has stationary increments with spectral representation E(X (t + 1 ) − X (t))(X (t + 2 ) − X (t)) = D(1 ; 2 ) Z ∞ (ei1 ! − 1)(e−i2 ! − 1) dF(!); = −∞

(2.1) where F(!) is a left-continuous real-valued non-decreasing function on R such that Z ∞ Z  2 ! dF(!) + dF(!) ¡ ∞; ∀ ¿ 0 (2.2) 0



(Yaglom (1986), pp. 400 – 401). Assuming that the spectral distribution function F(!) of (2.1) possesses a derivative f(!) given by the form (1.4), the requirement (2.2) is satis ed for ¡ 3=2 and  + ¿ 1=2: It is noted that, under the assumption of X (t) being nonstationary but having stationary increments, the range ¡ 1 of Eq. (1.4) is extended to ¡ 3=2 (with 0 ¡ ¡ 1=2 holding only when X (t) is stationary). Consequently, we shall consider a stationary increment process X (t) with structure function (2.1) and spectral density 1 c ; c ¿ 0; 1=2 ¡ ¡ 3=2; ¿0; ! ∈ R: (2.3) f(!) = |!|2 (1 + !2 ) The variance of its increments, D(; ); has the following behaviour. Lemma 2.1. For  ¿ 0 suciently small and  + ¡ 3=2; c ( + − 1=2) 2(+ )−1  : D(; ) ≈ 22(+ )−3

(2.4)

Proof. In view of Eqs. (2.1) and (2.3), Z ∞ sin2 (!=2) d! D(; ) = 8c |!|2 (1 + !2 ) 0  (1=2 − ) ( + − 1=2) = 8c 4 ()    2 × 1 − 2 F3 {1; 1=2 − }; {3=2 −  − ; 1=2; 1}; 4   ( + − 1=2) 2(+ )−1 2  F + {1; }; {1;  +

+ 1=2;  +

}; 2 3 2 22(+ ) 4

126

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

with ¡ 3=2 and  + ¿ 1=2; 2 F3 being a generalised hypergeometric function (see Formula 3, p. 379, of Prudnikov et al., 1986)  ( + − 1=2) 2(+ )−1  = 8c 22(+ )  1=2 − 2 1 (1=2 − ) ( + − 1=2) × − ::: ; (2.5) − 2 () 3=2 −  − 4 where the remaining terms in the above expansion have powers larger than 2 in : Consequently, for  suciently small and 1=2 ¡  + ¡ 3=2; the right-hand side of Eq. (2.5) behaves as [c ( + − 1=2)=22(+ )−3 ]2(+ )−1 : As given by Eq. (2.4), the process X (t) is locally self-similar (in the sense that it is self-similar for  suciently small). It is also noted that as |!| → 0, the spectral density (2.3) behaves as |!|−2 : For 1=2 ¡ ¡ 1; the process possesses short-range dependence (SRD). It displays LRD for 1 ¡ ¡ 3=2: On the other hand, the spectral density (2.3) behaves as |!|−2( +) for suciently large |!|. We will next show that, for a xed , the dimension of the level sets of almost all sample paths of X (t) decreases to 0 as  increases to 3=2 − : In other words, the essential support of the high values of X (t) is a fractal set with dimension D ↓ 0 as  ↑ (3=2 − ): As a result, the component (1 + !2 )− of Eq. (2.3) gives a description of the (second-order) intermittency of the process and the exponent  indicates the degree of this intermittency. Lemma 2.2. Let X (t); t¿0; be a zero-mean Gaussian process with stationary increments. If its spectral densityf(!) satisÿes the conditions |!|2 +1 f(!)¿C1

(2.6)

for some C1 ¿ 0; 0 ¡ ¡ 1 and for all suciently large |!|; and D(; )6C2 ||2

∀ ¿ 0

(2.7)

is satis ed for some C2 ¿ 0 and some 0 ¡  ¡ 1; then dim{t: 06t6T; X (t) = x}61 − 

(2.8)

for all x; almost surely. Proof. See Theorem 7.1 of Berman (1972). The condition  + ¡ 3=2 of Lemma 2.1 implies that the condition (2.6) is satis ed for f(!) of (2.3) with ↑ 1: On the other hand, except for the rst term, all the other terms in the expansion (2.5) are negative for ¿ 1=2;  + ¡ 3=2: Thus, c ( + − 1=2) 2(+ )−1  ; ∀ ¿ 0: 22(+ )−3 Hence, condition (2.7) is satis ed with =+ −1=2; and as a result, D=dim{t; X (t)= x}63=2 − − : It is seen that D ↓ 0 as  ↑ (3=2 − ) as to be proved. D(; )6

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

127

We now obtain the spectral density of the wavelet transforms of stationary increment processes. Let (t) ∈ L1 (R) with its Fourier transform b satisfying Z

b (0) = 0;

∞

0

The function is de ned as

| b (!)|2 d! ¡ ∞: !

(2.9)

is then a wavelet function. The continuous wavelet transform of X (t) Z

−1=2

Wa (t) = a

= a1=2

Z



∞

X (s)

−∞ ∞

−∞

s−t a

 ds

X (au + t) (u) du;

a ¿ 0; t ∈ R:

(2.10)

Now, using the spectral representation (2.1) and Fubini’s theorem, the covariance function of Wa (t) at a xed scale a is given by Z ∞Z ∞ E[X (au + t + )X (av + t)] E[Wa (t + )Wa (t)] = a −∞

−∞

× (u) (v) du dv  Z ∞ Z ∞ Z ∞ i(au+)! −iav! (e − 1)(e − 1)f(!) d! =a −∞

−∞

−∞

× (u) (v) du dv Z Z ∞ Z ∞ (ei(au+)! − 1) (u) du =a −∞

−∞

×f(!) d! Z ∞ e−i! | b (a!)|2 f(!) d! =a −∞

∞

−∞

(e−iav! − 1) (v)dv



(2.11)

using (2.9), which is independent of t. Consequently, the process Wa (t) is stationary with spectral density fWa (!) = a| b (a!)|2 f(!);

a ¿ 0; ! ∈ R:

(2.12)

3. Parameter estimation Assuming a spectral density of the form (2.3), it follows from Eq. (2.11) that the variance of the wavelet transform Wa (t) at a xed scale a is Z d! ; a¿0 R(0; a) = ca | b (a!)|2 2 |!| (1 + !2 ) R Z d! : (3.1) = ca2 | b (!)|2 2 |!| (1 + a−2 !2 ) R

128

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

Our method is based on the Haar wavelet function:  1; 06t ¡ 1=2;   2   sin!=4 i! (t) = −1; 1=26t ¡ 1; b (t) = e−i!=2 :  4 !=4   0; elsewhere; Noting that, for the Haar wavelet,  2 !2 sin!=4 1 | b (!)|2 6 ; 2 −2 2  |!| (1 + a ! ) 16 !=4 |!|2

(3.2)

the right-hand side of Eq. (3.2) is then independent of a and is integrable on R since 1=2 ¡ ¡ 3=2: Hence the dominated convergence theorem yields Z Z d! d! | b (!)|2 2 = | b (!)|2 2 lim a→∞ R |!| (1 + a−2 !2 ) |!| R = M1 ¡ ∞: Consequently, as a → ∞; R(0; a) ∼ cM1 a2 :

(3.3)

This result suggests the linear regression ln R(0; a) = c1 + 2 ln a + u;

(3.4)

u being white noise and a suciently large, for an estimation of . Eq. (3.1) can also be written as Z d! 2( +) | b (!)|2 2 2 : R(0; a) = ca |!| (a + ! 2 ) R Another application of the dominated convergence theorem gives, for 1=2 ¡ + ¡ 3=2; R(0; a) ∼ cM2 a2( +) as a → 0; R where M2 = R | b (!)|2 d!=|!|2( +) ¡ ∞: The result (3.5) suggests the linear regression: ln R(0; a) = c2 + 2( + )ln a + v;

(3.5)

(3.6)

v being white noise and a near 0, for an estimation of + : We shall estimate R(0; a) via the Allan variance, which is de ned as 1 VaA = E(Xa (k + 1) − Xa (k))2 ; k ∈ R; 2 R (k+1)a X (t) dt (Allan, 1966; Masry, 1991). where Xa (k) = (1=a) ka For a ¿ 0; k ∈ R; t ∈ R; denote a; k (t) = a−1=2 (a−1 t − k); −1=2 a; k (t) = a

(a−1 t − k);

where (t) = X[0;1] (t); the Haar scaling function, and (t) = X[0;1=2] (t) − X[1=2;1] (t); the Haar wavelet function, A (t) being the indicator function.

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

De ne

129

Z

ca (k) =

R

X (s)a; k (s) ds;

(3.7)

Z Wa (k) =

R

X (s)

a; k (s) ds:

(3.8)

Note the slight change of notation from Eq. (2.10) to Eq. (3.8), where k = t=a: This is more convenient in the following development. Lemma 3.1. VaA =

a−1 E(ca (k + 1) − ca (k))2 2

(3.9)

= a−1 E(W2a (k=2))2 :

(3.10)

Proof. From the de nition of the Haar scaling function, Z (k+1)a X (t) dt = a1=2 Xa (k); ca (k) = a−1=2 ka

from which Eq. (3.9) follows directly. Also, from the de nition of the Haar wavelet function, Z (k+1=2)a Z (k+1)a ! X (t) dt − Wa (k) = a−1=2 ka

=2

−1=2

(k+1=2)a

(ca=2 (2k) − ca=2 (2k + 1));

or equivalently, ca (k) − ca (k + 1) = 21=2 W2a (k=2):

(3.11)

The result (3.10) now follows from Eqs. (3.9) and (3.11). Table 1 Estimation of fBm

b

b + b

b

1.00

0.991 (0.011) 1.041 (0.010) 1.190 (0.009) 1.283 (0.008) 1.370 (0.008) 1.414 (0.007)

0.993 (0.001) 1.042 (0.001) 1.190 (0.001) 1.288 (0.002) 1.379 (0.003) 1.423 (0.003)

0.002

1.05 1.20 1.30 1.40 1.45

0.001 0.000 0.005 0.009 0.009

130

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

Table 2 Estimation of nonstationary multiplicative cascades

H

+ 2H

2b

2(b

+b )

b

0.14

0.5

1.14

0.6

1.34

0.14

0.7

1.54

0.14

0.8

1.74

0.14

0.9

1.94

1.624 (0.051) 1.824 (0.045) 2.019 (0.038) 2.209 (0.031) 2.392 (0.024)

0.057

0.14

1.510 (0.031) 1.706 (0.028) 1.901 (0.025) 2.092 (0.022) 2.275 (0.019)

0.44

0.5

1.44

0.6

1.64

0.44

0.7

1.84

0.44

0.8

2.04

0.44

0.9

2.24

1.729 (0.036) 1.927 (0.031) 2.120 (0.025) 2.306 (0.020) 2.483 (0.014)

0.063

0.44

1.604 (0.025) 1.799 (0.023) 1.992 (0.020) 2.180 (0.018) 2.358 (0.015)

0.64

0.5

1.64

0.6

1.84

0.64

0.7

2.04

0.64

0.8

2.24

0.64

0.9

2.44

1.828 (0.022) 2.024 (0.019) 2.214 (0.015) 2.396 (0.011) 2.565 (0.008)

0.067

0.64

1.695 (0.020) 1.889 (0.018) 2.079 (0.016) 2.263 (0.013) 2.435 (0.011)

0.94

0.5

1.94

0.6

2.14

0.94

0.7

2.34

0.94

0.8

2.54

0.94

0.9

2.74

2.053 (0.003) 2.244 (0.002) 2.424 (0.002) 2.591 (0.001) 2.735 (0.001)

0.068

0.94

1.916 (0.008) 2.106 (0.007) 2.288 (0.006) 2.457 (0.006) 2.606 (0.005)

0.059 0.059 0.058 0.058

0.064 0.064 0.063 0.063

0.067 0.067 0.067 0.065

0.069 0.068 0.067 0.065

Recall that the process {Wa (k)} is stationary in k. Lemma 3.1 implies that the variance of this process at scale 2a can be estimated as K−1 P b 2a) = aV baA = a (Xa (k + 1) − Xa (k))2 : R(0; K k=1

(3.12)

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

131

Fig. 1. A time series of 4015 observations of the US–Singapore exchange rate.

Lemma 3.2. Under certain conditions; baA (K)) = VaA ; E(V baA (K) − VaA )2 = o E(V



1 K 1−

 ;

0 ¡  ¡ 1;

(3.13)

i.e. VaA is mean-square consistent as K → ∞ with convergence rate 1=K: Proof. See Masry (1991), where the conditions for (3.13) to hold are speci ed. Summarising, we estimate from the regression (3.4) and +  from the regression (3.6), where R(0; a) is estimated from Eq. (3.12). The unbiasedness and mean-square consistency of the estimators are indicated by Lemma 3.2. 4. Simulation experiments In the rst experiment, we generate data which possess LRD but no intermittency, that is, ¿ 1 and  = 0: Five hundred sample paths of fBm BH (t); 06t6T = 8192, for H = 0:50; 0:55; 0:70; 0:80; 0:90 and 0:95 were generated using Mandelbrot’s mid-point

132

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

Fig. 2. A sample path of fractional Brownian motion generated from Mandelbrot’s mid-point displacement algorithm with Hurst index H = 0:65:

displacement algorithm: Put BH (0) √ = 0; BH (T ) = 1; generate a zero-mean Gaussian noise with standard deviation = 21−H − 1 and add it to [BH (0)+BH (T )]=2 to obtain BH (T=2); divide  by 2H and repeat between t = 0; T=2 for BH (T=4), and t = T=2; T for BH (3T=4); proceed similarly to smaller and smaller scales (e.g., 10 levels). The spectral density of fBm has the form c=|!|2H +1 ; thus, in view of (2.3), 2H + 1 = 2 ; b a) for 128 scale i.e. = H + 1=2: From T = 8192 sample points, we computed R(0; levels for each sample path. The parameter was estimated from (3.4) using the last 104 scale levels (from 25 to 128) and +  was estimated from (3.6) using the rst 23 scale levels (from 3 to 25). The estimation was repeated over 500 sample paths. The average estimates are reported in Table 1, with the standard errors in brackets below each estimate. It is seen that the method works very well, with estimates for both and  quite close to the given values. In the second experiment, data which display both LRD and intermittency were generated from multiplicative cascade models (Davis et al., 1994) in combination with the mid-point displacement algorithm. The multiplicative stationary cascades X (t) are constructed by initially setting X (t) = 1 on [0; T ] then subdividing this segment into two equal parts and computing w1 X (t) and w2 X (t) on each part, where w1 and w2 are unit-mean nonnegative lognormal random variables. This procedure is repeated for

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

133

Fig. 3. The high-pass ltered series of the exchange rate series (top) and the generated fBm series (bottom).

a large number of times. As noted in Davis et al. (1994), these stationary cascades have spectral density of the form |!|− with = 1 − log2 E(w2 ). This class of models was rst introduced to describe the eects of intermittency in the dissipation eld of fully developed turbulence. Two examples were given in Fig. 1a and b of Davis et al. (1994) for = 0:94 and = 0:77: Nonstationarity and LRD are then introduced into the data by applying the above procedure (of generating stationary multiplicative cascades) on sample paths generated from the mid-point displacement algorithm. The resulting process, which displays both LRD with index H (as in fBm) and intermittency with index (as in multiplicative cascades), has spectral density of the form f(!) ˙

1 : |!|2H +

(4.1)

In view of the forms (2.3) and (4.1), we get + 2H = 2 : Five hundred sample paths of length T = 16384 were generated; the Allan variance was computed for 150 scale levels for use in the regressions (3.4) and (3.6). The parameter was estimated using levels 20 –150, while +  was estimated using levels 3–15. The average estimates for

and  (over 500 sample paths) in a number of combinations of and H are reported in Table 2. The signi cant values of b  indicate intermittency in the data. The estimates for the LRD exponent + 2H show an upward bias for small (e.g. = 0:14), but the estimates are quite accurate for large (e.g. ¿ 0:5).

134

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

Fig. 4. The logarithm of Allan variance against log a of the wavelet transform of the exchange rate series.

Finally, we consider the time series of daily US–Singapore exchange rates from June 1983 through June 1994, totalling 4015 observations. The time series is plotted in Fig. 1. For a comparison, we also plot in Fig. 2 a time series of fBm with H = 0:65 generated from Mandelbrot’s mid-point displacement algorithm. While the behaviour of the two series appears remarkably similar, the plots of their high-pass ltered series in Fig. 3 show that the exchange rate series displays more intermittency. The plot of the logarithm of Allan variance of its wavelet transform in Fig. 4 shows a steeper slope in the range of small scales, meaning a signi cant positive value for  and hence indicating signi cant intermittency in the data. The estimates from our method are as follows. b 0:949 (0:020)

[ + 1:075 (0:024)

b  0:150

where b was obtained using scale levels 26 –128 and +  was estimated using scale levels 3–26. The estimate b close to 1 and a signi cant value of b  indicate that the exchange rate series follows a Brownian motion (H =1=2) but having a more pronounced bursty (small scale) behaviour. These numerical results con rm that the method in this paper is quite useful in studying the self-similarity=LRD and intermittency in the data. Another method is

V.V. Anh et al. / Journal of Statistical Planning and Inference 80 (1999) 123–135

135

reported in Gao et al. (1998) for the estimation of  and based on the principle of quasi-likelihood. Work on a comparison of the two methods is being undertaken. References Allan, D.W., 1966. Statistics of atomic frequency standards. Proc. IEEE 54, 210–230. Anh, V.V., Lunney, K.E., 1995. Parameter estimation of random elds with long-range dependence. Math. Comput. Modelling 21, 67–77. Anh, V.V., Angulo, J.M., Ruiz-Medina, M.D., 1999. Possible long-range dependence in fractional random elds. J. Statist. Plann. Inference 80, this issue. Anh, V.V., Huang, Y., Kloeden, P.E., 1997. Strong convergence of stochastic Taylor expansions of two-parameter random elds. Stochastic Anal. Appl. 15, 137–149. Beran, J., 1992. Statistical models for data with long-range dependence. Statist. Sci. 7, 404–427. Berman, S.M., 1972. Gaussian sample functions: Uniform dimension and Holder conditions nowhere. Nagoya Math. J. 46, 63–86. Davis, A., Marshak, A., Wiscombe, W., Cahalan, R., 1994. Multifractal characterizations of nonstationarity and intermittency in geophysical elds: observed, retrieved, or simulated. J. Geophys. Res. 99, 8055–8072. Gao, J., Anh, V.V., Heyde, C.C., Tieng, Q., 1998. Parameter estimation of stochastic models for fractal processes (submitted). Geweke, J., Porter-Hudak, S., 1983. The estimation and application of long-memory time series models. J. Time Series Anal. 4, 221–238. Gikman, I.I., 1976. Itˆo formula for two-parameter stochastic integrals. Theory of Random Processes 4, 40–48. Masry, E., 1991. Flicker noise and the estimation of the Allan variance. IEEE Trans. Inform. Theory 37, 1173–1177. Peters, E., 1994. Fractal Market Analysis. Wiley, New York. Prudnikov, A.P., Brychkov, Y.A., Marichev, O.I., 1986. Integrals and Series, vol. 1. Gordon and Breach, New York. Rozanov, Y.A., 1979. Stochastic Markovian elds. In: Krishnaiah, P.R. (Ed.), Developments in Statistics, vol. 2. Academic Press, New York, pp. 203–234. Rozanov, Y.A., 1982. Markov Random Fields. Springer, New York. Ruiz-Medina, M.D., Angulo, J.M., Anh, V.V., 1997. Fractional generalised random elds (submitted). Stein, E.M., 1970. Singular Integrals and Dierentiability Properties of Functions. Princeton University Press, Princeton, NJ. Stroock, D.W., Varadhan, S.R.S., 1979. Multidimensional Diusion Processes. Springer, New York. Yaglom, A.M., 1986. The Correlation Theory of Stationary and Related Random Processes, vol. 1. Springer, New York.