STATISTICAL INFERENCE ON SELF-SIMILAR AND INCREMENT STATIONARY

PROCESSES AND RANDOM FIELDS

By

Jeonghwa Lee

A DISSERTATION

Michigan State University

in partial fulfillment of the requirements

Submitted to

for the degree of

Statistics – Doctor of Philosophy

2019

ABSTRACT

STATISTICAL INFERENCE ON SELF-SIMILAR AND INCREMENT STATIONARY

PROCESSES AND RANDOM FIELDS

By

Jeonghwa Lee

This dissertation is about making statistical inference on self-similar and increment stationary
processes/random fields. Self-similarity and, more generally, fractality are seen in many objects in
nature which have similar features at different scales, for smaller scale and larger scale. The most
well known statistical model that has self-similarity is fractional Brownian motion (fBm). It has
been useful for its self-similarity, increment stationarity, and Gaussianity, and it naturally arises
as the scaling limits of random walk, having many applications in hydrology, telecommunication
network, finance, etc. Some extensions of fBm have been introduced including operator fractional
Brownian motion (OFBM)[19][34] and operator scaling Gaussian random field (OSGRF) [8].
OFBM and OSGRF are multivariate processes, random field, respectively, both have operator
self-similarity/operator scaling property, Gaussianity, and increment stationarity.

The first topic is about estimating Hurst parameter which is a measure for self-similarity in
statistical model. Hurst estimation is examined/developed in OFBM and OSGRF using wavelet
transform and discrete variation method, respectively. The asymptotics of the estimators are derived
in continuous sample path, discrete sample, and discrete noisy sample in OFBM. In OSGRF, the
asymptotics of the estimators are derived in different sampling methods, fixed domain/increasing
domain, and/or samples on the exact directions/ samples on the grid lines. The performance of the
estimators is examined through simulating OFBM and OSGRF, respectively, and the good choice
for scale parameter in wavelet function and discrete variation method is recommended.

The second topic is on measuring dependency between two random fields that are increment
stationary. The dependency is measured in spectral domain by defining and estimating coherence
in multivariate random field. The concept of coherence is originated from multivariate stationary
time series, and it measures correlation between two time series in spectral domain. Recently, the

definition is extended to multivariate random fields by [30]. In this dissertation, the concept of
coherence is extended to multivariate random field with stationary increments, its properties are
examined and its estimation method is developed. Especially, the concept and the estimator method
are applied to OFBM.

Copyright by
JEONGHWA LEE
2019

ACKNOWLEDGEMENTS

This dissertation research is partially supported by NSF grant DMS-1612885.
I would like to thank my advisor, Dr. Xiao, for his invaluable care about his students, his advice
and guidance in my research and future career. I would like to thank my friend Danise for being
good company, spending fun and good time with me here in Michigan. I would like to thank all the
members in Grad IV for our special time together in prayer and fun. I would also like to thank Sarah,
Jessie, Juna and Lori for our good time and our laughs. Finally, I would like to express my deep
appreciation to my family, especially my grandpa, for their love that is beyond my understanding.

v

TABLE OF CONTENTS

.
.

.
.

.
.

.
.

.
.

.
.

.
.

1.1

CHAPTER 1

LIST OF TABLES .
LIST OF FIGURES .

.
.
INTRODUCTION .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . .
Estimation methods of Hurst parameter
1.1.1 Wavelet method .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Discrete variations method . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Coherence of multivariate stationary processes and multivariate random fields . .
.
1.3 Operator fractional Brownian motion and operator scaling Gaussian random fields .
1.3.1 Operator fractional Brownian motion . . . . . . . . . . . . . . . . . . . . .
1.3.2 Operator scaling Gaussian random fields

1
3
3
4
6
8
8
. . . . . . . . . . . . . . . . . . 10

ix

CHAPTER 2 ESTIMATING HURST INDICES IN OPERATOR FRACTIONAL BROW-

NIAN MOTION .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1 Basic Properties of the Wavelet Coefficients of OFBM . . . . . . . . . . . . . . . 12
2.2 Wavelet estimation of Hurst index in OFBM . . . . . . . . . . . . . . . . . . . . . 14
2.2.1 Observing sample path of OFBM . . . . . . . . . . . . . . . . . . . . . . 14
2.2.2 Observing discrete sample paths from OFBM . . . . . . . . . . . . . . . . 20
2.2.3 Hurst estimator in discrete noisy data from OFBM . . . . . . . . . . . . . 24
. . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.4 Estimator for eigenvectors
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.3 Simulation Results

.

.

.

.

CHAPTER 3 COHERENCE OF MULTIVARIATE RANDOM FIELDS WITH STA-

.

.

.

.

.

.

.

.

.

.

.

.

.

. .

TIONARY INCREMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1 Preliminaries
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Coherence: definition and basic properties . . . . . . . . . . . . . . . . . . . . . . 42
3.3 Examples
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3.1 Linear model of coregionalization . . . . . . . . . . . . . . . . . . . . . . 44
3.3.2 Kernel transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3.3 Estimation of Coherence . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3.4 Operator fractional Brownian motion . . . . . . . . . . . . . . . . . . . . . 55
3.4 Simulation and estimation of spectra . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.4.1
. . . . . . . . . . . . . . . . . . . . . . 63
3.4.2 Correlated sample paths with diagonal H: Case II . . . . . . . . . . . . . . 65
3.4.3 Correlated sample paths with diagonizable H: Case III
. . . . . . . . . . . 66

Independent sample paths: Case I

CHAPTER 4 ESTIMATING HURST INDICES IN OPERATOR SCALING GAUS-

4.1 Estimation method .

SIAN RANDOM FIELD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.1.1 Estimation of eigenvectors of E . . . . . . . . . . . . . . . . . . . . . . . 69
Estimation of ratios of H(= 1) to eigenvalues of E . . . . . . . . . . . . . 76
4.1.2

.

.

vi

.
.

.
.

.
.

4.2.2

4.1.3

4.2.1

For H < 1 .
.

4.2 Simulation results .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
. 84
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Independent sample paths on exact directions . . . . . . . . . . . . . . . . 84
4.2.1.1
For H = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.2.1.2
For H < 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
sample surface on a grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
For H = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.2.2.1
For H < 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.2.2.2
CHAPTER 5 CONCLUSION AND DISCUSSION . . . . . . . . . . . . . . . . . . . . . 105
APPENDIX .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
BIBLIOGRAPHY .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

vii

LIST OF TABLES

. . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Table 2.1: Bias and standard error of ˜hi, j
Table 2.2: Bias and standard error of ˜hi, j
Table 2.3: Bias and standard error of ˜hi, j
Table 2.4: Bias and standard error of ˜hi, j
Table 2.5: Mean and standard error of ˆhi, j . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Table 2.6: Bias and standard error of ˜θ j

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

. . . . . . . . . . . . . . . . . . . . . . . . . . . 37

. . . . . . . . . . . . . . . . . . . . . . . . . . . 38

. . . . . . . . . . . . . . . . . . . . . . . . . . . 39

viii

LIST OF FIGURES

Figure 3.1: sample paths X1 and X2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Figure 3.2: The estimates of squared coherence of X1 and X2 . . . . . . . . . . . . . . . . 64

Figure 3.3: sample paths of X1 and X2 . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 65

Figure 3.4: The estimates of squared coherence of X1 and X2 . . . . . . . . . . . . . . . . 65

Figure 3.5: sample paths of X1 and X2 . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 66
Figure 3.6: θ = .2 ∗ π, h1 = .4, h2 = .8, coherence of X1,X2 . . . . . . . . . . . . . . . . . 67
Figure 3.7: θ = .2 ∗ π, h1 = .65, h2 = .8, . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Figure 3.8: θ = .2 ∗ π, h1 = .65, h2 = .8, . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Figure 4.1:

(a): a1 = .1, a2 = .3, (b): a1 = .2, a2 = .3, (c):a1 = .3, a2 = .8, (d):a1 =
.6, a2 = .8, (e): a1 = .8, a2 = .9 . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Figure 4.2: Estimates of cos θ/sin θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

Figure 4.3: Estimates of θ .

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

Figure 4.4: Estimates of a1 and a2 in a fixed domain . . . . . . . . . . . . . . . . . . . . . 92

Figure 4.5: Estimates of a1 and a2 in an increasing domain . . . . . . . . . . . . . . . . . 92

Figure 4.6:

(a): h1 = .6, h2 = .7, H = 75, (b): h1 = .6, h2 = .7, H = .85, (c):h1 =
.1, h2 = .3, H = .4, (d):h1 = .1, h2 = .3, H = .8 (e): h1 = .2, h2 = .3, H = .4 . . . 93

Figure 4.7: Estimates of θ .

.

.

Figure 4.8: Estimates of h1, h2 .

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Figure 4.9: Estimates of θ in a sample surface . . . . . . . . . . . . . . . . . . . . . . . .

. 100

Figure 4.10: Estimates of a1 and a2 in a sample surface with m = 1, 2, 3, 4 . . . . . . . . . . 100

Figure 4.11: Estimates of a1 and a2 in a sample surface with m = 5, 6, 7, 8 . . . . . . . . . . 101

Figure 4.12: Estimates of θ .

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

ix

Figure 4.13: Estimates of h1 and h2 in an fixed domain with m = 1, 2, 3, 4 . . . . . . . . . . . 102

Figure 4.14: Estimates of h1 and h2 in an fixed domain with m = 3, 4, 5, 6 . . . . . . . . . . . 102

Figure 4.15: Estimates of θ .

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Figure 4.16: Estimates of h1 and h2 in an fixed domain with m = 1, 2, 3, 4 . . . . . . . . . . . 103

Figure 4.17: Estimates of h1 and h2 in an fixed domain with m = 3, 4, 5, 6 . . . . . . . . . . . 104

x

CHAPTER 1

INTRODUCTION

Hurst exponent is a measure for self-similarity in time series/stochastic process. It is also called
fractal index or long-memory parameter as it measures roughness of path and/or correlatedness
among the variables in the path. The most well known statistical model that has Hurst exponent
which encompasses all three properties above is fractional Brownian motion. The fractional
Brownian motion (fBm) developed by Mandelbrot and Van Ness (1968)[36] has been studied by
many authors for its importance in modeling processes that have long memory and/or a statistical
self-similarity property.

Fractional Brownian motion (fBm), BH, is a Gaussian process with stationary increments and

covariance function

E(BH(t), BH(s)) =

(|t|2H + |s|2H − |t − s|2H),

1
2

where H ∈ (0, 1) is the Hurst index. Firstly, it is a self-similar process,
= {cH BH(t)}t∈R every c > 0,

{BH(ct)}t∈R

f .d

where f .d.
= denotes equality of all finite-dimensional marginal distributions, which can be easily
derived from the definition of fBm. Secondly, the larger the Hurst index, H ∈ (0, 1), is, the smoother
the sample path of BH is, even though for every H ∈ (0, 1,), the sample function BH(t) is nowhere
differentiable. Thirdly, if H ∈ (.5, 1), the increments of BH, {BH(i + 1) − BH(i)}i∈N, is a stationary
process with long-memory, i.e.

∞


i=1

∞


i=1

cov(BH(i + 1) − BH(i), BH(2) − BH(1)) = ∞,

which means that {BH(i+1)−BH(i)}i∈N has strong dependency and its correlation does not decrease
as fast as a conventional stationary process which has the absolutely summable autocovariace
function like when H ∈ (0, .5)

cov(BH(i + 1) − BH(i), BH(2) − BH(1)) < ∞.

1

Because of the usefulness of the Hurst parameter in capturing the above important properties, there
have been many applications of fBm and an extensive literature on methodology of estimating
Hurst parameter. For example, fBm has been used in image generation and interpolation, texture
classification and the modeling of burst errors in communication channels, 1/ f noise in oscillators
and current noise in metal films and semiconductor devices. Some extensions of fBm have been
introduced including operator fractional Brownian motions [19][34] and operator scaling Gaussian
random fields [8].

A large body of literature studies estimation of the Hurst index, H, or memory parameter, d,
for fractional Brownian motion and more general processes. One way is to use discrete variation
of sample paths [2][10][13][14], and the other method uses wavelet transform [4][22][26][42][43]
[48][50]. For stationary processes, Fourier transform can be useful, or if it is non-stationary then
one uses this method after differencing the process, ∆k X, and making it stationary [22][25][46][47].
Overview and comparison of the two methods, wavelet and Fourier methods, can be found in [22].
In this dissertation, Hurst estimation is developed for an operator fractional Brownian motion
(OFBM) and operator scaling Gaussian random field (OSGRF) using the wavelet transform and
discrete variation method, respectively. In this case, the self-similarity index is a matrix H whose
eigenvalues determine many of their statistical properties. An estimator of the Hurst parameters in
OFBM using wavelet method is proposed in [4] in an increasing domain. In [4], eigenvalue method
was used as the Hurst matrix H was assumed to be diagonizable having eigenvalues h1, h2, and
only hurst index h1 was taken into account for the convergence rate of the estimators. And it was
mentioned that increasing the value of h2 − h1 can increase estimator performance. (see Remarks
4.1,4.4 in [4].)

In Chapter 2, the estimator of the Hurst indices of OFBM is studied in a fixed domain using
the similar method developed in [4]. The asymptotics of the estimator is derived in a continuous
sample path, discrete sample, and discrete noisy sample, respectively. It is revealed that not only
h1, h2, but the dependence structure within the components, i.e. the covariance matrix Γ(1, 1) of
BH(1), and the choice of the scale parameter of wavelet function also affect the performance of the

2

estimator.

In [32], parametric estimation method was developed for estimating Hurst parameters in OS-
GRFs using a fast Fourier transform approximation to harmonizable representation of OSSRFs. In
Chapter 4, a different method is proposed for estimating the Hurst parameters in OSGRF, and its
performance is investigated in different sampling methods.

The wavelet and discrete variation method in fBm is introduced in Section 1.1, and the models
of OFBM and OSGRF will be explained in Section 1.3. The methodology and the results will be
shown in Chapter 2 and Chapter 4 for OFBM and OSGRF, respectively.

The concept of coherence originates from multivariate stationary time series, and it measures
correlation between two time series in the spectral domain. Recently, the definition extended
to multivariate stationary random fields by [30]. The definition of coherence in stationary time
series and random fields are introduced in Section 1.2. In Chapter 3, the concept of coherence is
extended to multivariate random fields with stationary increments, its properties are examined and
its estimation method is developed. Particularly, the concept and the method are applied to OFBM.

1.1 Estimation methods of Hurst parameter

There are several methods for estimating Hurst parameter. Among them are wavelet method

and discrete variation method which will be used in Chapter 2 and Chapter 4, respectively.

1.1.1 Wavelet method

The wavelet transform of BH is defined as

dj,k :=

Þ

R

ψj,k(t)BH(t)dt,

where ψj,k(t) = 2j/2ψ(2jt − k), and where ψ is a wavelet function which has the following
properties.
(W-1) ψ has compact support and

Þ

R

ψ2(t)dt = 1

3

(W-2) ψ has M vanishing moments:

tmψ(t)dt = 0 for all m = 0, ..., M − 1.

R

Þ
j,k) = 2−j(2H+1)Þ

E(d2

Then for any fixed scale j, the wavelet coefficients {dj,k}k is mean-zero stationary Gaussian process
with variance

Note thatÞ

R ψ(t)ψ(s)|t − s|2Hdtds = E(d2

0,k) and the covariance between d0,k and d0,k(cid:48) decreases

ψ(t)ψ(s)|t − s|2Hdtds.

R

as fast as cov(d0,k, d0,k(cid:48)) ∼ |k − k(cid:48)|2H−2M as k − k(cid:48) → ∞.
Since

log d2

j ∼ −2jH log 2 + log E(d2

0,k),

the estimator of H is the least square estimate on log-regression of d2
the sample mean of {d2

j,k; k = 0, 1, ..., n}.

j on j = J1, ..., J2 where d2

j is

where

J2

J2


j=J1

ˆH = −1
2

J2


wj log2 d2
j

wj = 0,

jwj = 1.

j=J1

j=J1

1.1.2 Discrete variations method
This sections is mainly from [15]. Define a = {aq, q = 0, 1, ..., (cid:96)} as a filter of length (cid:96) + 1 with
order p ≥ 1, that is a vector with (cid:96) + 1 components satisfying

(cid:96)


q=0

Define

(cid:96)


q=1

qpaq (cid:44) 0.

q jaq = 0 for j = 0, ..., p − 1 and

(cid:96)


q=0

H(i) :=
Ba

aqBH(i − q).

4

The covariance and correlation functions are given by

E(Ba

H(j)Ba

H(i + j)) = −1
2

aqar|q − r + i|2H

(cid:96)


q,r=0

for i, j ∈ Z. By filtering BH with a, the correlation of Ba
|i|2H−2p as i → ∞. For m ≥ 1, define dilated filters am of a filter a by

H decreases as cor(Ba

H(j)Ba

H(j + i)) ∼

ai/m

0

am
i =

m(cid:96)


q,r=0

H )2(cid:17)
(cid:16)(Bam

if i/m is integer

otherwise

(cid:96)


q,r=0

for i = 0, ..., m(cid:96). That is, am is a filter with length m(cid:96) + 1 with order p. Then it is derived that

E(Bam

H (j)Bam

H (i + j)) = −1
2

and

am
q am

r |q − r + i|2H = −1
2

aqar|mq − mr + mi|2H,

= m2H πa

H(0),

E
H )2 is the sample mean of {(Bam
H(j))2) and (Bam

H(0) := E((Ba

where πa
estimator of H, ˆH, is the ordinary least square estimate from the regression of log(Bam
other words, since

H (j))2, j = 1, ..., n}. Then the
H )2 on m. In

the estimator ˆH is

H(0),

log(Bam

H )2 ∼ 2H log m + log πa
(cid:16)

H )2(cid:17)

log(Bam

A(cid:48)
2||A||2

m=M1,...,M2 ,

ˆH =


M2

m=M1

log m and M1 ≤ m ≤ M2. In particular, if m = 2M1, ..., 2M2

where Am = log m−
then ˆH can be expressed as

M2−M1+1

1

where

ˆH =

M2


M2


r=M1

1
2

wr log2 (Ba2r
H )2

,

M2


wr = 0,

rwr = 1.

r=M1

r=M1

5

1.2 Coherence of multivariate stationary processes and multivariate random

fields

This section is mainly from 11.6 in [11]. Let Xt = (Xt,1, Xt,2)(cid:48) be a bivariate stationary time

series with mean zero and covariances ci j(h) = E(Xt+h,i Xt, j) satisfying

∞

∞


h=−∞

1
2π

h=−∞

|ci j(h)| < ∞ i, j = 1, 2.

e−ihλC(h) =(cid:169)(cid:173)(cid:173)(cid:171) f11(h)

f21(h)

(cid:170)(cid:174)(cid:174)(cid:172)

f12(h)
f22(h)

Then the matrix

f(λ) =

is called the spectral density matrix or spectrum of Xt where C(h) is the 2 × 2 matrix with (i, j)
element ci j(h), i, j = 1, 2. Especially, fii(·) is called the spectral density of the univariate series
{Xt,i}, i = 1, 2, and f12(·) is called the cross spectrum or cross spectral density of {Xt,1} and {Xt,2}.
Since cii(·) is symmetric around zero,
fii(·) is real-valued and symmetric around zero, but not
fi j(·), i (cid:44) j. It is followed from the above definition that

Þ π

−π

C(h) =

eihλ f(λ)dλ

Þ
(−π,π] eitλdZi(λ)

and Xt has spectral representation,

Xt,i =

where

and

fi j(λ)dλ = E(dZi(λ)dZ j(λ)),

i, j = 1, 2,

E(dZi(λ)dZ j(µ)) = 0 for λ (cid:44) µ and i, j = 1, 2.

fi j(·) = fji(·) and therefore f(λ) = f ∗(λ) where * denotes complex conjugate
It follows that
transpose. Also, for any a = (a1, a2)(cid:48) ∈ C2, a∗ f(λ)a is the spectral density of {a∗Xt}. Therefore,
a∗ f(λ)a ≥ 0 and the matrix f(λ) is non-negative definite.

6

The correlation between dZ1(λ) and dZ2(λ) is called coherence, γ12(λ), at frequency λ,

γ12(λ) = f12(λ)/[ f11(λ) f22(λ)]1/2.

Note that 0 ≤ |γ12(λ)|2 ≤ 1 for −π ≤ λ ≤ π, and γ12(λ) close to one indicates a strong linear
relationship between dZ1(λ) and dZ2(λ).

The above notion of coherence in multivariate stationary processes was extended to multivariate
random fields in [30]. Suppose X(t) = (X1(t), ..., Xp(t))(cid:48) ∈ Cp is a mean-zero p-variate weakly
statinary random field on t ∈ Rd with a matrix-valued covariance function C(h) = (Ci j(h))p
i, j=1
where Ci j(h) = Cov(Xi(t + h)Xj(t)). For complex-valued stationary processes, Cov(Xi(t)Xj(s)) =
E(Xi(t)Xj(s)) so that Ci j(h) = Cji(−h).
Theorem 1.2.1 (Cramér (1940)). A matrix valued function C : Rd → Cp×p, C = (Ci j)p
non-negative definite if and only if

is a

i, j=1

Ci j(h) =

Rd

eiω(cid:48)h fi j(ω)dω

Þ

Þ

for i, j = 1, ..., p such that the matrix f(ω) = ( fi j(ω))p

is nonnegative definite for all ω ∈ Rd .
Like before, fi j(ω) are called spectral and cross-spectral densities for the marginal and cross-

i, j=1

covariance functions Ci j(h), and fi j(ω) = fji(ω). Note that with covariance function Ci j(·)

fi j(ω) =

1
(2π)d

Rd

e−iω(cid:48)hCi j(h)dh

and the coherence function is defined

γ(ω) =

(cid:112) f11(ω) f22(ω)

f12(ω)

for ω ∈ Rd, and fii(ω) > 0. If fii(ω) = 0, then γ(ω) = 0. The squared coherence satisfies
0 ≤ |γ(ω)|2 ≤ 1 for all ω by Theorem 1.2.1 and |γ(ω)|2 close to one indicates X1 and X2 has strong
linear relationship at particular frequency bands.

7

1.3 Operator fractional Brownian motion and operator scaling Gaussian

random fields

Operator fractional Brownian motions and operator scaling Gaussian random fields are exten-
sions of fBm with operator self-similar/scaling property in multivariate stochastic processes and
random field.

1.3.1 Operator fractional Brownian motion

Operator fractional Brownian motions (OFBM) arise naturally as the scaling limit of various
multivariate random walks with (long-range) dependent steps on lattices normalized by linear
operators, from asymptotic analysis of nonstationary fractionally integrated multivariate random
sequences, from central limit theorem in Hilbert spaces, or from modeling heavy traffic in queuing
systems. See Pitt [44], Marinucci and Robinson [37], Rackauskas and Suquet [45], Düker [20, 21],
Delgado, R [17], Majewski [35], just to mention a few.

Let X = {X(t), t ∈ Rd} be a p-variate random field. We say that X is an OFBM with exponent
H if it is a Gaussian field with mean zero and stationary increments, and satisfies the following
operator scaling (operator self-similar) property: For any c > 0,

{X(ct), t ∈ Rd} f .d.

= {cHX(t)},

where H is a linear operator on Rd and cH = exp(H ln c) =
∞

k=1 logk(c)Hk/k! for c > 0. For
convenience, we will not distinguish the operator H from its associated matrix relative to the
standard basis of Rd. One can refer to the monographs [38, 41] for more on operator stable
distributions.

(1.1)

Denote by λH and ΛH the minimum and the maximum of the real parts of the eigenvalues of
H, respectively. Mason and Jurek [38] showed that 0 ≤ λH ≤ ΛH ≤ 1 and that X is non-trivial
if and only if λH > 0. It turns out [39] that many sample path properties of OFBM are reflected
through the real parts of the eigenvalues of the operator H.

Recently, there has been much interest in studying OFBM in the literature.

8

Maejima and Mason [34] established a time domain construction of OFBM , while Mason and
Xiao [39] gave a spectral domain construction à la Itô-Yaglom, for general d ≥ 1. When d = 1,
Didier and Pipiras [19] found both time and spectral domain representations for OFBM which is
more general than those given in [34, 39].

An Rp-valued stochastic process {XH(t)}t∈R is said to be OFBM if it is an increment stationary,
zero-mean, and operator self-similar Gaussian process, i,e. when its law scales according to a matrix
exponent H, i.e. for c > 0, {XH(ct)}t∈R
= {cHXH(t)}t∈R. When the operator H is diagonal, it is
called the multivariate fractional Brownian motion.

f .d

In the spectral domain, with the mild assumption that eigenvalues of H satisfy

0 < Re(hk) < 1,

k = 1, 2, ..., p.

(1.2)

Þ

(cid:2)((t − u)H−(1/2)I

)M−(cid:3) B(du),

OFBM admits the following representation:

Þ

eitx − 1

XH(t) =

(x−(H−(1/2)I)

C + x−(H−(1/2)I)

C) ˜B(dx)

+

R

−

ix

(1.3)
for some p × p complex-valued matrix C , where x± = max{±x, 0} and ˜B(dx) is a multivariate
complex-valued Gaussian measure such that ˜B(−dx) = ˜B(dx) and E ˜B(dx) ˜B(dx)∗ = dx. In the
time domain, in addition to (2), if Re(hk) (cid:44) 1/2, k = 1, 2, ..., p, then OFBM has the following
representation:

− (−u)H−(1/2)I

+

)M+ + ((t − u)H−(1/2)I

− (−u)H−(1/2)I

+

R

(1.4)
where M−, M+ are p × p matrices with real-valued entries and B(du) is a multivariate real-valued
Gaussian measure. The representation (4) is obtained from (3) by taking Fourier transformation of
deterministic kernel in (3).

−

−

Unlike the univariate case, OFBM does not have the following property in general.

E(XH(t)XH(s)∗) = E(XH(s)XH(t)∗),

(1.5)

where * represents Hermitian transposition. In OFBM, time reversibility of the process is equivalent
to (1.5). Didier and Pipiras (Theorem 5.1, Corollary 5.1 in [19]) provided some conditions on C

9

and M in (1.3) and (1.4) respectively for OFBM to be time-reversible. With this condition, one
has the following equation analogous to fBm.

E(XH(t)XH(s)∗) =

1
2

(|t|H Γ(1, 1)|t|H∗

+ |s|H Γ(1, 1)|s|H∗

+ |t − s|H Γ(1, 1)|t − s|H∗),

(1.6)

where Γ(1, 1) = EXH(1)XH(1)∗.

In Chapter 2, OFBM with d = 1, p = 2 is used with scaling matrix H, and its eigenvalues (Hurst

indices) h1, h2. In Chapter 3, OFBM with d ≥ 1, p = 2 is used.

1.3.2 Operator scaling Gaussian random fields
A scalar valued random field {X(t)}t∈Rd is called operator-scaling if for some d × d matrix E with
positive real parts of the eigenvalues and some H > 0 we have

{X(cEt)}t∈Rd

f .d.

= {cH X(t)}t∈Rd

for all c > 0.

(1.7)

If E is the identity matrix, then X has self-similar property with Hurst index H, {X(ct)}Rd
{cH X(t)}Rd, and fractional Brownian field is well known class of such field. If E is diagonizable
= {cH/λi Xi(t)}t∈R for any c > 0, where
matrix, Hui = λiui, i = 1, ..., d, then {Xi(ct)}t∈R
Xi(t) = X(tui). As this model has different Hurst indices H/λi along the different coordinates ui, it
is expected to be useful to capture such features observed in nature object. Among the applications
of this model is ground water hydrology, [7].

f .d.
=

f .d.

In [8], harmonizable representation and moving average representation of operator scaling

stable random fields (OSSRFs) are defined as below: Harmonizable representation is

Xψ(t) = Re

(ei(cid:104)t,ξ(cid:105) − 1)ψ(ξ)−H−q/αWα(dξ), t ∈ Rd,

where ψ is a continuous, Et-homogeneous function such that ψ(ξ) (cid:44) 0 for ξ (cid:44) 0, 0 < α ≤ 2, and
H ∈ (0, a1), a1 is the smallest real parts of the eigenvalue of E. Moving average representation is

Þ

Rd

Þ

Xϕ(t) =

Rd

(ϕ(t − y)H−q/α − ϕ(−y)H−q/α)Zα(dy), t ∈ Rd,

10

where ϕ be an E-homogeneous, (β, E)−admissible function, 0 < β, 0 < α ≤ 2, and 0 < H < β. ϕ
is E-homogeneous if ϕ(cEt) = cϕ(t), for all c > 0 and t ∈ Rd \{0}, and is called (β, E)−admissible
if for any 0 < A < B, there exists positive constant C such that, for A ≤ ||y|| ≤ B, τ(t) ≤ 1 ⇒
|ϕ(t + y) − ϕ(y)| ≤ Cτ(t)β. It turns out that the harmonizable representation is more flexible in
the class of possible functions ϕ whereas the moving average representation is more restrictive.
However they both satisfy (1), have stationary increments and are continuous in probability.

In [9], explicit covariance functions of operator scaling Gaussian random fields are provided

making a fast and exact method of simulation available in these classes. They define a function

τE(t) =(cid:0) d


|(cid:104)t, ui(cid:105)|2ai(cid:1)1/2

i=1

, t ∈ Rd,

(1.8)

where 1/ai = λi > 0, ui, i = 1, 2, ..., d, are eigenvalues and eigenvectors of a diagonizable matrix
E, and find out that for H ∈ (0, 1], υE,H(t) = τE(t)2H is the semi-variogram funtion of a centered
Gaussian random field that has stationary increments and satisfies (1.7). Recall semi-variogram of
X is defined by υE,H(h) = 1

2 E(X(t + h) − X(t))2.

Throughout this dissertation, var, cov means variance and covariance respectively, vec means

vectorization of a matrix, especially for symmetric matrix vec((cid:169)(cid:173)(cid:173)(cid:171)a b

b d

(cid:170)(cid:174)(cid:174)(cid:172)) = (a, b, d)(cid:48), →d indicates

convergence in distribution, C and c are a generic constant matrix and scalar respectively, which
are not related to j, n, H. For two matrices A = (ai j), B = (bi j), A ≤ B means ai j ≤ bi j for all i, j.
a ∼ b means a/b → 1, an = Op(bn) means ab/bn is bounded in probability which also implies an
converges in probability if bn → 0.

11

ESTIMATING HURST INDICES IN OPERATOR FRACTIONAL BROWNIAN MOTION

CHAPTER 2

2.1 Basic Properties of the Wavelet Coefficients of OFBM

Let XH(t) be an operator fractional Brownian motion satisfying (1.1-1.3),(1.5-1.6). For nota-
tional convenience, XH(t) is denoted as Xt = (Xt,1, Xt,2) hereafter. Suppose we observe Xt in a
fixed domain t ∈ [0, 1]. The assumptions below are given to wavelet function and OFBM.

ASSUMPTION (W1) :
ψ is a wavelet function with two vanishing moments and has compact support [0, S].

ASSUMPTION (OFBM1):
Time reversibility holds, i.e., {Xt} f .d.
ASSUMPTION (OFBM 2):
Xt is a bivariate OFBM with scaling matrix H = Pdiag(h1, h2)P−1, 0 < h1 ≤ h2 < 1, i.e. H is

= {X−t}.

diagonizable with P =(cid:169)(cid:173)(cid:173)(cid:171)cos θ − sin θ

cos θ

sin θ

For j ∈ Z and for some large n greater than 2j, define the set Sj = {ki, j : ki, j = 2j
Denote the random wavelet coefficients of Xt by

n i, i = 1, ..., n}.

(cid:170)(cid:174)(cid:174)(cid:172) .
Þ

dj,ki, j

=

R

ψj,ki, j(t)Xtdt

ki, j ∈ Sj,

where ψj,ki, j(t) = 2j/2ψ(2jt − ki, j). The sequence {dj,ki, j} is centered bivariate Gaussian and
stationary with respect to the location parameter i for fixed j. It has the following properties.
(P1)

{dj,ki, j ; i ∈ [1, ..., n]} f .d.

= {2−j(H+1/2)d0,ki, j ; i ∈ [1, ..., n]}.

12

(cid:48)

1)

(P

{dj,ki, j, dj+1,ki, j+1; i ∈ [1, ..., n]} f .d.

= {2−j(H+1/2)d0,ki, j, 2−j(H+1/2)d1,ki, j+1; i ∈ [1, ..., n]}.

(P2)

where c(ψ, H) = 1/2Þ ψ(t)ψ(s)|t − s|H Γ(1, 1)|t − s|Hdsdt.

var(dj,ki, j) = 2−j(H+1/2)c(ψ, H)2−j(H+1/2)

,

) = 2−j(H+1/2)cov(d0,ki, j, d0,ki(cid:48), j

)2−j(H+1/2) , and it decays hyperboli-
(P3) cov(dj,ki, j, dj,ki(cid:48), j
cally fast (each element of the covariance matrix decays hyperbolically fast) as |ki, j − ki(cid:48), j| → ∞
due to two vanishing moments of ψ ([13],[49]).

(P4) |cov(d0,ki, j, d0,ki(cid:48), j

)(cid:96),p| ≤ c(1 + |ki, j − ki(cid:48), j|)2(h2−2),

(cid:96), p = 1, 2, h1 ≤ h2.

) = 2−j(H+1/2)cov(d0,ki, j, d1,ki(cid:48), j+1

)2−j(H+1/2) and

(cid:48)

(P

4) cov(dj,ki, j, dj+1,ki(cid:48), j+1
Þ
cov(d0,ki, j, d1,ki(cid:48), j+1
)
Þ

=

1
2
1
2

=

ψ(t)ψ(s)|2−1t − s + ki, j − 2−1ki(cid:48), j+1|H Γ(1, 1)|2−1t − s + ki, j − 2−1ki(cid:48), j+1|Hdtds
ψ(t)ψ(s)|2−1t − s + ki, j − ki(cid:48), j|H Γ(1, 1)|2−1t − s + ki, j − ki(cid:48), j|Hdtds,

therefore

|cov(d0,ki, j, d1,ki(cid:48), j+1

)(cid:96),p| ≤ c(1 + |ki, j − ki(cid:48), j|)2(h2−2)

,

(cid:96), p = 1, 2, h1 ≤ h2.

(P5) By the stationary increments of Xt and wavelet assumption (W1), {(dj,ki, j, ..., dj+L,ki, j+L); i =
1, ..., n} is a multivariate stationary process for a fixed j, L ∈ Z.

13

where dj,k is a vector of wavelet coefficients of sample path of OFBM,

k∈Sj

Q j = 

Þ
dj,k =(cid:169)(cid:173)(cid:173)(cid:171)
Þ

dj,kd(cid:48)

j,k/n,

R Xs,1ψj,k(s)ds
R Xs,2ψj,k(s)ds

(cid:170)(cid:174)(cid:174)(cid:172) .

2.2 Wavelet estimation of Hurst index in OFBM

2.2.1 Observing sample path of OFBM

Throughout this chapter, assumptions (W1, OFBM1, OFBM2) are assumed. Assume further that a
continuous sample path of Xt is observed for t = [0, 1]. Then, the estimation strategy is as follows.
Define a 2 × 2 sample covariance matrix Q j in the following way:

As n increases, Q j gets closer to its expectation,

Q j ∼ 2−j(H+1/2)c(ψ, H)2−j(H+1/2)

.

It is shown later that the eigenvalues of Q j satisfy ρ(Q j)i ∼ ci2j(2hi+1), i = 1, 2, for some constants
ci, where ρ(Q j)i are the eigenvalues of Q j. The estimator hi, j of hi is the log regression of ρ(Q j)i
on j.

(2.1)

where

j


(cid:96)=jL

1
2

hi, j =

w(cid:96) = 0,

j


w(cid:96) log2 ρ(Q(cid:96))i − 1
2,
j

P−12−j(H+1/2)d0,k =(cid:169)(cid:173)(cid:173)(cid:171)2−j(h1+1/2)

(cid:96)w(cid:96) = 1.

0

(cid:96)=jL

(cid:96)=jL

(cid:170)(cid:174)(cid:174)(cid:172) P−1d0,k .

0

2−j(h2+1/2)

Note that the eigenvalues of Q j are the same as those of P−1Q j P, and

Since our estimators use eigenvalues of the matrix Q j and ρ(Q j) = ρ(P−1Q j P), it is assumed
that the matrices P, P−1 were premultiplied on both sides of the matrix Q j and assume that

14

H = diag(h1, h2). i.e. ρ(Q j) = ρ(P−1Q j P) so we analyze P−1Q j P instead. Therefore, from now
on, we assume H = diag(h1, h2) and

var(dj,k) =(cid:169)(cid:173)(cid:173)(cid:171)2j(h1+1/2)

0

(cid:170)(cid:174)(cid:174)(cid:172) P−1c(ψ, H)P(cid:169)(cid:173)(cid:173)(cid:171)2j(h1+1/2)

0

0

2j(h2+1/2)

(cid:170)(cid:174)(cid:174)(cid:172)

0

2j(h2+1/2)

without loss of generality.
Proposition 2.2.1. Let 2j−N → m, where m is a constant, m ∈ [0, 1] , as n → ∞. Then as n → ∞,

where Vj =(cid:169)(cid:173)(cid:173)(cid:171)2j(h1+1/2)




0
{E[(d0,kd(cid:48)

2j/2{vec(Vj(Q j − EQ j)Vj)} →d N(0, Σj),
0

(cid:170)(cid:174)(cid:174)(cid:172) and the elements of Σj are

2j(h2+1/2)

2j
n

0,k(cid:48))(cid:96),(cid:96)(cid:48)]E[(d0,k(cid:48)d(cid:48)

0,k)p,p(cid:48)] + E[(d0,kd(cid:48)

0,k(cid:48))(cid:96),p(cid:48)]E[(d0,k(cid:48)d(cid:48)

lim
n→∞
k,k(cid:48)∈Sj
(cid:96), p, (cid:96)(cid:48)p(cid:48) = 1, 2.
Proof. Note that {d0,k; k ∈ Sj} is a stationary Gaussian process for fixed j, and let

0,k)p,(cid:96)(cid:48)]},

f(cid:96),p(d0,k) := d0,k,(cid:96)d0,k,p − E[d0,k,(cid:96)d0,k,p],
Then, f(cid:96),p has a generalized Hermite rank 2 (see (2.2) in [5]), and

(cid:96), p = 1, 2.

d0,k,(cid:96)d0,k,p − E[d0,k,(cid:96)d0,k,p]

n

:= 


k∈Sj

f(cid:96),p(d0,k)

n

.

2j(h(cid:96)

+1/2)2j(hp+1/2)(Q j − EQ j)(cid:96),p

Note that for (cid:96), p = 1, 2,

where r(cid:96),p(k) = cov(d0,0,(cid:96), d0,k,p), by (P
Theorem 4 in [5],

n1/2


k∈Sj

(cid:20)

1
n

lim
n→∞

|r(cid:96),p(k)|2 < ∞,

k∈Sj
(cid:48)
4) and the fact that Sj = { i2j

(cid:18)2j/2d0,k
= 2j/2


n1/4

f1,1

(cid:19)

(cid:19)
(cid:18)2j/2d0,k
(cid:18) f1,1(d0,k)

n1/4

, f1,2

, f2,2
f1,2(d0,k)

,

n

,

(cid:18)2j/2d0,k

n1/4
f2,2(d0,k)

(cid:19)(cid:21)
(cid:19)

k∈Sj
converges to multivariate normal distribution.

n

(cid:3)

n , i = 1, 2, ..., n.}. Then, by

f .d.

= 




k∈Sj

2j
n

n

15

Corollary 2.2.2. Let 2j−N → m. (see Proposition 2.2.1 for m) Then as n → ∞,

vec(Vj(Q j − EQ j)Vj)

2−1/2vec(Vj−1(Q j−1 − EQ j−1)Vj−1)

...

2−jL/2vec(Vj−jL(Q j−jL − EQ j−jL)Vj−jL)

→d N(0, Σ∗
j),

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

is a block matrix where the (m, m)-th block is Σj−m+1, and (m, n)-th block, m > n, consists

where Σ∗
j
of
n→∞ 2(m−n)/22(hl +hp+1)(m−n)
lim




k∈Sj−m+1,k(cid:48)∈Sj−n+1

{E[(d0,kd(cid:48)

m−n,k(cid:48))l,l(cid:48)]E[(d0,k(cid:48)d(cid:48)

m−n,k)p,p(cid:48)]

2j
n

+E[(d0,kd(cid:48)

m−n,k(cid:48))l,p(cid:48)]E[(d0,k(cid:48)d(cid:48)

m−n,k)p,l(cid:48)]},

l, p, l(cid:48)

, p(cid:48) = 1, 2.

Proof. Without loss of generality, let jL = j − 1. By (P

(cid:48)

4),

(dj,k,1, dj,k,2, dj−1,k(cid:48),1, dj−1,k(cid:48),2)k∈Sj,k(cid:48)∈Sj−1
is a vector valued stationary Gaussian process for a fixed j. Define f j
proposition. Since both of their generalized rank are 2, and by (P

4)

(cid:48)

l,p and f j−1

l,p

as in the above

lim
n→∞

2j
n

cov(d0,0,(cid:96), d0,k,p)2 < ∞, lim
n→∞

2j
n

cov(d1,0,(cid:96), d0,k,p)2 < ∞,

for (cid:96), p = 1, 2. Therefore, the result is derived from Theorem 9 in [5].

(cid:3)

Proposition 2.2.3. If H is diagonizable and has real eigenvalues, then for i = 1, 2,

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

2j/2




k∈Sj




k∈Sj



O(2−2j(h2−h1)),

0,

if h1 (cid:44) h2
if h1 = h2,

| log ρ(EQ j)i + j(2hi + 1) log 2 + ci| =

where c1 = a, c2 = d − b2
a .

16

Proof. To calculate the eigenvalues of EQ j, we only need to calculate the eigenvalues of the matrix,

(cid:169)(cid:173)(cid:173)(cid:171)2−j(h1+1/2)

0

0

2−j(h2+1/2)

0

b d

(cid:169)(cid:173)(cid:173)(cid:171)2−j(h1+1/2)
(cid:169)(cid:173)(cid:173)(cid:171)a b
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172) = Pc(ψ, H)P−1

(cid:169)(cid:173)(cid:173)(cid:171)a b

b d

(cid:170)(cid:174)(cid:174)(cid:172) ,

0

2−j(h2+1/2)

where

by Assumption(OFBM2) and (P2).
Let a(j) = 2−j(2h1+1)a, d(j) = 2−j(2h2+1)d, and b(j) = 2−j(h1+1/2)2−j(h2+1/2)b. Then the eigen-
values ρ(EQ j)i, i = 1, 2 of EQ j are

(a(j) + d(j)) ±(cid:113)(a(j) + d(j))2 − 4(a(j)d(j) − b2(j))
ρ(EQ j)i = 2−j(2h+1)(a + d) ±(cid:112)(a + d)2 − 4(ad − b2)

2

.

2

If h1 = h2 = h, the result follows since

.

(2.2)

where θ j ∈ [0, 1]. Therefore

Now assume h1 (cid:44) h2. By Taylor’s expansion,

(cid:113)(a(j) − d(j))2 + 4b2(j) =
ρ(EQ j)1 = 2−j(2h1+1)(cid:16)
ρ(EQ j)2 = 2−j(2h2+1)(cid:16)

d −

a +

,

4b2(j)

(cid:113)(a(j) − d(j))2 +
2(cid:113)(a(j) − d(j))2 + θ j4b2(j)
(cid:113)(a2−2j(h1−h2) − d)2 + θ j4b22−2j(h1−h2)
(cid:113)(d2−2j(h2−h1) − a)2 + θ j4b22−2j(h2−h1)

4b2

4b2

4

4

(cid:17)
(cid:17)

.

.

Similarly,

The result is derived by using the Mean value theorem and the fact that 22j(h2−h1) → ∞.

17

(2.3)

(2.4)

(cid:3)

Let

g1(a, b, d) := 2−j(2h1+1)(cid:16)
g2(a, b, d) := 2−j(2h2+1)(cid:16)

a +

d −

4b2

(cid:113)(a2−2j(h1−h2) − d)2 + θ j4b22−2j(h1−h2)
(cid:113)(d2−2j(h2−h1) − a)2 + θ j4b22−2j(h2−h1)

4b2

4

4

(cid:17)
(cid:17)

,

(2.5)

.

(2.6)

Proposition 2.2.4. If H is diagonizable and has real eigenvalues, then
2j/2(log ρ(Q j)i − log ρ(EQ j)i) →d N(0, Σρ)
Σg(1), and where Σ is from Proposition 2.2.1, and g(1) is the vector of

for i = 1, 2, where Σρ = g(1)(cid:48)
the first derivatives of (g1, g2).
Proof. Let

Vj EQ jVj =(cid:169)(cid:173)(cid:173)(cid:171)a b
(cid:170)(cid:174)(cid:174)(cid:172) ,
(cid:169)(cid:173)(cid:173)(cid:171)a b
(cid:170)(cid:174)(cid:174)(cid:172) = Pc(ψ, H)P−1

b d

b d

Note that

(cid:170)(cid:174)(cid:174)(cid:172) .

VjQ jVj =(cid:169)(cid:173)(cid:173)(cid:171)aj bj
2j−1

(cid:170)(cid:174)(cid:174)(cid:172) law=

d0,kd(cid:48)

bj dj

k=0

0,k/2j−1.

(cid:169)(cid:173)(cid:173)(cid:171)aj bj

bj dj

and

Note that aj, bj, dj are random variables which converge to a, b, d, respectively, with the convergence
rate of 2(j−1)/2 by Proposition 2.2.1. From the previous Proposition 2.2.3 and delta method, the
result is derived, since log ρ(Q j)i − log ρ(EQ j)i = gi(aj, bj, dj) − gi(a, b, d), which has derivative
for each element. By Proposition 2.2.1, the result is derived with Σρ = g(1)(cid:48)
Σg(1), where Σ is
from Proposition 2.2.1, and g(1) is the vector of the first derivatives of (g1, g2) with three variables
(cid:3)
a, b, d.

Corollary 2.2.5. For i = 1, 2,

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

2j/2

log ρ(Q j)i − log ρ(EQ j)i
log ρ(Q j−1)i − log ρ(EQ j−1)i

..

log ρ(Q j−jL)i − log ρ(EQ j−jL)i

18

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

→d N(0, Σρ∗),

where Σρ∗ is a (jL + 1) × (jL + 1) matrix whose ((cid:96), m)-th element is Σρ∗
Σ∗
((cid:96),m) is the ((cid:96), m)-th block matrix in Corollary 2.2.2.
Proof. The result follows from Propositon 2.2.4 and Corollary 2.2.2.

(cid:96),m

= g(1)(cid:48)

Σ∗
((cid:96),m)g(1), and

Define

(cid:3)

(2.7)

(2.8)

(cid:3)

hi, j = −1
2

w(cid:96) log2 ρ(Q(cid:96))i − 1
2,

j

j


(cid:96)=jL

i, j = −1
hE
2
Proposition 2.2.6. As j → ∞, for i = 1, 2,

(cid:96)=jL

w(cid:96) log2 ρ(EQ(cid:96))i − 1
2 .

2j/2(hi, j − hE

i, j) →d N(0, σh),
where σh = W(cid:48)Σρ∗W, and W = (wj/(2 log 2), ..., wjL/(2 log 2))(cid:48)
Proof. The result is derived from Corollary 2.2.5.

Proposition 2.2.7. For i = 1, 2,

i, j − hi| = Op(ci,1
|hE
ci

where f1(j) = 2−j(2h2−2h1). If h2 − h1 → 0,

( f1(j) − f1(jL))),

|hE
i, j − hi| =

ci,1
ci

(h2 − h1)2−2jL(h2−h1) + o((h2 − h1)2−2jL(h2−h1))

for c1,1 = 2b2/a2, c2,1 = 2b2d/a3.
Proof. From (2.3),(2.4) in Proposition 2.2.3,

ρ(EQ j)i = 2−j(2hi+1)(ci + ci,12−2j(h1−h2) + o(2−2j(h1−h2))).

Therefore

log ρ(EQ j)i = −j(2hi + 1) log 2 + log ci +
( f1(j) − f1(jL)).

i, j = hi − ci,1
hE
ci

ci,1
ci

2−2j(h2−h1) + o(2−2j(h2−h1))

(2.9)

(2.10)

19

If h2 − h1 is small,

2−2j(h2−h1) = 2−2jL(h2−h1) + 2−2jL(h2−h1)(−2j(h2 − h1) + 2jL(h2 − h1)) log 2

+ o(2−2jL(h2−h1)(−2j(h2 − h1) + 2jL(h2 − h1)))

and

jU


i, j − hi = −1
hE
2
= −ci,1
ci

w(cid:96) log2 ρ(EQ(cid:96))i − .5 − hi

(cid:96)=jL
(h2 − h1)2−2jL(h2−h1) + o(−(h2 − h1)2−2jL(h2−h1)).

2.2.2 Observing discrete sample paths from OFBM

(cid:3)

With the discrete sample path of OFBM Xn
i
following:

= (Xn

i,1, Xn

i,2) = (Xi/n,1, Xi/n,2), i = 1, ..., n, define the

(cid:16)Þ k/2j +(q+1)/2N

k/2j +q/2N

=

S2N−j−1

j = 


q=0
˜d(cid:48)n
˜dn
j,k
j,k
n

k∈Sj

˜dn
j,k

˜Qn

(cid:17)Xn

k2N−j +q

ψj,k(t)dt

(2.11)

(2.12)

for j ≤ log2 n and 2N = n. Note that for fixed j, { ˜dn
mean zero and the covariance between ˜dn
similarly to (P4). For these reasons, the analogous to Proposition 2.2.1 holds for ˜Q j .
Proposition 2.2.8. Let 2j−N → m. (m in Proposition 2.2.1.) Then as n → ∞,

j,k}k is a stationary Gaussian process with
j,k(cid:48) decays hyperbolically fast as |k − k(cid:48)| → ∞

j,k and ˜dn

2j/2{vec(Vj( ˜Qn

j − E ˜Qn

j)Vj)} →d N(0, ˜Σj),

where the elements of ˜Σj are defined in (2.15).
If j <<

log2 n, then

2h1
.5+2h1

2j/2{vec(Vj( ˜Qn

j − EQ j)Vj)} →d N(0, ˜Σj).

20

Proof. Note that

˜dn
j,k

= 
S2N−j−1

q=0
Let dn∗
ψ0,0(t)dt
0,k
decaying as |k − k(cid:48)| → ∞, and the limit of dn∗

q=0

ψj,k(t)dt

(cid:17)Xn
(cid:16)Þ (q+1)2j/2N
(cid:17)X2N−j

q2j/2N

=

f .d.

q=0

k/2j +q/2N

(cid:16)Þ k/2j +(q+1)/2N
S2N−j−1

= 2−j(H+1/2) S2N−j−1

(cid:16)Þ (q+1)2j/2N
d0,k,

Sm−1

q2j/2N

dn∗
0,k

q=0

qm

=

k2N−j +q

0,kd(cid:48)n∗
. Then E(dn∗
0,k exists in almost sure sense.

(cid:16)Þ (q+1)m

(cid:17)Xm

km+q

ψ(t)dt

if m = 0
if m ∈ (0, 1].

k2N−j +q

ψ0,0(t)dt

(cid:17)X2N−j

k2N−j +q

(2.13)

.

(2.14)

0,k(cid:48)) is hyperbolically

d∗
0,k := lim
n




Therefore, the result follows with the elements of ˜Σj from

lim
n→∞

k,k(cid:48)∈Sj

2j
n

{E[(d∗

0,kd(cid:48)∗

0,k(cid:48))(cid:96),(cid:96)(cid:48)]E[(d∗

0,k(cid:48)d(cid:48)∗

0,k)p,p(cid:48)] + E[(d∗

0,kd(cid:48)∗

0,k(cid:48))(cid:96),p(cid:48)]E[(d∗

0,k(cid:48)d(cid:48)∗

0,k)p,(cid:96)(cid:48)]},
(2.15)

(cid:96), p, (cid:96)(cid:48), p(cid:48) = 1, 2. Note that

Vj E ˜Qn

jVj − Vj EQ jVj = E(dn∗
S2N−j−1


(cid:16)Þ (q+1)2j/2N

Þ (q(cid:48)+1)2j/2N
0,kd(cid:48)n∗

(cid:48)
0,k) − E(d0,kd
0,k)

=

q(cid:48),q=0

q2j/2N

q(cid:48)2j/2N

ψ0,0(t)ψ0,0(s)G j,n(t, s, q, q(cid:48))dtds

(cid:17)

where

and

G j,n = |t − s|H Γ(1, 1)|t − s|H(cid:48) − |(q − q(cid:48))2j/2N|H Γ(1, 1)|(q − q(cid:48))2j/2N|H(cid:48)

(G j,n(t, s, q, q(cid:48)))1,1 < 2(2j/2N)2h1
(G j,n(t, s, q, q(cid:48)))2,2 < 2(2j/2N)2h2
(G j,n(t, s, q, q(cid:48)))1,2 < 2(2j/2N)h1+h2

21

for q2j/2N ≤ t ≤ (q + 1)2j/2N, q(cid:48)2j/2N ≤ s ≤ (q(cid:48) + 1)2j/2N . Therefore

(Vj E ˜Qn
(Vj E ˜Qn
(Vj E ˜Qn

jVj − Vj EQ jVj)1,1 < c(2j/2N)2h1
jVj − Vj EQ jVj)2,2 < c(2j/2N)2h2
jVj − Vj EQ jVj)1,2 < c(2j/2N)h1+h2 .

If j <<

2h1
.5+2h1

log2 n, 2j/2(2j/2N)2h1 → 0, the result is derived.

Analogous to Corollary 2.2.2 is provided below for discrete case.

(2.16)

(2.17)

(2.18)

(cid:3)

Corollary 2.2.9. Let 2j−N → m. ( see Proposition 2.2.1 for m) Then as n → ∞,

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

2j/2

vec(Vj( ˜Qn
2−1/2vec(Vj−1( ˜Qn

j − E ˜Qn
j)Vj)
j−1 − E ˜Qb

j−1)Vj−1)

→d N(0, ˜Σ∗
j),

− E ˜Qn

j−jL

)Vj−jL)

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

is a block matrix where the (m, m)-th block is ˜Σj−m+1, and (m, n)-th block, m > n, consists

where ˜Σ∗
j
of
n→∞ 2(m−n)/22(h(cid:96)
lim

...
j−jL

2−jL/2vec(Vj−jL( ˜Qn


0,kd(cid:48)∗

+E[(d∗

+hp+1)(m−n)

{E[(d∗

2j
n

m−n,k)p,p(cid:48)]

m−n,k(cid:48))(cid:96),(cid:96)(cid:48)]E[(d∗

m−n,k(cid:48))(cid:96),p(cid:48)]E[(d∗

k∈Sj−m+1,k(cid:48)∈Sj−n+1

0,kd(cid:48)∗
0,k(cid:48)d(cid:48)∗

0,k(cid:48)d(cid:48)∗
, p(cid:48) = 1, 2.
(cid:48)
As in the previous case, P−1, P are assumed to be premultiplied by each side of the matrix ˜Qn
j
i . Therefore

and see the convergence rate of the estimator. This is the same as multiplying P−1 by Xn
it is assumed that H = diag(h1, h2). Let
1
2

w(cid:96) log2 ρ( ˜Qn

m−n,k)p,(cid:96)(cid:48)]},

(cid:96))i − 1
2,

˜hi, j =

(cid:96), p, (cid:96)

˜hE
i, j =

1
2

w(cid:96) log2 ρ(E ˜Qn

(cid:96))i − 1
2,

i = 1, 2, where
j

(cid:96)=jL

data is as follows.

w(cid:96) = 0,
j

(cid:96)=jL

(cid:96)w(cid:96) = 1. Now the convergence rate of the estimator of discrete

22

j

j


(cid:96)=jL

(cid:96)=jL

(2.19)

(2.20)

Theorem 2.2.10. If the domain is fixed interval in R, H is diagonizable where h1, h2 ∈ (0, 1) then,

n )2h1,(2j

n )2h1 + (2j

n )2h2,(2j

i, j − hi = O(ci,1
˜hE
ci

( f1(j) − f1(jL)) + O((c + c∗)(cid:48)(f2(j) − f2(jL)))

(2.21)
c(cid:48)
i
ci
is a vector of first derivatives of ci(a, b, d) with respect to a, b, d, and c∗ is a vector of first

for i = 1, 2, where f1(j) = 2−j(2h2−2h1), f2(j) = ((2j
where c(cid:48)
i
derivatives of ci,1
ci
Proof. Set aj, bj, dj from v(cid:48)
jvj as in Proposition 2.2.4. The con-
vergence rates of |aj − an
j| are (2j/2N)2h1,(2j/2N)2h2 and (2j/2N)2h1 +
(2j/2N)2h2, respectively by (2.16-2.18). Therefore similar methods to Propositions 2.2.3,2.2.7,
and by (2.9),(2.10),

2−2jL(h2−h1) with respect to a, b, d.
j from v(cid:48)

j EQ jvj, an

j | and |bj − bn

n )2h2)(cid:48), c =

j |, |dj − dn

j E ˜Qn

j , bn

j, dn

i, j − hE
˜hE

i, j = O((c + c∗)(cid:48)(f2(j) − f2(jL)).

i, j − hi = O( f1(j) − f1(jL)) from Proposition 2.2.7, the result follows.

Since hE
Theorem 2.2.11. As j → ∞, for i = 1, 2,

(cid:3)

2j/2(˜hi, j − ˜hE

i, j) →d N(0, ˜σh),

where ˜σh = W(cid:48) ˜Σρ∗W, and ˜Σρ∗ is analogous with Σρ∗ from ˜Σ instead of Σ.
Remark 1 (1) Op(( f1(j) − f1(jL))) is getting smaller as |h2 − h1| is getting smaller or j is getting
larger.
(2) Op( f2(j) − f2(jL)) is getting smaller as h1 is increasing or j is getting smaller.
(3) If b2
ad
and c−1
2

(cid:170)(cid:174)(cid:174)(cid:172) is close to one, then the bias of ˜h2, j is large, especially for larger j, since f2(j)

are large. Since

b d

in(cid:169)(cid:173)(cid:173)(cid:171)a b
(cid:170)(cid:174)(cid:174)(cid:172) = P(cid:48)c(ψ, H)P
Þ

(cid:169)(cid:173)(cid:173)(cid:171)a b

b d

ψ(t)ψ(s)(cid:169)(cid:173)(cid:173)(cid:171)|t − s|h1

0

=

1
2

(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48)E(XH(1)XH(1)(cid:48))P(cid:169)(cid:173)(cid:173)(cid:171)|t − s|h1

0

(cid:170)(cid:174)(cid:174)(cid:172) dtds,

0

|t − s|h2

0

|t − s|h2

23

where M := P(cid:48)E(XH(1)XH(1)(cid:48))P. This implies that if the
therefore for fixed h1, h2, b2
determinant of the matrix P(cid:48)E(XH(1)XH(1)(cid:48))P is close to zero, the bias of the estimator ˜h2, j will
be large.

M11M22

= c

ad

M2
12

2.2.3 Hurst estimator in discrete noisy data from OFBM

Suppose we observe

i = 1, ..., n, where the ξn
Xt = (Xt,1, Xt,2) is operator fractional Brownian motion. Define the following terms:

i,1, ξn

(2.22)

i,1
ξn
i,2

i,1
Y n
i,2

i,1
Xn
i,2

(cid:170)(cid:174)(cid:174)(cid:172) +(cid:169)(cid:173)(cid:173)(cid:171) ξn
(cid:170)(cid:174)(cid:174)(cid:172) =(cid:169)(cid:173)(cid:173)(cid:171) Xn
(cid:169)(cid:173)(cid:173)(cid:171) Y n
(cid:16)Þ k/2j +(q+1)/2N
(cid:16)Þ k/2j +(q+1)/2N
(cid:16)Þ k/2j +(q+1)/2N

i,2 are centered noise terms. Y is observed in the interval [0, 1],
S2N−j−1

S2N−j−1

S2N−j−1


k2N−j +q,1
ξn
k2N−j +q,2

ψj,k(t)dt

ψj,k(t)dt

k/2j +q/2N

k/2j +q/2N

(cid:170)(cid:174)(cid:174)(cid:174)(cid:172) ,

k2N−j +q

,

(cid:170)(cid:174)(cid:174)(cid:172) ,
(cid:17)Yn
(cid:17)(cid:169)(cid:173)(cid:173)(cid:173)(cid:171) ξn
(cid:17)2(cid:169)(cid:173)(cid:173)(cid:171)σn

(cid:170)(cid:174)(cid:174)(cid:172) .

ρn
1
ρn σn
2

ˆdn
j,k

=

en
j,k

=

vn
j,k

=

q=0

q=0

q=0

(2.23)

(2.24)

(2.25)

Note that ˆdn
j,k

= ˜dn
j,k

+ en

j,k and vn

j,k is the covariance matrix of en
j,k.

k/2j +q/2N

ψj,k(t)dt

ˆdj,k ˆd(cid:48)

j,k − vn
n

j,k

,

w(cid:96) log2 ρ( ˆQn

(cid:96))i − 1
2

for i = 1, 2,

w(cid:96) log2 ρ(E ˆQn

(cid:96))i − 1
2

for i = 1, 2.

ˆQn

j = 

j

j


k∈Sj
1
2

(cid:96)=jL

(cid:96)=jL

ˆhi, j =

ˆhE
i, j =

1
2

Two different assumptions are made for the noise term.

24

ASSUMPTION A.
The noise terms {ξn
= (ξn
i
] = 0, En
(cid:96) = 1, 2, En
H,σ
r, r ∈ [−1, 1] and supi,n En

[ξn
i,(cid:96)

H,σ

i,1, ξn

H,σ

i,2)}i are mutually independent and independent of X. Moreover, for
2) =
)2] = σn
[(ξn
i,(cid:96)
)4] < ∞.
[(ξn
i,(cid:96)

i,2] = ρn, limn σn

= σ(cid:96), limn ρn/(σn

[ξn
i,1ξn

(cid:96) , En

1 σn

H,σ

(cid:96)

ASSUMPTION B.
= (ξn
The noise terms{ξn
i
[(ξn
0. Moreover, En
i,(cid:96)
and c(cid:48)
n, respectively. En
i
supi,n En

[(ξn
i,(cid:96)

H,σ

H,σ

i,2)}i are martingale increments and independent of X, En
i,1, ξn
i/n], En
)2|Fn
[ξn
i,2|Fn
i,1ξn
)4] < ∞ for (cid:96) = 1, 2.

|Fn
i/n] =
i/n] are constants σn
i,p|Fn
(cid:96) , cn
i,(cid:96),
= σ(cid:96), limn ρn/(σn
2) = r, r ∈ [−1, 1] and
1 σn

[(ξn
i,(cid:96)
i/n] = ρn, limn σn

i/n] and En

[(ξn
i,(cid:96)

)3|Fn

[ξn
i,(cid:96)

)2ξn

H,σ

H,σ

H,σ

H,σ

(cid:96)

Assumption A is a special case of Assumption B. Asymptotic distribution of the estimator will

be proved under Assumption B. Note that ˆQn

j − ˜Q j =
3

j

u=1 rn,(u)
=

rn,(2)

j

/n, with
j,k(dj,k)(cid:48)
en

=

=


k

j

j

.

,

k

k

rn,(1)
rn,(3)

j,k)(cid:48)
+ dj,k(en

j,k)(cid:48) − vn
(en
j,k)(en
j,k,
j,k)(cid:48)
j,k)(cid:48)
j,k(bn
j,k(en
+ en
bn
, rn,(3)
Since the expectation of rn,(1)
are zero, E ˆQn
j
estimator ˆhi, j, we get the same result as in Theorem 2.2.10.
Proposition 2.2.12. If the domain is a fixed interval in R, H is diagonizable where h1, h2 ∈ (0, 1)
then, for i = 1, 2,

j . Therefore for the bias of the

, rn,(2)

= E ˜Qn

j

j

j

( f1(j) − f1(jL)(cid:17)

+ O(cid:16)(c + c∗)(cid:48)(f2(j) − f2(jL))(cid:17)

i, j − hi = O(cid:16) ci,1

ˆhE

ci

.

(2.26)

To find the variance of the estimator ˆhi, j, note that
/n + rn,(2)

j + rn,(1)

ˆQn
j = ˜Qn

j

j

/n + rn,(3)

j

/n.

(2.27)

The following result is obtained in a similar way to Proposition 2 in [26].

25

Proposition 2.2.13. Under Assumption B

|( ˆQn

|( ˆQn

j − ˜Qn
j − ˜Qn

+1)),
j)(cid:96),(cid:96)| = Op(n−22j/2 ∨ n−1/22−j(h(cid:96)
j)(cid:96),p| = Op(n−22j/2 ∨ n−1/22−j(h1+1)),

for (cid:96) (cid:44) p, (cid:96), p = 1, 2.

Proof.

E[(rn,(1)

j

)2
(cid:96),p] ≤ c

≤ c

× E

k





(cid:20)(cid:18)(cid:169)(cid:173)(cid:173)(cid:173)(cid:171)

k

n
2j

ξn
k2N−j +q,1
22j
n4 .

≤ cn

j,ken
j,k

E[(en

S2N−j−1


(cid:48) − vn

j,k)2
(cid:96),p]

(cid:16)Þ k/2j +(q+1)/2N

q=0

k/2j +q/2N

(cid:17)4

ψj,k(t)dt

ξn
k2N−j +q,1

2

ξn
k2N−j +q,1

ξn
k2N−j +q,2

ξn
k2N−j +q,2

ξn
k2N−j +q,2

2

(cid:170)(cid:174)(cid:174)(cid:174)(cid:172) −(cid:169)(cid:173)(cid:173)(cid:171)σn

ρn
1
ρn σn
2

(cid:19)2

(cid:21)

(cid:96),p

(cid:170)(cid:174)(cid:174)(cid:172)

j,k)(cid:48) − vn

j,k)(en

The first inequality is derived since (en
second inequality is from the Burkholder-Davis inequality since (en
j,k)(en
martingale increments. The third inequality is derived by the fact that ξn

j,k are uncorrelated when |k − k(cid:48)| ≥ S. The
j,k)(cid:48) − vn
j,k is a sum of
i,l has bounded fourth
/n is

n . Therefore, each element of the matrix rn,(1)

moment andÞ k/2j +(q+1)/2N

ψj,k(t)dt ≤ c 2j/2

j

k/2j +q/2N

of order n−22j/2.

j

j

For rn,(2)

/n, use the Cauchy-Swartz inequality and uncorelatedness of en
/n ≤ 2−jn−1/2(n−HC + Cn−H) and rn,(2)

/n and rn,(3)
j,k when |k−
k(cid:48)| > S, we obtain rn,(3)
/n ≤ 2−j/2n−1/2(C2−j(H+1/2) +
2−j(H+1/2)C). Since 2−jn−1/2(n−HC + Cn−H) ≤ 2−j/2n−1/2(C2−j(H+1/2) + 2−j(H+1/2)C), the
results are derived.

j

j

Theorem 2.2.14. If the domain is a fixed interval in R, H is diagonizable then, for i = 1, 2,

(cid:3)

ˆhi, j − ˆhE

i, j = Op(ci, j,n),

26

where

ci, j,n = Op(n−22j(2hi+3/2) ∨ n−1/22jhi ∨ 2−j/2).

(2.28)

j − E ˆQn

Proof. By (2.3),(2.4), the order of |( ˆQn
E ˆQn
and in (2.4) it converges to a. Set ˜an
The convergence rate of | ˜an
j − ˆan
j| are 2j(2h2+1)|( ˆQn
and | ˜bn
j − ˆbn
therefore | ˆQn
j − E ˆQn
j | = | ˆQn
Proposition 2.2.13.

j)1,1| matters for i = 1 and the orders of |( ˆQn
j −
j)(cid:96),p|, (cid:96), p = 1, 2, matter for i = 2, respectively, since in (2.3) the denominator goes to infinity
j as in Proposition 2.2.4.
j − ˆdn
j |
j, and
j |, the result follows from Proposition 2.2.8 and
(cid:3)

j , ˆbn
j)1,1|, the convergence rates of | ˜dn
= E ˜Qn

j from ˜Qn
j | is 2j(2h1+1)|( ˆQn
j − ˜Qn
j − ˜Qn

j)2,2| and 2j(h2+h1+1)|( ˆQn
j | + | ˜Qn

j)1,2|. Since E ˆQn

j , ˆan
j − ˜Qn

j from ˆQn

j − E ˜Qn

j , ˜bn

j, ˜dn

j − ˜Qn

j, ˆdn

j

Remark 2 Proposition 2.2.12 and Theorem 2.2.14 reveal that with the presence of the noise, the
bias of the estimator ˆhi, j remains the same as in the case when discrete sample path is given without
noise, i.e. the bias of ˆhi, j is the same as that of ˜hi, j . However, the standard error is much larger
as it is given in (2.28), especially when the Hurst indices are big. Theorem 2.2.14 reveals that the
bigger the Hurst indices are, the smaller the j should be chosen in order to reduce the variance of
the estimator. But even so, the performance of the estimator is not good, since for small j, 2−j/2 is
not small. This means that the larger the Hurst indices are, the worse the estimator performs.

In practice, it is hard to calculate vn

is hard to obtain

j,k since either

(cid:16)Þ k/2j +(q+1)/2N
(cid:17)2
(cid:170)(cid:174)(cid:174)(cid:172) is unknown or both. Next, two ways are provided to replace vn
Þ ψ2
j,k(t)dt(cid:169)(cid:173)(cid:173)(cid:171)σn
(cid:169)(cid:173)(cid:173)(cid:171)σn
(cid:170)(cid:174)(cid:174)(cid:172) = 1

ρn
1
ρn σn
2

ψj,k(t)dt

k/2j +q/2N

n

ρn
1
ρn σn
2

j,k . Note that vn

j,k is

(cid:170)(cid:174)(cid:174)(cid:172) as n → ∞.

or(cid:169)(cid:173)(cid:173)(cid:171)σn

ρn
1
ρn σn
2

not changed according to k and it converges to 1
n

27

Replace vn

j,k by vn∗ := 1

n

(cid:170)(cid:174)(cid:174)(cid:172) = 1

n

(cid:169)(cid:173)(cid:173)(cid:171)σn

ρn
1
ρn σn
2

(cid:170)(cid:174)(cid:174)(cid:172) , and define

ρn
1
ρn σn
2
j,k − vn∗
ˆdj,k ˆd(cid:48)
n

,

ˆQn(cid:48)

Þ ψ2
j,k(t)dt(cid:169)(cid:173)(cid:173)(cid:171)σn
j = 

j

j


k∈Sj
(cid:48)
i, j = −1
ˆh
2

(cid:96)=jL

w(cid:96) log2 ρ( ˆQn(cid:48)

(cid:96) )i − 1
2

for i = 1, 2,

for i = 1, 2.

i, j = −1
ˆh(cid:48)E
2

w(cid:96) log2 ρ(E ˆQn(cid:48)

(cid:96) )i − 1
2
is replaced by vn∗, then for i = 1, 2,

(cid:96)=jL

Theorem 2.2.15. If vn
j,k

i, j = Op(ci, j,n)

i, j − hi = O(cid:16) ci,1

ˆh(cid:48)
i, j − ˆh(cid:48)E
ˆh(cid:48)E

( f1(j) − f1(jL)(cid:17)

+ O(cid:16)(c + c∗)(f2(j) − f2(jL))(cid:17)

+ O(2j/n2).

ci

Proof. Since

(2.29)

(2.30)

(2.31)

(2.32)

(2.33)

(cid:3)

(2.35)

(2.36)

(2.37)

/n + rn,(2)

j

/n + rn,(3)

j

/n

and

(cid:17)
j,k − vn∗) + rn,(1)

j

,

j + (vn

j,k − vn∗ = O(cid:16) 2j

ˆQn(cid:48)
j = ˜Qn
vn
E ˆQn(cid:48)

n2

j + O(2j/n2).

j = E ˜Qn
i, j − hi + O(2j/n). Since ˆQn(cid:48)

j

i, j − hi = ˜hE

Therefore ˆh(cid:48)E
random variable, the order of ˆh(cid:48)

Next we define different estimator by using the difference between ˆQn(cid:48)

i, j − ˆh(cid:48)E

i, j is the same as the order of ˆhi, j − ˆhE
i, j .
j − ˆQn(cid:48)
j−1.

= ˆQn

j − vn

j,k

+ vn∗, and vn

(2.34)
j,k − vn∗ is not a

ˆQn(cid:48)(cid:48)
j = ˆQn(cid:48)
(cid:48)(cid:48)
i, j = −1
ˆh
2

i, j = −1
ˆh(cid:48)(cid:48)E
2

j − ˆQn(cid:48)
j

j−1,
w(cid:96) log2 ρ( ˆQn(cid:48)(cid:48)
j


w(cid:96) log2 ρ(E ˆQn(cid:48)(cid:48)

(cid:96)=jL

(cid:96) )i − 1
2

(cid:96) )i − 1
2

(cid:96)=jL

for i = 1, 2,

for i = 1, 2.

28

Note that

j−1 = 


k∈Sj

ˆQn(cid:48)
j − ˆQn(cid:48)

− 


k∈Sj−1

ˆdj,k ˆd(cid:48)
n

j,k

ˆdj−1,k ˆd(cid:48)
n

j−1,k

.

Therefore, ˆh

(cid:48)(cid:48)
i, j is the estimator that can be obtained without knowing vn

j,k or vn∗.

Theorem 2.2.16. For i = 1, 2,
i, j = Op(c(cid:48)

(cid:48)(cid:48)
i, j − ˆh(cid:48)(cid:48)E
ˆh
ˆh(cid:48)(cid:48)E

i, j − hi = O(cid:16) ci,1

ci

i, j,n),

( f1(j) − f1(jL)(cid:17)

+ O(cid:16)(c + c∗)(f2(j) − f2(jL))(cid:17)

+ O(2j+1/n2),

where

c(cid:48)
i, j,n = Op(2n−22j(2hi+3/2) ∨ 2n−1/22jhi ∨ 2−j/2).

(2.38)

Proof. Note that
ˆQn(cid:48)(cid:48)
j = ˜Qn

j − ˜Qn

j−1 + vn

j,k − vn

j−1,k

+ rn,(1)

j

/n − rn,(1)

j−1 /n + rn,(2)

j

/n − rn,(2)

j−1 /n + rn,(3)

j

/n − rn,(3)
j−1 /n

and

by (2.33). Since

E ˆQn(cid:48)(cid:48)

j = E ˜Qn
= E ˜Qn

j − E ˜Qn
j − E ˜Qn

j,k − vn∗ − vn

j−1 + vn
j−1 + O(2j+1/n2)

+ vn∗

j−1,k

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

O((2j/2N)2h1)
O((2j/2N)h1+h2)
O((2j/2N)2h2)

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)
(cid:169)(cid:173)(cid:173)(cid:171)2−j(h1+1/2)

0

(cid:170)(cid:174)(cid:174)(cid:172) ,

0

2−j(h2+1/2)

vec(Vj(E ˜Qn

j − E ˜Qn

j−1)Vj − Vj(EQ j − EQ j−1)Vj) =

and

EQ j − EQ j−1 =(cid:169)(cid:173)(cid:173)(cid:171)2−j(h1+1/2)

0

(cid:170)(cid:174)(cid:174)(cid:172)

(cid:169)(cid:173)(cid:173)(cid:171)a(cid:48) b(cid:48)

b(cid:48) d(cid:48)

(cid:170)(cid:174)(cid:174)(cid:172)

0

2−j(h2+1/2)

29

where a(cid:48) = a(1 − 22h1+1), b(cid:48) = b(1 − 2h1+h2+1), d(cid:48) = d(1 − 22h2+1), using the same method as in
i, j − hi with a(cid:48), b(cid:48)d(cid:48) instead
Proposition 2.2.7 and Theorem 2.2.10, ˆh(cid:48)(cid:48)E
of a, b, d. Since

i, j − hi has the same order as ˜hE

Vj(( ˜Qn
by Corollary 2.2.9, and rn,(1)
twice as much as the order in Proposition 2.2.13 by (2.27), the result for ˆh
same way as in Theorem 2.2.15.

j−1))Vj = OP(2j/2)
j−1 /n + rn,(3)

j−1) − (E ˜Qn
j−1 /n + rn,(2)

j − ˜Qn
/n − rn,(1)

/n − rn,(2)

j − E ˜Qn

j

j

/n − rn,(3)

j

j−1 /n has the order
(cid:48)(cid:48)
i, j − ˆh(cid:48)(cid:48)E
i, j follows the
(cid:3)

by vn∗ increases the bias of the
Remark 3 It is revealed from Theorem 2.2.15 that replacing vn
j,k
i, j, by using ˆQn(cid:48)
estimator by 2j/n2 which is small for small j. From ˆh(cid:48)
i, j to ˆh(cid:48)(cid:48)
j , there is a small
increase in bias but significant increase in standard error as it is seen by comparing Theorem
2.2.15 and Theorem 2.2.16.

2.2.4 Estimator for eigenvectors

i,1, Xn

= (Xn

i,2) = (Xi/n,1, Xi/n,2) is observed for
Assume a discrete sample path of OFBM Xn
i
i = 1, ..., n. In Sections 2.2.1-2.2.3, to calculate the eigenvalues, it was assumed that eigenvector
matrices P, P(cid:48) were premultiplied both sides of Q j, ˜Qn
j . Now, to calculate eigenvectors, we no
j)i), i = 1, 2, be the eigenvector and eigenvalue of the
longer have that assumption. Let (˜un
matrix ˜Qn
j . By ASSUMPTION (OFBM2), one needs to know θ to obtain eigenvectors. To estimate
1, j,2). Since
θ, the eigenvector corresponding to the bigger eigenvalue is chosen, i.e. ˜un
1, j
1, j ≈ (cos θ, sin θ), the esitimator of θ is defined as follows:
˜un
(cid:17)
(cid:16) ˜un

i, j, ρ( ˜Qn

1, j,1, ˜un

= ( ˜un

j, ˆQn

˜θn
j = arctan

1, j,2
˜un
1, j,1

j be the estimator of θ calculated in the same way as above with the matrix E ˜Qn
j .

Let θE
Theorem 2.2.17. As j → ∞,
i)

2j(h2−h1+1/2)( ˜θn

j − θE

j ) →d N(0, σ).

30

ii)

j − θ = O( bn
j
θE
an
j

2j(h1−h2)).

Proof. Note that by Proposition 2.2.8,

j = 2−j(2h1+1)P(cid:169)(cid:173)(cid:173)(cid:171)
j = 2−j(2h1+1)P(cid:169)(cid:173)(cid:173)(cid:171)

˜Qn

E ˜Qn

˜an
j

j2−j(h2−h1)
˜bn

j2−j(h2−h1)
˜bn
j 2−2j(h2−h1)
˜dn
j2−j(h2−h1)
bn
j2−j(h2−h1) dn
j 2−2j(h2−h1)
bn

an
j

j − an

j

= Op(2−j/2), ˜bn

j − bn

j

and ˜an
not change even if a constant is multiplied by the matrix. Therefore ˜θn
from

= Op(2−j/2), ˜dn

j − dn

j

= Op(2−j/2). Eigenvectors of matrix does
j are the ones calculated

j , θE

j = P(cid:169)(cid:173)(cid:173)(cid:171)
j = P(cid:169)(cid:173)(cid:173)(cid:171)

E ˜Qn

j/an

j/ ˜an
˜bn

j/ ˜an
˜Qn

j 2−j(h2−h1)
j/ ˜an
˜bn
1
j 2−2j(h2−h1)
j 2−j(h2−h1)
j / ˜an
˜dn
j 2−j(h2−h1)
j/an
bn
1
j 2−2j(h2−h1)
j 2−j(h2−h1) dn
j /an
j/an
j − E ˜Qn
j/ ˜an
respectively. Therefore the result for i) follows since ˜Qn
j 2−j(h2−h1)
j/an
j 2−2j(h2−h1)
j /an

j 2−j(h2−h1) − bn
j/ ˜an
j 2−2j(h2−h1) − dn
j / ˜an
˜dn

(cid:169)(cid:173)(cid:173)(cid:171)2j(1/2+h2−h1)

2j(1/2+2(h2−h1))

j/an
bn

(cid:170)(cid:174)(cid:174)(cid:172)

0

0

j →d 0 with

(cid:170)(cid:174)(cid:174)(cid:172) →d N(0, Σ).

ii) follows since

with

(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48)
(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48)

(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48)
(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48)

,

,

(cid:169)(cid:173)(cid:173)(cid:171) ˜bn
j − P(cid:169)(cid:173)(cid:173)(cid:171)1 0

0 0

j/an

j/an
j /an
dn

j 2−j(h2−h1)
j 2−2j(h2−h1)

E ˜Qn

(cid:169)(cid:173)(cid:173)(cid:171) bn

(cid:170)(cid:174)(cid:174)(cid:172) P(cid:48) → 0
(cid:170)(cid:174)(cid:174)(cid:172) → 0.

31

(cid:3)

q,q(cid:48)=0

q2j/2N

bn
j
j
bn
j dn
j

S2N−j−1


Þ (q+1)2j/2N
(cid:169)(cid:173)(cid:173)(cid:171)an
(cid:170)(cid:174)(cid:174)(cid:172) =
where Dq,q(cid:48) =(cid:169)(cid:173)(cid:173)(cid:171)|(q − q(cid:48))/2N−j|h1

0

Þ (q(cid:48)+1)2j/2N

q(cid:48)2j/2N

0

|(q − q(cid:48))2N−j|h2

1

2 ψ(t)ψ(s)Dq,q(cid:48)P(cid:48)E(XH(1)XH(1)(cid:48))PDq,q(cid:48)dtds,

(cid:170)(cid:174)(cid:174)(cid:172) .

(2.39)

= cM1,2/M1,1 for

has in terms of both bias and standard error.

Remark 4 (1) The bigger j is, the better precision ˜θn
j
Therefore j = log2 n is the best choice.
(2) The bigger h2 − h1 is, the better the estimator ˜θn
error.
(3) For the standard error of the estimator ˜θ j, j and h2 − h1 play a role, as they increase, the
standard error decreases, whereas for the bias of the estimator ˜θ j, j, h2 − h1 and bn
affect on
the performance. Since by (2.13),(2.14)

is, for it has smaller bias and smaller standard

j/an

j

j

Therefore the matrix M := P(cid:48)E(XH(1)XH(1)(cid:48))P also plays a role, since bn
fixed H, j, n.

j/an

j

2.3 Simulation Results

Operator fractional Brownian motion was simulated with the circulant embedding method,
[12], [28], with t = i/n, i = 1, ..., n, n = 213, θ = .2π. For wavelet function, the second derivative
of Gaussian pdf was used, which has the second vanishing moment. The simulation was repeated
60 times, resulting in 60 independent discrete sample paths of OFBM and (˜h1, j,r, ˜h2, j,r)60
r=1. The
R package “wmtsa" was used for wavelet transform of sample paths. In Table 2.1 and Table 2.2,

C =(cid:169)(cid:173)(cid:173)(cid:171).3 .2

.1 .4

(cid:170)(cid:174)(cid:174)(cid:172) was used, and C =(cid:169)(cid:173)(cid:173)(cid:171).3 .1

.1 .3

(cid:170)(cid:174)(cid:174)(cid:172) , C =(cid:169)(cid:173)(cid:173)(cid:171).1 .3

.2 .4

(cid:170)(cid:174)(cid:174)(cid:172) for Table 2.3 and Table 2.4, respectively.

Table 2.1 to Table 2.4 show the results comparing the eigenvalue method in this chapter with
different sets of j and element-wise method, which uses the first element and the second element
of the paths separately to estimate Hurst indices. The bias and standard error of (˜h1, j,r, ˜h2, j,r)60
r=1
were recorded for different sets of h1, h2 ∈ (0, 1).

32

In Table 2.1, it is seen that the estimation works better when h1, h2 ∈ (.5, 1) than when
h1, h2 ∈ (0, .5). This was expected from Theorem 2.2.10 since f2(j) and f2(j) − f2(jL) are smaller
when h1 > .5. It is also noticeable that, for smaller hi, smaller j works better than when larger j
was used, which comes from the fact that for smaller hi, f2(j) gets larger dominating other term in
(2.2.21), and to make f2(j) smaller one needs to use smaller j so that 2j/n becomes smaller. From
Table 2.2, it is clear that the eigenvalue method works well even when the difference between h2
and h1 is as small as .05.

Table 2.3 shows similar patterns as in Table 2.1 as it shows that, for smaller hi, ˜hj with smaller j
has smaller bias than when larger j was used. For larger hi, ˜hi with larger j performs better, which
comes from the fact that, for large hi, f1(j) needs to be smaller in (2.2.21). This can be seen as a
trade off between f1(j) and f2(j) in (2.2.21) and which term gets larger when hi becomes smaller
or larger.

In Table 2.4, one noticeable difference from previous tables is that the performance of ˜h2 is
significantly worse than previous cases especially when h1 is small. This shows the delicate relation
between h1 and h2 and the dependency within variables. It arises from the fact that, in this case

(cid:170)(cid:174)(cid:174)(cid:172) , c2 is close to zero, which results in large bias in ˜h2 especially when h1 is small

with C =(cid:169)(cid:173)(cid:173)(cid:171).1 .3

large. This phenomena was expected from Theorem 2.2.10 and

.2 .4

since it makes both f2(j) and c−1
2
Remark 1. Note that for (cid:96), p = 1, 2,

M(cid:96),p = (P(cid:48)E[XH(1)XH(1)(cid:48)]P)(cid:96),p = (P(cid:48)C C (cid:48)P)(cid:96),pm(cid:96),p

where

m(cid:96),p =

Þ |eitx − 1|

|x|

|x|−(h(cid:96)

+hp+1)dx.

Therefore for fixed h1, h2, the determinant of M is determined by that of P(cid:48)C C (cid:48)P and with

(cid:170)(cid:174)(cid:174)(cid:172) , det(P(cid:48)C C (cid:48)P) = 4, whereas with C =(cid:169)(cid:173)(cid:173)(cid:171).3 .2

.1 .4

(cid:170)(cid:174)(cid:174)(cid:172) ,(cid:169)(cid:173)(cid:173)(cid:171).3 .1

.1 .3

(cid:170)(cid:174)(cid:174)(cid:172) , det(P(cid:48)C C (cid:48)P) = 100, 64,

C =(cid:169)(cid:173)(cid:173)(cid:171).1 .3

.2 .4
respectively.

33

Throughout Tables 2.1- 2.4, it is seen that with j = 3, jL = .8, the estimator ˜hi, j performs well
in terms of both bias and standard error for any range of hi ∈ (0, 1), whereas for the estimator with
element-wise method, the bias for h1 and standard error for h1, h2 are small as those of ˜hi, j, but the
bias for h2 is close to h2 − h1. This implies that for element-wise method, both the estimators for
h1, h2 are close to h1.

In Table 2.5, the noise term was added to discrete sample paths with t = i/n, i = 1, ..., n, n =

213, C =(cid:169)(cid:173)(cid:173)(cid:171).3 .2

(cid:170)(cid:174)(cid:174)(cid:172) , θ = .2π, and ξn

i

iid∼ N(0, v), v =(cid:169)(cid:173)(cid:173)(cid:171).2 .1

.1 .2

(cid:170)(cid:174)(cid:174)(cid:172) . The simulation was repeated 60 times

.1 .4
i, j,r)60

j , ˆbn

i, j, ˆh(cid:48)(cid:48)

i, j,r)60

r=1,(ˆh(cid:48)(cid:48)

i, j . This is from the fact that the variances of ˆan

to obtain (ˆh(cid:48)
r=1, and the sample mean and the standard error of the estimators were
recorded. Table 2.5 shows the results for hi < .5 and j = 8, jL = 6. Neither the larger set of j nor
larger hi could be used for the estimation of the Hurst indices since in these cases many “NA" were
produced for ˆh(cid:48)
j, ˆdn
j are large as expected in
Proposition 2.2.13, Theorems 2.2.14-2.2.16, making the eigenvalues of ˆQn(cid:48)
j negative which cannot
be used for the logarithm for the estimation, i.e. ci, j,n, c(cid:48)
i, j,n are too large in these cases to estimate
hi. From Table 2.5, it is seen that the means of ˆh(cid:48)
i, j, ˆh(cid:48)(cid:48)
i, j are relatively close to hi when compared
to the estimator ˜hi, j in the presence of noise. With the noise term, the mean of ˜hi, j is far from hi as
seen in the last column in Table 2.5. Also note that the standard error of ˆh(cid:48)
i, j is smaller than that
of ˆh(cid:48)(cid:48)

i, j as it is expected from Theorems 2.2.15-2.2.16.

For Table 2.6, the result was recorded for the eigenvector estimator, ˜θ j, when C =(cid:169)(cid:173)(cid:173)(cid:171).3 .1
(cid:170)(cid:174)(cid:174)(cid:172)

.1 .3
with t = i/n, i = 1, ..., n, n = 213 as before but varying θ ∈ {π/12, π/6, π/4}. As it was predicted
in Theorem 2.2.17, it is seen in Table 2.6 that, the larger the difference between two Hurst indices,
h2 − h1, the better the estimator performs in terms of both bias and standard error. Also, to estimate
θ, it is the best to use the largest j possible, in this case it is 13, as it has the smallest bias and smallest
standard error among the other possible j. It is also noticeable from Table 2.6 that the standard
error is affected by j and h2 − h1, whereas the bias is affected by j, h1 − h2, and θ. From Table
2.6, the bias is getting smaller as θ moves from π/12 to π/4. This is from the fact that the matrix
M = P(cid:48)E(XH(1)XH(1)(cid:48))P is changed and it affects the bias of the ˜θ j as mentioned in Remark 3.

34

For θ = π/12, π/6, π/4, M1,2/M1,1 = 0.217, 0.107, 510−17, respectively.

From the results above, the following estimation methods are recommended.

To estimate h1, h2:
i) When there is no error and the determinant of E(XH(1)XH(1)∗) is not small: The set with larger
j is better to estimate hi when hi is large (e.g. hi ≥ .4) since in this case both bias and standard
error are small. The set with wide range of j is recommended when hi is small, since in this case
standard error can be made small and bias is not too large.
ii) When there is no error and the determinant of E(XH(1)XH(1)∗) is small: For large hi, the set of
large j is recommended for it makes both bias and standard error small. For small hi, either the set
of small j or the set of wide range of j is suggested, but there is a trade off since the first set has
large standard error with smaller bias and the second set has smaller standard error with large bias.
iii) When there is error in the process: It is recommend to use the set of smaller j, but it cannot
be too small. As hi is increasing, the estimation would be worse even with the smaller j, and a
different estimation method is needed.

To estimate θ: Pick largest j possible, i.e. j = log2 n.

35

Table 2.1: Bias and standard error of ˜hi, j

j

13–11

11–8

8–6

6–4

11–6

Element
-wise(11-6)

h1 = .1
h2 = .3

h1 = .2
h2 = .3

h1 = .2
h2 = .4

h1 = .4
h2 = .9

h1 = .25
h2 = .9

h1 = .7
h2 = .8

h1 = .7
h2 = .9

˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2

.264(.02)
.109(.02)

.156(.02)
.103(.01)

.164(.02)
.073(.02)

.063(.02)
.019(.02)

.127(.02)
.029(.02)

.019(.02)
.018(.02)

.020(.02)
.013(.02)

.067(.03)
.025(.04)

.021(.03)
.025(.03)

.021(.03)
.013(.03)

.011 (.03)
-.001(.03)

.028(.03)
.010(.03)

.008(.03)
.002(.03)

.004(.03)
.012(.04)

.032(.10)
.009(.10)

-.003(.09)
.014(.12)

.007(.10)
-.003(.09)

-.004(.10)
.018(.11)

.013(.09)
.040(.10)

-.002(.09)
.021(.06)

.026(.09)
-.014(.09)

-.003(.11)
.021(.13)

.031(.14)
-.012(.13)

-.005(.11)
.022(.16)

-.058(.14)
-.375 (.35)

-.057(.13)
-.484(.33)

-.260(.21)
-.067(.25)

-.269(.20)
-.291(.27)

.050(.03)
.016(.04)

.008(.03)
.020(.04)

.015(.03)
.005(.03)

.054(.03)
.218(.03)

.019 (.03)
.097(.04)

.029(.04)
.199(.03)

.009(.04)
.005(.04)

.009(.04)
.496(.04)

.024(.04)
.020(.03)

.006(.03)
.010(.03)

.013(.03)
.001(.04)

.015(.03)
.657(.03)

.007(.04)
.087 (.04)

.017(.04)
.184(.04)

36

Table 2.2: Bias and standard error of ˜hi, j

j

13–11

11–8

8–6

6–4

11–6

Element
-wise(11-6)

h1 = .2
h2 = .25

h1 = .35
h2 = .4

h1 = .7
h2 = .75

h1 = .85
h2 = .9

.162(.02)
.129(.02)

.086(.02)
.073(.02)

.024(.01)
.016(.02)

.008 (.02)
.012(.02)

.020(.04)
.015(.03)

.005(.03)
.010(.03)

.012(.03)
-.003(.03)

-.005 (.02)
.001 (.03)

.004(.12)
.072(.11)

.011(.10)
.000(.09)

-.010(.11)
.013(.10)

-.004(.11)
.057(.09)

-.028(.10)
.008(.12)

-.074(.15)
.029(.16)

-.156(.22)
-.097(.18)

-.388(.25)
-.167 (.19)

.015(.03)
.032(.04)

.008(.04)
.011(.04)

.003(.03)
.003(.03)

-.004(.04)
.020(.04)

.021(.04)
.052(.03)

.006(.03)
.044(.04)

.005(.04)
.041(.04)

.003(.04)
.040(.04)

˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2

37

Table 2.3: Bias and standard error of ˜hi, j

j

13–11

11–8

8–6

6–4

11–6

Element
-wise(11-6)

h1 = .1
h2 = .3

h1 = .2
h2 = .3

h1 = .2
h2 = .4

h1 = .4
h2 = .9

h1 = .25
h2 = .9

h1 = .7
h2 = .8

h1 = .7
h2 = .9

˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2

.264(.01)
.105(.02)

.161(.02)
.108(.02)

.163(.02)
.067(.02)

.067 (.02)
.013(.01)

.132(.02)
.013(.02)

.026(.02)
.013(.01)

.022 (.02)
.007(.01)

.069(.03)
.012(.03)

.021(.03)
.006(.04)

.021(.04)
.015(.03)

.022 (.04)
.005(.04)

.027(.04)
.007(.03)

-.012(.03)
.004(.03)

.009(.03)
.017(.03)

.039(.13)
.026(.10)

.028(.09)
.021(.10)

.028(.11)
.037(.10)

.002 (.09)
.044(.09)

-.013(.11)
.012(.12)

.030(.11)
.045(.10)

.009(.10)
-.009(.09)

.017(.14)
.007(.11)

-.009(.15)
-.018(.12)

.004(.17)
.032(.16)

-.062(.12)
-.471 (.36)

-.049(.15)
-.463(.33)

-.179(.23)
-.061(.29)

-.286(.20)
-.241(.29)

.053(.04)
.018(.04)

.021(.03)
.011(.04)

.024(.04)
.025(.04)

.018(.04)
.019(.04)

.012(.04)
.009(.05)

-.003(.04)
.020(.04)

.009(.04)
.009(.04)

.040(.03)
.219(.03)

.020(.04)
.098(.03)

.019(.03)
.198 (.04)

.011(.04)
.501(.04)

.018(.04)
.664(.03)

.006(.04)
.088(.04)

.003(.05)
.184 (.04)

38

Table 2.4: Bias and standard error of ˜hi, j

j

13–11

11–8

8–6

6–4

11–6

Element
-wise(11-6)

h1 = .1
h2 = .3

h1 = .2
h2 = .3

h1 = .2
h2 = .4

h1 = .4
h2 = .9

h1 = .25
h2 = .9

h1 = .7
h2 = .8

h1 = .7
h2 = .9

˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2
˜h1
˜h2

.261(.02)
.268(.02)

.159(.02)
.203(.02)

.161(.02)
.243(.01)

.071 (.02)
.159(.02)

.131(.02)
.217(.01)

.022(.02)
.030(.02)

.022(.01)
.053(.02)

.065(.04)
.233(.03)

.013(.03)
.076(.03)

.021(.04)
.113(.03)

.007 (.03)
.036(.03)

.010(.04)
.062(.04)

.002(.04)
-.009(.04)

.012(.04)
.003(.03)

.026(.09)
.105(.11)

.029(.10)
.009(.10)

.011(.11)
.032(.09)

.055 (.12)
.030(.09)

.030(.09)
.021(.10)

-.017(.10)
.035(.10)

-.005(.11)
.017(.10)

.025(.11)
-.032(.19)

-.036(.13)
-.024(.15)

-.015(.10)
-.090(.18)

-.060(.15)
-.769 (.51)

.027(.10)
-.870(.42)

-.229(.26)
-.294(.24)

-.222(.21)
-.752(.42)

.044(.03)
.179(.03)

.015(.04)
.051(.04)

.014(.03)
.078(.03)

.023(.04)
.036(.04)

.011(.04)
.044(.03)

-.006(.03)
.004(.04)

.008(.03)
.007(.04)

.060(.04)
.213(.04)

.035(.03)
.096(.03)

.042(.03)
.199 (.03)

.014(.04)
.493(.04)

.017(.04)
.656(.04)

.017(.04)
.081(.04)

.015 (.03)
.172(.03)

Table 2.5: Mean and standard error of ˆhi, j
j(6-8) Qj − Q j−1 with v(6-8) without v(6-8)
ˆh1
ˆh2
ˆh1
ˆh2
ˆh1
ˆh2

.161 (.10)
.069(.09)

.073 (.09)
.133(.10)

.184(.10)
.382(.17)

.290(.10)
.474(.22)

.219(.09)
.025(.11)

.080(.09)
.265(.12)

.083(.17)
.238(.20)

.179(.17)
.340(.25)

.282(.14)
.510(.32)

h1=.1
h2=.3

h1=.2
h2=.4

h1=.3
h2=.45

39

θ

π/12

π/6

π/4

j

13
11
9

13
11
9

13
11
9

Table 2.6: Bias and standard error of ˜θ j
h1 = .25
h2 = .9

h1 = .4
h2 = .9

h1 = .2
h2 = .3

h1 = .2
h2 = .4

.147(.003)
.189(.010)
.230(.023)

.072(.003)
.092(.009)
.111(.019)

.052 (.001)
.082 (.004)
.113 (.012)

.026 (.001)
.041(.004)
.056(.009)

.004(e-05)
.009(5e-04)
.018(2e-03)

.002(9e-05)
.004(4e-04)
.008(2e-03)

.001(2e-05)
.003(e-04)
.007(8e-04)

4e-04(2e-05)
.001(e-04)
.003(6e-04)

h1 = .1
h2 = .3

.040(.001)
.071(.004)
.108(.011)

.020(.001)
.035(.003)
.051(.009)

h1 = .7
h2 = .8

.172(.004)
.203(.009)
.234(.026)

.084(.004)
.101(.009)
.119(.018)

h1 = .7
h2 = .9

.067 (.002)
.091(.004)
.121(.013)

.033 (.001)
.045 (.004)
.060 (.010)

-6e-05(.001)
-2e-04(.003)
e-03(.009)

3e-04(.003)
e-03(.008)
.002(.016)

-e-04 (.001)
-e-04(.003)
.001(.010)

5e-06(8e-05)
3e-05(4e-04)
6e-05(e-03)

3e-07(2e-05)
9e-06(e-04)
-2e-05(6e-04)

5e-04(.003)
-2e-04(.009)
-.001(.020)

-2e-04 (.001)
-e-04(.004)
-.002(.009)

40

CHAPTER 3

COHERENCE OF MULTIVARIATE RANDOM FIELDS WITH STATIONARY

INCREMENTS

3.1 Preliminaries

Let X = {X(t), t ∈ Rd} be a p-variate random field with X(t) = (X1(t), . . . , Xp(t))(cid:48) ∈ Cp as a
column vector. Throughout we assume that X is a second order random field, that is E[|Xi(t)|2] < ∞
for all t ∈ Rd and 1 ≤ i ≤ p. Here |a| denotes the modulus of a ∈ C.

We say that X has stationary increments (or X is intrinsically stationary) in the weak sense, if

E[X(t) − X(t − r)] = m(r)

for all t, r ∈ Rd, and

E[(X(t + h) − X(t + h − r1))(X(t) − X(t − r2))(cid:48)] = D(h; r1, r2)

for all t, h, r1, r2 ∈ Rd. Here m(r) is called the mean vector of the field X and the matrix-valued
function D(h; r1, r2) is called the structure function of the field X, respectively. We refer to
[27, 29, 52] for historical accounts of the terminology. Here, we will adapt Yaglom’s result for
generalized multivariate random fields to our setting.
Theorem 3.1.1 ([52, Thm. 7]). Let {X(t), t ∈ Rd} be a second order p-variate random field with
stationary increments, mean vector m(r) and structure function D(h; r1, r2). Then there exist a
p × d matrix A with complex entries, a p × p matrix of complex measure F on Rd∗ = Rd \ {0},
and a p × p block matrix U with d × d block matrix Ui j such that

and

Þ

Rd∗

D(h; r1, r2) =

m(r) = A r,

eiλ·h(1 − e−iλ·r1)(1 − eiλ·r2)F(dλ) + Ur1 · r2.

41

Here, each entry Fi j(dλ) of the measure F(dλ) satisfies

Þ

Rd∗

Þ

Rd∗

|λ|2
1 + |λ|2 Fi j(dλ) < ∞,

(3.1)

and F(S) is Hermitian non-negative definite for all Borel set S ∈ Rd∗; the d × d matrices satisfies

d
Ui j = U ∗
, where ∗ denotes the Hermitian conjugate. Let ui j(k, (cid:96)) be the entries of Ui j, then
k,(cid:96)=1 ui j(k, (cid:96))αik α j(cid:96) ≥ 0 for any αik, i = 1, ..., p, k = 1, ..., d.


p

i, j=1

ji

The matrix-valued measure F is called the spectral measure of X which plays a crucial role in
the sequel. There is a general stochastic representation theorem for X. For the sake of simplicity,
we assume throughout that A = 0, U = 0, X(0) = 0. In such a case, the representation of X is
simpler to state.
Theorem 3.1.2 ([52, Thm. 9]). Let {X(t), t ∈ Rd} be a second order p-variate random field with
stationary increments, mean vector m(r) ≡ 0 and structure function D(h; r1, r2) as in Theorem
3.1.1 such that A = 0, U = 0 and X(0) = 0. Then there exists a vector-valued complex random
measure Z(dλ) = (Z1(dλ), . . . , Zp(dλ))(cid:48), such that

X(t) =

(eit·λ − 1)Z(dλ),

(3.2)

where E[Zi(S1)Z j(S2)] = Fi j(S1 ∩ S2) for any Borel S1, S2 ⊂ Rd∗.

Reciprocally, if F , the control measure of Z, satisfies (3.1), then the random field defined as in

(3.2) is a second order random field with stationary increments.

Afterward, Z is referred to as the random spectral measure of X.

3.2 Coherence: definition and basic properties

Similarly to the stationary case considered by Kleiber [30], the coherence between the com-
ponent processes of a multivariate random field with stationary increments can be defined as
follows.

Let Fi j(dλ) be the i j-th entry of the spectral measure F(dλ). Assume that Fi j (cid:28) Leb with

non-vanishing density fi j(λ), where Leb is the Lebesgue measure.

42

Definition 3.2.1. The coherence of component i and component j of the field X at the frequency
λ is defined to be

γi j(λ) :=

(cid:112) fii(λ) fj j(λ) .

fi j(λ)

(3.3)

We remark that, if {Y(t), t ∈ Rd} is a second order stationary p-variate random field with spectral
measure F(dλ) whose entries are necessarily finite measures, by the multivariate extension of
Bochner’s theorem (Cramér [16]), then the random field {X(t), t ∈ Rd} defined by X(t) = Y(t)−Y(0)
has stationary increments with the same spectral measure F(dλ). It follows from Kleiber [30] that
(3.3) is the same as the coherence between the components Yi and Yj of Y.

Let us describe some general properties of the coherence function. Since F(S) is Hermitian
non-negative definite, we know |γ(λ)| ≤ 1. Values of |γ(λ)| near unity indicates strong linear
relationship between Xi and Xj at particular frequency bands.A straightforward way to construct a
legitimate spectral density matrix is the following.
Proposition 3.2.1. Let f ∈ L1(Rd∗, |λ|2(1 + |λ|2)−1dλ) be non-vanishing and C = (ci j) be a p × p
Hermitian non-negative definite matrix. Define fi j(λ) = ci j f(λ). Then the measure F(dλ) =
( fi j(λ)dλ) is Hermitian non-negative definite, and

γi j(λ) =

ci j√
ciicj j

.

When the spectral densities fi j are constructed as products of square integrable functions with

respect to the measure |λ|2(1 + |λ|2)−1dλ, the coherence is constant.
Proposition 3.2.2. For each 1 ≤ i ≤ p, let fi : Rd → [0,∞) be non-vanishing and satisfy

|λ|2
1 + |λ|2 f 2

i (λ)dλ < ∞.

(3.4)

Define fi j(λ) = fi(λ) fj(λ) for 1 ≤ i, j ≤ p. Then the measure F(dλ) = ( fi j(λ)dλ) is Hermitian
non-negative definite, and

γi j(λ) = 1.

43

Þ

Rd∗

Proof. That F defined as such is Hermitian is clear. The non-negative definiteness follows from
the fact that a matrix defined as C = a(cid:48)a for any column vector a ∈ Rp is non-negative definite.
Finally, the integrability condition (3.1) is satisfied by the Cauchy-Schwartz inequality and the
(cid:3)
condition (3.4) on fi.
Remark 1. Let d = 1, p = 2 and f1(λ) = |λ|−1/2−α and f2(λ) = |λ|−1/2−β, 0 < α, β < 1 then the
condition (3.4) is satisfied if and only if 0 < α + β < 1.

Since the Hadamard product of non-negative definite matrices is non-negative definite, we have
the following. Recall that the Hadamard product A ◦ B = (ai j bi j) where ai j, bi j are i j-th entry of
A and B.

Proposition 3.2.3. Let A be a Hermitian non-negative definite matrix and F(dλ) be a Hermitian
non-negative definite matrix of complex measures with spectral densities fi j. Then A ◦ F(dλ) is
non-negative definite and the coherence is

γi j(λ) =

3.3 Examples

(cid:112)aii fii(λ)aj j fj j(λ) .

ai j fi j(λ)

3.3.1 Linear model of coregionalization
In this subsection, we consider the linear model of coregionalization (LMC). Let W = {(W1(t), W2(t))(cid:48), t ∈
Rd} be a second order bivariate random field with stationary increments, A = (ai j) be a 2 × 2
matrix, and define X = A W. Denote the spectral density matrix of W by g = (gi j).

Proposition 3.3.1. The bivariate random field X, defined as above, has stationary increments and
its spectral density matrix f is equal to A gA ∗.
In particular, if the field W has uncorrelated
components, then

f11 = |a11|2g11 + |a12|2g22,

f22 = |a21|2g11 + |a22|2g22

44

and

f12 = f21 = a11a21g11 + a12a22g22.

If we assume further that a11 = a22 = 1, then the coherence between the components of X is of
modulus 1 if and only if a21a12 = 1.
Proof. For any t, h, r1, r2 ∈ Rd, denote X(t) = (X1(t), X2(t))(cid:48), W(t) = (W1(t), W2(t))(cid:48), and let
DX(h; r1, r2) and DW(h; r1, r2) be the structure functions of X and W respectively. We have that

(cid:104)(X(t + h) − X(t + h − r1))(X(t + h) − X(t + h − r1))(cid:48)(cid:105)
(cid:104)(W(t + h) − W(t + h − r1))(W(t + h) − W(t + h − r1))(cid:48)(cid:105) A ∗

DX(h; r1, r2) = E

= A E

= A DW(h; r1, r2)A ∗

,

which completes the proof of the proposition.

(cid:3)

3.3.2 Kernel transform

Now we focus on the effect of convolution to the coherence. Consider a second order complex-
valued random field {X1(t), t ∈ Rd} that has stationary increments in the weak sense. Denote by
f1 and Z1 the spectral density and the random spectral measure of X1, respectively. Assume that
f1(λ) is everywhere non-zero for λ ∈ Rd∗ and let K : Rd → R be a continuous symmetric kernel
that belongs to L2(Rd, Leb). Define

Þ

Rd

X2(t) :=

K(u − t)X1(u)du,

t ∈ Rd .

Proposition 3.3.2. The random field {X2(t), t ∈ Rd} is of second order and has stationary in-
crements. Moreover, {X(t) := (X1(t), X2(t))(cid:48), t ∈ Rd} is a second order bivariate random field
with stationary increments, and its spectral density and random spectral measure are respectively
written as



1

(cid:98)K(λ)

(cid:98)K(λ)
|(cid:98)K(λ)|2

 f1(λ), Z(dλ) = (1,(cid:98)K(λ))(cid:48)Z1(dλ),

f(λ) =

λ ∈ Rd∗,

(3.5)

45

where the Fourier transform,(cid:98)K(λ) =Þ

Rd eit·λK(u)du, is a real valued function. In particular, for

1 ≤ i, j ≤ 2, the coherence functions γi j(λ) = 1 for all λ ∈ Rd∗.
Proof. We first show that {X2(t), t ∈ Rd} has stationary increments. It is enough to calculate the
variogram. Note that, for any s, t ∈ Rd, by Theorem 3.1.2 and the Fubini theorem, we have that

E

Rd

= E

K(u)(X1(u + s + t) − X1(u + s))du

(cid:16)|X2(t + s) − X2(s)|2(cid:17)
(cid:32)(cid:12)(cid:12)(cid:12)(cid:12)Þ
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Þ
(cid:32)Þ
ei(u+t+s)·λ(cid:16)
1 − e−it·λ(cid:17)
= E(cid:169)(cid:173)(cid:171)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)2(cid:170)(cid:174)(cid:172)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Þ
ei(t+s)·λ(cid:16)
1 − eit·λ(cid:17)(cid:98)K(λ)Z1(dλ)
= E(cid:169)(cid:173)(cid:171)
Þ
(1 − cos(t · λ)) |(cid:98)K(λ)|2 f1(λ)(dλ),

K(u)

= 2

Rd∗

Rd∗

Rd

(cid:12)(cid:12)(cid:12)(cid:12)2(cid:33)

Z1(dλ)

(cid:33)

du

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)2(cid:170)(cid:174)(cid:172)

Rd∗

where (cid:98)K is the Fourier transform of K. Hence, {X2(t), t ∈ Rd} has stationary increments and
its spectral density is f2(λ) := |(cid:98)K(λ)|2 f1(λ). Furthermore, we calculate its structure functions as

follows. From the definition and Theorem 3.1.1, we get that

D22(h; r1, r2) = E

(cid:104)(X2(t + h) − X2(t + h − r1))(X2(t) − X2(t − r2))(cid:105)
Þ
Þ
eiλ·(u+h−v)(cid:16)
1 − eiλ·r2(cid:17)
1 − e−iλ·r1(cid:17)(cid:16)
eiλ·h(cid:16)
1 − eiλ·r2(cid:17)(cid:98)K(λ)(cid:98)K(−λ) f1(λ)dλ.

K(u)K(v)D11(u + h − v; r1, r2)dudv
K(u)K(v)

(cid:32)Þ
1 − e−iλ·r1(cid:17)(cid:16)

Rd∗

Rd

Rd

Þ
Þ
Þ

=

=

=

Rd

Rd

Rd∗

(cid:33)

f1(dλ)

dudv

As for the cross structure functions, D12(h; r1, r2) and D21(h; r1, r2), using the same methods,

we can obtain that,

D12(h; r1, r2) = E

(cid:104)(X1(t + h) − X1(t + h − r1))(X2(t) − X2(t − r2))(cid:105)
1 − eiλ·r2(cid:17)(cid:98)K(−λ) f1(dλ),
eiλ·h(cid:16)

1 − e−iλ·r1(cid:17)(cid:16)

Þ

=

Rd∗

46

and

D21(h; r1, r2) = E

(cid:104)(X2(t + h) − X2(t + h − r1))(X1(t) − X1(t − r2))(cid:105)
1 − eiλ·r2(cid:17)(cid:98)K(λ) f1(dλ).
eiλ·h(cid:16)

1 − e−iλ·r1(cid:17)(cid:16)

Þ

=

Because the kernel K is real valued symmetric function, we know that(cid:98)K is a real valued function
and(cid:98)K(λ) = (cid:98)K(−λ). Hence, {X(t), t ∈ Rd} is a bivariate random field with stationary increments

Rd∗

and its spectral density is just the f(λ) in (3.5). Finally, by the definition and Theorem 3.1.2 again,
we know that

Þ
Þ
Þ
Þ

Rd

Rd

Rd∗

Rd∗

X2(t) =

=

=

=

Þ

ei(u+t)·λ − 1

K(u − t)X1(u)du
K(u)
Rd∗
eit·λ − 1
eit·λ − 1

(cid:17)
(cid:16)
Þ
(cid:17)(cid:98)K(λ)Z1(dλ) +
(cid:17)(cid:98)K(λ)Z1(dλ) + X2(0),

Z1(dλ)

(cid:16)
(cid:16)

Rd∗

(cid:16)(cid:98)K(λ) −(cid:98)K(0)(cid:17)

Z1(dλ)

which implies that the random spectral measure of {X(t), t ∈ Rd} is just the Z in (3.5).

(cid:3)

The following result extends Theorem 2 of Kleiber and Nychka [31] and illustrates the role
of coherence in predictive estimation of {X2(t), t ∈ Rd} in terms of kernel smoothed process of
{X1(t), t ∈ Rd}.
Theorem 3.3.3. Suppose that {X(t) = (X1(t), X2(t))(cid:48), t ∈ Rd} is a complex-value zero mean
, λ ∈
bivariate field with stationary increments and its spectral density matrix, f(λ) =
Rd∗, is everywhere non-zero. Let Z = (Z1, Z2)(cid:48) be the random spectral measure of {X(t), t ∈ Rd}.
Then the continuous symmetric square integrable function, K : Rd → R, that minimizes the least

fi j(λ)(cid:17)2

i, j=1

(cid:16)

squares error, E(|X2(t) −Þ
Þ
Þ

1
(2π)d

K(t) =

Rd∗

=

1
(2π)d

Rd∗

Rd K(u − t)X1(u)du|2), is

e−it·λ ·(cid:169)(cid:173)(cid:173)(cid:171) Re
(cid:115)

e−it·λ

(cid:16)(1 − eit·λ) f12(λ)(cid:17)
(cid:32)(cid:16)
1 − eit·λ(cid:17)

f11(λ)

f22(λ)
f11(λ) Re

+(cid:98)K(0) cos(t · λ)(cid:170)(cid:174)(cid:174)(cid:172) dλ
(cid:115)
γ12(λ) + eit·λ(cid:98)K(0)

(cid:33)

dλ.

(3.6)

f11(λ)
f22(λ)

47

Rd K(u − t)X1(u)du, is

(cid:33)(cid:35)2

,

λ ∈ Rd∗ .

(3.7)

E

Re

f11(λ)
f22(λ)

Proof. By the Fubini theorem, we know that,

(cid:34)
(cid:32)(cid:16)
1 − eit·λ(cid:17)
(cid:101)f2|1(λ) = f22(λ)
(cid:32)(cid:12)(cid:12)(cid:12)(cid:12)X2(t) −
Þ
Þ

Furthermore, the spectral density of the predictor,(cid:101)X2(t) :=Þ
(cid:115)
γ12(λ) + eit·λ(cid:98)K(0)
(cid:12)(cid:12)(cid:12)(cid:12)2(cid:33)
(cid:16)
X1(u + t)X1(v + t)(cid:17)
X2(t)X1(u + t)(cid:17)
Þ

T1 := E(X2(t)X2(t)),
T2 :=

K(u)K(v)E
K(u)E

Theorem 3.1.2 yields that,

Rd
T3 = T4 :=

K(u − t)X1(u)du

Þ
Þ
(cid:16)

where,

du.

(cid:17)

(cid:16)

Rd

Rd

Rd

dudv,

(cid:33)

(cid:0)eit·λ − 1(cid:1) Z2(dλ)

:= T1 + T2 − T3 − T4,

T1 = E

(cid:32)Þ
Þ

Rd∗

(cid:16)

Þ

=

= 2

eit·λ − 1

(cid:17)(cid:16)

Z2(dλ) ·
e−it·λ − 1

(cid:17)

Rd∗
eit·λ − 1

Rd∗
f22(λ)dλ

Similarly, by Theorem 3.1.2 and the properties of(cid:98)K, we can obtain that

Rd∗

(1 − cos(t · λ)) f22(λ)dλ.

Þ
Þ
Þ
Þ

Rd

Rd

Rd

Rd∗

T2 =

=

=

=

Rd

dudv

K(u)K(v)E

Þ
X1(u + t)X1(v + t)(cid:17)
(cid:16)
(cid:32)Þ
Þ
(cid:16)
(cid:17)
(cid:32)Þ
Þ
(cid:17)(cid:16)
(cid:16)
(cid:16)((cid:98)K(λ))2 − 2(cid:98)K(λ)(cid:98)K(0) cos(t · λ) + ((cid:98)K(0))2(cid:17)

ei(u+t)·λ − 1

ei(u+t)·λ − 1

K(u)K(v)E

K(u)K(v) ·

Z1(dλ)

Rd∗

Rd∗

Rd

Rd

e−i(v+t)·λ − 1

f11(λ)dλ,

(cid:33)

dudv

Þ

Rd∗

(cid:0)ei(v+t)·λ − 1(cid:1) Z1(dλ)
(cid:17)

(cid:33)

f11(λ)dλ

dudv

48

and

T3 = T4 =

=

=

Þ
Þ
Þ
Þ

Rd

Rd

Rd

du

K(u)E

K(u)E

(cid:16)
X2(t)X1(u + t)(cid:17)
(cid:32)Þ
Þ
(cid:16)
(cid:17)
(cid:32)Þ
(cid:16)
(cid:17)(cid:16)
1 − e−it·λ(cid:17)
(cid:16)(cid:98)K(λ) −(cid:98)K(0)eit·λ(cid:17) ·(cid:16)

eit·λ − 1

eit·λ − 1

Z2(dλ)

K(u) ·

Rd∗

Rd∗

Rd∗

e−i(u+t)·λ − 1

(cid:33)

du

(cid:0)ei(u+t)·λ − 1(cid:1) Z1(dλ)
(cid:17)

(cid:33)

f21(λ)dλ

du

f21(λ)dλ.

=

Rd∗

Note that, the least squares error is a functional of(cid:98)K. So, in order to optimize it, considering the
1 − e−it·λ(cid:17)
(cid:16)

functional derivative, we can get the following equation,

f21(λ) = 2

f21(λ) +

λ ∈ Rd∗,

f11(λ),

(cid:16)

1 − eit·λ(cid:17)
(cid:98)K(λ) =

Re

(cid:16)(cid:98)K(λ) −(cid:98)K(0) cos(t · λ)(cid:17)
+(cid:98)K(0) cos(t · λ).

(cid:16)(1 − eit·λ) f12(λ)(cid:17)

which implies that the minimizer of the least squares error is just

f11(λ)

Finally, the spectral density of the estimator(cid:101)X2(t) follows by Proposition 3.3.2.
(cid:16)
fi j(λ)(cid:17) p
t ∈ Rd and 1 ≤ i ≤ p, define Yi(t) =Þ

It is natural to consider the p-variate random field related to the kernel transformations. Let
X be a p-variate field with stationary increments admitting a spectral density matrix, f(λ) =
, λ ∈ Rd∗, that is everywhere non-zero. Denote its coherence matrix by γf(λ). Let
{Ki, 1 ≤ i ≤ p} be a family of real valued symmetric continuous kernels in L2(Rd, Leb). For
Rd Ki(u − t)Xi(u)du and Y(t) = (Y1(t), . . . , Yp(t))(cid:48). Then, by

i, j=1

(cid:3)

the similar calculations as above, we can get the following result,

Proposition 3.3.4. Under the aforementioned notations, the p-variate random field Y = {Y(t), t ∈
Rd} has stationary increments and its spectral density matrix is given by

(cid:104)(cid:98)Ki(λ) fi j(λ)(cid:98)Kj(λ)(cid:105) p

i, j=1

g(λ) =

,

λ ∈ Rd∗ .

(3.8)

In particular, γg(λ) = γf(λ) for all λ ∈ Rd∗, here γg is the coherence matrix of Y.

49

3.3.3 Estimation of Coherence

To estimate coherence, it is natural to use periodogram. Here we assume that we observe real-valued
bivariate intrinsic stationary process on fixed domain, i.e. Xn(t) = X(t/n) = (X1(t/n), X2(t/n)), t =
1, 2, ..., n. Let X0

n(t) := Xn(t), and

(3.9)

(3.10)

(3.11)

for some k < n, k ∈ Z, and

Then we can estimate the coherence of Xj

n(t) = Xn(t) − Xn(t − 1)
X1
n(t) = Xn(t) − 2Xn(t − 1) + Xn(t − 2)
X2
n(t) = Xj−1
Xj
n


(t) − Xj−1

(t − 1)

n

n

n(t)e−itλ
Xj

D j,n(λ) =

t=1

I j,n(λ) = n−1D j,n(λ)D j,n(λ)∗
12 (λ), by
(cid:113)
12 (λ)
I j,n
11 (λ)I j,n
I j,n

12 (λ) =

n, γ

j,n

j,n

ˆγ

22 (λ) .

for j = {3, ..., k}

for j = {1, ..., k}.

It is easy to see that

I j,n(λ) = I j−1,n(λ)(1 − e−iλ)(1 − eiλ)

and therefore

ˆγ

j,n

12 (λ) =

Lemma 3.3.5.

(cid:113)

12 (λ)
I j,n
11 (λ)I j,n
I j,n

22 (λ)

(cid:113)

=

(λ)

I j+1,n
12
(λ)I j+1,n
I j+1,n
11

22

= ˆγ

j+1,n
12

(λ)

for j = {1, ..., k}.

(λ)

Denote spectral density of {Xj

for j = {1, ..., k}.

ˆγ

j+1,n
12

j,n
12 = ˆγ
n(t), t = 1, ..., n} as f j

n (λ), and its coherence as
, λ ∈ [−π, π].

(cid:113)

j,n

12 (λ) =

γ

n (λ)12
f j
n (λ)11 f j
f j

n (λ)22

50

Remark 5 (1) Note that an intrinsic stationary process and the differenced of the process have the
same coherence. Let ξt = X(t + 1) − X(n) where X(t) satisfies Theorem 3.1.2 with the assumtion
that Fi j (cid:28) Leb with non vanishing density fi j(λ). Then {ξt} is stationary with spectral density :
i j(λ) = (1 − e−iλ)(1 − eiλ) fi j(λ). Therefore coherence of ξ is
f ξ
fi j(λ)(1 − e−iλ)(1 − eiλ)

(cid:112) fii(λ) fj j(λ) = γi j(λ)

fi j(λ)

=

i j(λ) =
γξ

(cid:113)

i j(λ)
f ξ
ii (λ) f ξ
f ξ

j j(λ)

(cid:113)

=

fii(λ) fj j(λ)(1 − e−iλ)2(1 − eiλ)2

(2) The same thing can be said for discretely sampled series in a fixed domain. Let X(t) satisfies
Theorem 3.1.2 and assume that Fi j (cid:28) Leb with non vanishing density fi j(λ), and Xj
are
:= n(1 − e−iλ)j(1 −
defined as (3.9-3.11). Then spectral density of Xj
(cid:96)=−∞ fi j(n(λ + 2π(cid:96)))

(λ)i j := (1 − e−iλ)j+1(1 − eiλ)j+1n
∞

eiλ)j
∞

n and Xj+1
n

n, Xj+1
n

n (λ)i j

are f j

n

(cid:96)=−∞ fi j(n(λ + 2π(cid:96))) and f j+1
(cid:113)

for λ ∈ (−π, π). Therefore
n (λ)12
f j
n (λ)11 f j
f j

12 (λ) =

n (λ)22

(cid:113)

j,n

=

γ

(λ)12
f j+1
n
(λ)22 f j+1

n

f j+1
n

= γ

j+1,n
12

(λ),

λ ∈ (−π, π)

(3.12)

(λ)22

But it is well known that Ii,n is not consistent estimator. Therefore we will use smoothed (cross)
periodogram. Let {X(t), t ∈ R} be a second order bivariate stationary random field satisfying
d = 1, (3.2) and assume that Fi j (cid:28) Leb with non vanishing density fi j(λ) with

f11(λ) ∼ d11|λ|−1−2α1
f22(λ) ∼ d22|λ|−1−2α2
f(cid:96)p(λ) ∼ d(cid:96)p(eiθ1 I(λ>0) + eiθ2 I(λ<0))|λ|−1−2α(cid:96)p

(3.13)

(3.14)

(3.15)

when |λ| → ∞ for some constants d(cid:96)p, θ1, θ2 where 0 < α(cid:96), αp < 1, α(cid:96)p = (α(cid:96) + αp)/2,(α(cid:96)(cid:96) = α(cid:96))
for (cid:96), p = 1, 2. For j ≥ 1

E(I j,n(λ)) =n−1

E

n((cid:96))Xj

e−i((cid:96)−m)λ

n(m)∗(cid:17)

n


(cid:16)Xj

(cid:96),m=1

51

and

E(Xj

n((cid:96))Xj

n(m)∗) =

=

Þ
Þ
ei((cid:96)−m) x
[−π,π] ei((cid:96)−m)x

R

∞


q=−∞

n(e−ix − 1)j(eix − 1)j f(x)dx

n|e−ix − 1|2j f(n(x + 2πq))dx.

Therefore for any fixed λ ∈ [−π, π], λ (cid:44) 0, as n → ∞ the spectral density of Xj

n, f j

n (λ), for

j ≥ 1 has the following limit when max{2α1, 2α2} − 2j < 0, by [32].

∞


q=−∞

n (λ)(cid:96)p ∼ d(cid:96)p
f j

|e−iλ − 1|2j

|n(λ + 2πq)|1+2α(cid:96)p

n

for (cid:96), p = 1, 2.

Therefore we have for (cid:96), p = 1, 2,

∞


q=−∞

n2α(cid:96)p f j

n (λ)(cid:96)p → d(cid:96)p

|e−iλ − 1|2j

|(λ + 2πq)|1+2α(cid:96)p

:= g(λ)j
(cid:96)p

for max{2α1, 2α2} − 2j < 0.

(3.16)

Now let mj = min{ j : max{2α1, 2α2} − 2j < 0} and define a set SJ = { j : 0 ≤ j ≤ c + mj} for
some constant c. Note that with j ∈ SJ, g(·)j

(cid:96)p becomes integrable function. (see [32].)

Define the smoothed cross-periodogram as in [32] by

ˆf j,h(2πn−1J) = 


I∈Tn

Wh(I)I j,n(2πn−1(J + I))

for

j ∈ SJ

(3.17)

where

Wh(I) =


L∈Tn Kh(2πn−1L)

Kh(2πn−1I)

and K is a symmetric continuous function which satisfies Kh(x) = K( x
0, J ∈ Tn, Tn = {(cid:98)−n − 1/2(cid:99), ..., (cid:98)n − n/2(cid:99)}. From [32], it is known that

h)I{|x|≤h}, K(0) > 0, K(x) ≥



nη



n2α1 ˆf j,h
n2α12 ˆf j,h
n2α2 ˆf j,h

11 (2πn−1J1) − g(λ1)j
12 (2πn−1J1) − g(λ1)j
22 (2πn−1J1) − g(λ1)j

11

12

22

...

n2α2 ˆf j,h

22 (2πn−1Jr) − g(λr)j

22

52

→ N(0, Σ)

(3.18)

where

h = Cn−ν, η = (1 − ν)/2, for some C > 0, 0 < ν < 1,
max{2α1, 2α2} − 2j < 0,

(3.19)

(3.20)

and limn→∞ 2πn−1Jr = λr (cid:44) 0. Elements of Σ are

n2ηcov

lim
n

(cid:96)p (2πn−1Jq), n2α(cid:96)(cid:48)p(cid:48) ˆf j,h

(cid:16)



=

n2α(cid:96)p ˆf j,h
{(2π/C) ¯K2/ ¯K2
{(2π/C) ¯K2/ ¯K2
0

1}g(λq)j
1}g(λq)j

(cid:96)(cid:96)(cid:48)g(λq)j
pp(cid:48)
(cid:96)p(cid:48)g(λq)j
(cid:96)(cid:48)p

(cid:96)(cid:48)p(cid:48)(2πn−1Jq(cid:48))(cid:17)

if Jq = Jq(cid:48),
if Jq = −Jq(cid:48),
if λq (cid:44) ±λq(cid:48)

where ¯Kp =Þ

[−1,1] K(x)pdx, and 1 ≤ (cid:96), (cid:96)(cid:48) ≤ p, p(cid:48) ≤ 2. Define

˜γ

j,n

12 (λ) =

(cid:113) ˆf j,h

12 (λ)
ˆf j,h
11 (λ) ˆf j,h

22 (λ)

Lemma 3.3.6. Under (3.19) and limn→∞ 2πn−1J = λ (cid:44) 0

for

j ∈ SJ .

(3.21)

| ˜γ

12 (2πn−1J) − ˜γ
j1,n

12 (2πn−1J)| = O(h) = O(n−ν), λ (cid:44) 0, λ ∈ (−π, π),
j2,n
Proof. Let j1 < j2 and assume j2 satisfies (3.20). From (3.17), for (cid:96), p = 1, 2,

(cid:96)p (2πn−1J)|1 − e−i2πn−1J|2+
(cid:96)p (2πn−1J) = ˆf j1,h
ˆf j2,h
(cid:96)p (2πn−1J)
ˆf j2,h
As h = Cn−ν, 2πn−1K ≤ Cn−ν and therefore

1 − |1 − e−i2πn−1J|2(j2−j1)
|1 − e−i2πn−1(J+K)|2(j2−j1)

(cid:18)

1 − |1 − e−i2πn−1J|2(j2−j1)

|1 − e−i2πn−1(J+K)|2(j2−j1) = O(n−ν),

j1, j2 ∈ SJ .

(cid:19)

.

which results in

(cid:96)p (2πn−1J) = n2α(cid:96)p ˆf j1,h

n2α(cid:96)p ˆf j2,h
(cid:96)p (2πn−1J) converges to g(λ)(cid:96)p as it was shown in [32]. With (3.21), the result is

(cid:96)p (2πn−1J)|1 − e−i2πn−1J|2(j2−j1) + O(n−ν),

since n2α(cid:96)p ˆf j2,h
derived. For the other cases, the result is derived in the similar way.

(cid:3)

53

In (3.16), g(λ)j

12 was defined for j satisfying (3.20). Notice that for any j1, j2 satisfying (3.20),

(cid:113)

g(λ)j2
12
11g(λ)j2
g(λ)j2

22

,

(cid:113)

=

g(λ)j1
12
11g(λ)j1
g(λ)j1
(cid:113)

22

therefore we define

˜γ12(λ) :=
for any j satisfying (3.20). Then for any j ∈ SJ

(cid:113)

j,n

12 (λ) =

γ

n (λ)12
f j
n (λ)11 f j
f j

n (λ)22

=

,

22

g(λ)j
12
11g(λ)j
g(λ)j
nβ12−1 f j∗
n (λ)11nβ22−1 f j∗

n (λ)12

nβ11−1 f j∗

(cid:113)

→ ˜γ12(λ)

n (λ)22

from (3.12) and (3.16) where j∗ satisfies (3.20).
Remark 6 Note that ˜γ12(λ) defined above is slightly different than (3.3). This is because we estimate
coherence in a fixed domain.

Theorem 3.3.7. Let {X(t), t ∈ R} be a second order bivariate stationary random field with d = 1,
(3.1),(3.4), (3.19), and ν > 1/3. For limn→∞ 2πn−1J = λ (cid:44) 0, λ ∈ (−π, π)
j ∈ SJ

12 (2πn−1J) − ˜γ12(λ)) → N(0, Σγ)

nη( ˜γ

j,n

where

(1)
12 (λ),
(1)
12 (λ) is a vector of first derivative of ˜γ12 with respect to g(λ)j

(cid:48)(λ)Σ(λ) ˜γ

Σγ = ˜γ

(1)
12

˜γ

(cid:18)

(cid:113)
− g(λ)j
12
g(λ)j
2
22

(1)
12 (λ) =
˜γ

{g(λ)j

11}−3/2,

12, g(λ)j

22, i.e.

(cid:113)

1
g(λ)j
11g(λ)j

22

11, g(λ)j
(cid:113)
,− g(λ)j
12
g(λ)j
11

2

{g(λ)j

22}−3/2

(cid:19)(cid:48)

and Σ(λ) is from (3.18) for r = 1, λ1 = λ.
Proof. Assume j satisfies (3.20). By (3.18), asymptotic normality of ˜γ
For j which does not satisfy (3.20), Lemma 3.3.6 gives the result since n−ν+η → 0.

j,n

12 (λ) − ˜γ12(λ) is derived.
(cid:3)

Remark 7 The above results show that one doesn’t need to know the spectral density on the tail,
βlp, to estimate coherence since ˜γ

works well as much as ˜γ

2,n
12 .

0,n
12

54

3.3.4 Operator fractional Brownian motion
Let X = {X(t), t ∈ Rd} be a p-variate random field. We say that X is an OFBM with exponent
H if it is a Gaussian field with stationary increments, and satisfies the following operator scaling
property: For any c > 0,

{X(ct), t ∈ Rd} f .d.

= {cHX(t)},

(3.22)

where H is a linear operator on Rd and cH = exp(H ln c) for c > 0. Denote by λH and ΛH the
minimum and the maximum of the real parts of the eigenvalues of H, respectively.

We need some notations. Let Z be a complex-valued Gaussian measure with the control

measure Leb such that, for every Borel set A ⊂ Rd with finite Lebesgue measure,

ReZ(−A) = ReZ(A) and ImZ(−A) = −ImZ(A), a.s.

The construction of such complex valued Gaussian measure is given in Appendix. Let Z1, . . . , Zp
be independent copies of Z and Z = (Z1, . . . , Zp)(cid:48). Note that Z(c ·) f .d.
= cd/2Z(·) for any c > 0.
The following Theorem 3.3.9 is essentially Theorem 3 in [44] and Theorem 3.1 in [39]. Since there
is a slight difference in the representation of [39] and that of Theorem 3.1.2, we shall prove it for
completeness.

We start with a lemma. To simplify the presentation we only consider the case p = 2 in this

subsection. Without loss of generality we assume that H is in its real canonical form.

Lemma 3.3.8. Define the random measure

(cid:18) 1

(cid:19) H+dI/2

Z(dλ),

λ ∈ Rd∗ .

M(dλ) :=

h 0

0 h

|λ|

 with 0 < h < 1, then the spectral density of M is |λ|−(2h+d)I2, where I2 is

i) (H1): If H =

the identity operator on R2.

55

|λ|−(2h2+d)

0

0

P

0

h2

0

0

 ,

0

h2

0

ii) (H2): If H =

|λ|−(2h2+d)

.

h1

0

iii) (H3): If H = P
PP(cid:48) = I, then the spectral density of M is

 with 0 < h1, h2 < 1, then the spectral density of M is
|λ|−(2h1+d)
 .
 P(cid:48) with 0 < h1, h2 < 1, and P = (pi j) is orthogonal matrix, i.e.
h1
|λ|−(2h1+d)
 P(cid:48)
h 1
 with 0 < h < 1, then the spectral density of M is
|λ|−(2h+d)(1 + ln2 1|λ|)
 with 0 < h < 1 and ±β ∈ R (which correspond to the conjugate
h −β
Þ

v) (H5): If H =
eigenvalues h ± iβ), then the spectral density of M is |λ|−(2h+d)I2.
Theorem 3.3.9. Let C be an invertible linear operator on Rd and H be a linear operator on Rd
with 0 < λH, ΛH < 1. Define
X(t) =

|λ|−(2h+d) ln 1|λ|
|λ|−(2h+d)

(3.23)
Then X = {X(t), t ∈ Rd} is a Gaussian field with stationary increments and satisfies the operator
scaling property (3.22) with exponent H.

|λ|

Rd∗

λ ∈ Rd∗ \ {λ ∈ Rd : |λ| ≥ 1};

|λ|−(2h+d) ln 1|λ|

(cid:19) H+dI/2

C Z(dλ),

t ∈ Rd .

iv) (H4): If H =

(cid:18) 1

(eit·λ − 1)

0 h

β

h

Proof. First, it is not hard to check that Yaglom’s condition (3.1) is satisfied for each entry of the
spectral density matrix, the form (3.23) thus implies that X is a second order bivariate field with
stationary increments by Theorem 3.1.2.

56

Second, since e(λ) = eit·λ − 1 is Hermitian in the sense that e(−λ) = e(λ), and the matrix C is
real-valued, we get that the stochastic integral (3.23) is real valued. The Gaussianity follows from
the fact that Z is a Gaussian measure.

Finally, the operator-scaling property (3.22) follows from the scaling property of the measure
(cid:3)

M defined in Lemma 3.3.8, see [39] for details.

Combining (3.23) and Lemma 3.3.8, we can obtain easily the coherence of OFBM.

Proposition 3.3.10. i) Suppose that C = I2. If H =

h 1

0 h

 with 0 < h < 1, then

γ12(λ) =

λ ∈ Rd∗ \ {λ ∈ Rd : |λ| ≥ 1};

,

ln 1|λ|
1 + ln2 1|λ|

(cid:113)
(cid:112)(cid:104)P(cid:48)(λ)e1, P(cid:48)(λ)e1(cid:105)(cid:104)P(cid:48)(λ)e2, P(cid:48)(λ)e2(cid:105),
 ,
 p11 p12|λ|−(h2−h1)

p21 p22|λ|−(h2−h1)

(cid:104)P(cid:48)(λ)e1, P(cid:48)(λ)e2(cid:105)

P(λ) :=

λ ∈ Rd∗,

h 0

0 h

 with 0 < h < 1 or

for H in shape of (H1), (H2) or (H5) in Lemma 3.3.8, we have that γ(λ) = 0. If H is of (H3), then

γ12(λ) =

where

and e1 = (1, 0)(cid:48) and e2 = (0, 1)(cid:48).
ii) Suppose that C = (ci j) is an invertible 2 × 2 matrix.

If H =

β

h

h −β
h1

0

 with 0 < h < 1 and β (cid:44) 0, then
 with 0 < h1 < h2 < 1, then

γ12(λ) =

h2

0

H =

if H =

(cid:112)(cid:104)C (cid:48)e1, C (cid:48)e1(cid:105)(cid:104)C (cid:48)e2, C (cid:48)e2(cid:105),

(cid:104)C (cid:48)e1, C (cid:48)e2(cid:105)

λ ∈ Rd∗;

γ12(λ) =

(cid:112)(cid:104)C (cid:48)(λ)e1, C (cid:48)(λ)e1(cid:105)(cid:104)C (cid:48)(λ)e2, C (cid:48)(λ)e2(cid:105),

(cid:104)C (cid:48)(λ)e1, C (cid:48)(λ)e2(cid:105)

λ ∈ Rd∗,

57



where

if H = P

where

if H =

C(λ) :=

c11

c12

0

h2

γ12(λ) =

h1

0

c21|λ|−(h2−h1) c22|λ|−(h2−h1)

P(λ)e1, C (cid:48)
(cid:104)C (cid:48)
P(λ)e1, C (cid:48)
P(λ)e1(cid:105)(cid:104)C (cid:48)

 P(cid:48) with 0 < h1 < h2 < 1, P = (pi j), then
(cid:113)(cid:104)C (cid:48)
 with 0 < h < 1, then1
h 1
11 + (c11 ln |λ| − c12)2 ·(cid:113)
(cid:113)

0 h
c11c21 ln2 |λ| − (c11c22 + c12c21) ln |λ| + c11c21 + c12c22
21 + (c21 ln |λ| − c22)2
c2

P(λ)e2(cid:105)
P(λ)e2, C (cid:48)

CP(λ) := [PC(λ)P] ;

P(λ)e2(cid:105),

γ12(λ) =

c2

 ;

λ ∈ Rd∗,

λ ∈ Rd∗\{λ ∈ Rd : |λ| ≥ 1}.

,

Remark 8 γ12(λ) tends to 1 as |λ| → 0 or ∞ in either {(H3)} or {(H2, H4) and C (cid:44) I}.

Let {X(t), t ∈ [0, 1]} be operator fractional Brownian motion observed in a fixed domain whose

spectral representation is (3.23) with d = 1. For Xj

∞


q=−∞

n (λ) =
f j

|e−iλ − 1| j

|n(λ + 2πq)|H+I/2

n

n defined as (3.9-3.11), its spectral density is
C C (cid:48)

|e−iλ − 1| j

j ∈ SJ .

|n(λ + 2πq)|H+I/2,

(3.24)

We further restrict H as one of (H1), (H2), and (H3) for which we consider the following cases.

Case I: C is diagonal matrix and (H1) or (H2)
Case II: C is not diagonal and (H1) or (H2)
Case III: (H3)

For (H3) and C = I,

∞


q=−∞

n (λ) = P
f j

|e−iλ − 1|2j
|n(λ + 2πq)|2Λ+I P(cid:48)

,

n

1Need to be checked again.

58

(3.25)

h1

0

where Λ =

0

h2

 . As in proposition 3.3.10 i), we have the following:
(cid:113)(cid:104)P j,n(cid:48)(λ)e1, P j,n(cid:48)(λ)e1(cid:105)(cid:104)P j,n(cid:48)(λ)e2, P j,n(cid:48)(λ)e2(cid:105),
p11 p12q j,n(λ)

(cid:118)(cid:117)(cid:116) ∞


(cid:104)P j,n(cid:48)(λ)e1, P j,n(cid:48)(λ)e2(cid:105)

p21 p22q j,n(λ)

(cid:30) ∞


|e−iλ − 1|2j

|e−iλ − 1|2j

n

|n(λ + 2πq)|2h2+1

q=−∞

n

|n(λ + 2πq)|2h1+1

q=−∞

j,n

12 (λ) =

γ

P j,n(λ) =

q j,n(λ) =

λ ∈ (−π, π)

(3.26)

(3.27)

(3.28)

Now we are ready to derive the following result. For notational simplicity, q(λ) means q j,n(λ).

Lemma 3.3.11. Let X(t) be operator fractional Brownian motion defined as in Theorem 3.3.9 with
Case III, observed in a fixed domain. For λ (cid:44) 0, λ ∈ (−π, π)
i) C = I and (H3)

j,n

12 (λ)| − 1 = O(n2(h1−h2)),
|γ

j ∈ SJ .

ii) C (cid:44) I and (H3)

(cid:18) b12

b11

nh1−h2 ∨ b22
b11

n2(h1−h2)(cid:19)

,

j ∈ SJ,

|γ

j,n

12 (λ)| − 1 = O

where B = (bi j) = P(cid:48)C C (cid:48)P.
Proof. i) By remark 5 (2), let j = 2 without loss of generality. Since h1, h2 ∈ (0, 1),

|e−iλ − 1|2j
|(λ + 2πq)|2h1+1 (cid:16) 1,

λ (cid:44) 0, λ ∈ (−π, π).

∞


|e−iλ − 1|2j

q=−∞

|(λ + 2πq)|2h2+1

q=−∞
Therefore q(λ) (cid:16) nh1−h2 . By (3.26-3.28),

(cid:30) ∞


(cid:113)(p2

j,n

12 (λ) =

γ

p11p21 + p12p22q2(λ)
21 + p2

12q2(λ))(p2

11 + p2

22q2(λ)) .

59

ii) Similarly, from (3.24)

(cid:114)

j,n

12 (λ) =

γ

p11p21 + b12
b11
p11p12q(λ) + b22
b11

p12p22q2(λ)
(p12p21 + p11p22)q(λ) + b22
b11
p21p22q(λ) + b22
21 + 2 b12
b11
b11

12q2(λ))(p2
p2

(p2
11 + 2 b12
b11

.

22q2(λ))
p2

(cid:3)

j,n

12 (λ) ≡ 0 for Case I.

, C(cid:96),k = (C C (cid:48))(cid:96),k, (cid:96), k = 1, 2. This is from (3.24)

j,n

C12√
C11C22

(cid:118)(cid:117)(cid:116) ∞


Remark 9 (1) In a similar way to proposition 3.3.10, γ
12 (λ) (cid:16) cγ where cγ =
(2) For Case II, γ
∞

and the fact that for λ (cid:44) 0, λ ∈ (−π, π),
|e−iλ − 1|2j
|(λ + 2πq)|2h1+1 (cid:16)
C12
∞

q=−∞
since from (3.24),

|(λ + 2πq)|2h2+1

|e−iλ − 1|2j

q=−∞

∞


q=−∞

j,n

12 (λ) =

γ

(cid:114)
C11C22
∞
12 (λ) ≡ cγ.

j,n

q=−∞

|e−iλ−1|2j

|(λ+2πq)|h1+h2+1


∞

q=−∞
|e−iλ−1|2j
|(λ+2πq)|2h2+1

q=−∞

.

|e−iλ−1|2j
|(λ+2πq)|2h1+1

Especially for (D1), γ
Let us continue to assume that (H3) and C = I hold, since results for other cases follow

|e−iλ − 1|2j

|(λ + 2πq)|h1+h2+1,

similarly. Let

n


n(t) = P(cid:48)Xj
n(t),
Yj
n(t)e−itλ,
Y (λ) =
Yj
D j,n
Y (λ) = n−1D j,n(wj)D j,n(λ)∗
I j,n

t=1

,

j ∈ SJ .

Define the smoothed cross-periodogram as in (3.17) by

Y (λ) = 


ˆf j,h

K∈Tn

Wh(K)I j,n

Wh(K)I j,n(λ).

ˆf j,h(λ) = 


K∈Tn

By (3.25), it is easy to see that the spectral density of Y j

n is

Y (λ) =
f j

|e−iλ − 1|2j

n

|n(λ + 2πq)|2Λ+I

Y (λ)
∞


q=−∞

60

and

ˆf j,h(λ) = P ˆf j,h

Y (λ)P(cid:48)

.

(3.29)

Since the spectral density of Y satisfies (3.13-3.15) with α1 = h1, α2 = h2, d11 = d22 = 1, d12 =
d21 = 0, the result (3.18) is applied for ˆf j,h

Y (λ),



nη

nh1 ˆf j,h

Y,11(λ1) − gY(λ1)j

11

...

nh2 ˆf j,h

Y,22(λr) − gY(λr)j

22

 → N(0, ΣY),

∞


q=−∞

|e−iλ − 1|2j
|(λ + 2πq)|2Λ+I

where

gY(λ)j =

and ΣY is defined as in (3.18) except that g is replaced by gY . The coherence estimator of X is
defined as

and, by (3.29), ˆf j,h(λ) converges to

(3.30)

(3.31)

(3.32)

(cid:113) ˆf j,h
 gY(λ)j

j,n

12 (λ) :=
˜γ

12 (λ)
ˆf j,h
11 (λ) ˆf j,h

22 (λ),

P

.

0

0

0
n−2h2

0
n−2h2

n−2h1

 P(cid:48)
n−2h1
(cid:112)(cid:104) ˜P(cid:48)(λ)e1, ˜P(cid:48)(λ)e1(cid:105)(cid:104) ˜P(cid:48)e2, ˜P(cid:48)(λ)e2(cid:105) = ˜γ12(λ),
p11 p12q(λ)

 , q(λ) =

n−2h2gY(λ)j
22
n−2h1gY(λ)j
11

(cid:104) ˜P(cid:48)(λ)e1, ˜P(cid:48)(λ)e2(cid:105)

p21 p22q(λ)

.

j,n

12 (λ) =

γ

˜P(λ) =

Note that by (3.16) and (3.24), in OFBM, γ

j,n

12 (λ) = ˜γ12. By (3.25-3.29) and gY defined above,

where

Theorem 3.3.12. Let X(t) be operator fractional Brownian motion defined as in Theorem 3.3.9.
For limn→∞ 2πn−1J = λ (cid:44) 0, λ ∈ (−π, π), j ∈ SJ with
Case I:

nη( ˜γ

j,n

12 (λ) − ˜γ12(λ)) → N(0, Σγ),

61

where

Case II:

˜γ12(λ) ≡ 0.

nη( ˜γ

j,n

12 (λ) − ˜γ12(λ)) → N(0, Σγ),

where ˜γ12(λ) (cid:16) cγ and

(cid:112)(C C (cid:48))1,1(C C (cid:48))2,2

(C C (cid:48))1,2

.

cγ =

Case III: i) C = I and (H3)

nη+2(h2−h1)( ˜γ

j,n

12 (λ) − ˜γ12(λ)) → N(0, Σγ),

where

ii) C (cid:44) I and (H3)

where

| ˜γ12(λ)| − 1 = O(n2(h1−h2)).

nη+τ( ˜γ

j,n

12 (λ) − ˜γ12(λ)) → N(0, Σγ),

| ˜γ12(λ)| − 1 = O(τ(n, C , h2, h2)) and τ(n, C , h2, h2) =

b12
b11

nh1−h2 ∨ b22
b11

n2(h1−h2)

.

Proof. For CaseI and II the results follow from Theorem 3.3.7 and Remark 9.
For CaseIII − i), by (3.29), ˜γ

12 (λ) can be written in the same way as in (3.26) with

j,n

p11 p12 ˜q j,n(λ)
 ,
(cid:118)(cid:117)(cid:117)(cid:117)(cid:116) ˆf j,h

p21 p22 ˜q j,n(λ)
Y,22(λ)
Y,11(λ) .
ˆf j,h

˜Pi,n(λ) =

˜q j,n(λ) =

From (3.30), (3.32), it follows that

nη+2(h2−h1)( ˜q j,n(λ) − q(λ)) → N(0, σ2
q),

62

where

(cid:113)

q = q(1)(cid:48)ΣY(λ)q(1)
σ2

,

gY(λ)j

22{gY(λ)j

11}−3/2, 1/2{gY(λ)j

q(1) = (−1/2
11}−1/2)(cid:48), and ΣY(λ) is from (3.30)
with r = 1, λ1 = λ. Therefore by Lemma 3.3.11, and (3.31-3.32), the result is obtained with
(1)
12 (λ) is the first derivative of ˜γ12 with q(λ). In CaseIII − ii), the result
Σγ = (σq ˜γ
(cid:3)
follows in the same way.

(1)
12 (λ))2, where ˜γ

22}−1/2{gY(λ)j

Remark 10 (1) As the coherence reflects the correlation between X1 and X2 on the spectral domain,
depending on C and H, the dynamics of X1 and X2 have different features. In Case I, the coherence
function is zero, which implies that it has zero correlation, and in Case II the coherence function
depends on C C (cid:48). In Case III, the coherence is getting close to either 1 or -1 as sample size n is
increasing on a fixed domain, and the larger the difference in Hurst indices h2 − h1 is, the closer
the coherence is to 1 or -1, which implies that X1 and X2 have very strong correlation.
12 (λ) = ˜γ12. In Case I and II, ˜γ12(= γ
(2) Note that by (3.16) and (3.24), in OFBM, γ
function of j, λ, but not n. However, in Case III, ˜γ12(= γ
according to sample size n. More specifically, it increases to one as the sample size grows.

12 (λ)) is a
12 (λ)) is a function of n, j, λ, and it changes

j,n

j,n

j,n

3.4 Simulation and estimation of spectra

In this section, we provide simulation results on the squared coherence function of operator
fractional Brownian motion discussed in Section 3.3.4.
In practice, squared of the coherence
function is more often used as it is real-valued in [0, 1], and if its value is near 1 that indicates strong
linear relationship between X1 and X2 at particular frequency bands.

3.4.1

Independent sample paths: Case I

In this case, where C and H are diagonal matrices, two processes are independent. As such, the
coherence is zero in every frequency. We set H as diagonal matrix with h1 = .4, h2 = .8, C = I, and
simulate the sample paths with t = i/n, i = 0, .., n, n = 1000. As the sample paths are independent,

63

Figure 3.1: sample paths X1 and X2

2,n

each of them is governed by the Hurst index h1, or h2 as seen in Figure 3.1. Figure 3.2 shows
the estimates of squared coherence of X1 and X2, ( ˜γ
12 )2, and its theoretical coherence
function ( ˜γ12)2. The result in Figure 3.2 reflects well Theorem 3.3.12, ii). The estimates of
squared coherence, ( ˜γ
n, and two times differenced series X2
n,
respectively, are close to each other except on the frequency near zero and they are close to the
theoretical coherence function ( ˜γ12)2, which is zero function.

12 )2, from original series X0

12 )2,( ˜γ

0,n

0,n

12 )2,( ˜γ

2,n

Figure 3.2: The estimates of squared coherence of X1 and X2

64

3.4.2 Correlated sample paths with diagonal H: Case II

If either C or H is not diagonal matrix, then sample paths are correlated. Let us first assume that
H is diagonal but C is not, which is Case II. Then, as H is diagonal, the roughness of sample paths
are still governed by either h1 or h2. But the two processes are not independent anymore. Figure

3.3 shows the simulated sample paths with C =(cid:169)(cid:173)(cid:173)(cid:171)1 3

t = i/n, i = 0, .., n, n = 1000.

2 1

(cid:170)(cid:174)(cid:174)(cid:172) and H =(cid:169)(cid:173)(cid:173)(cid:171)h1

0

(cid:170)(cid:174)(cid:174)(cid:172) , h1 = .4, h2 = .8 and

0

h2

Figure 3.3: sample paths of X1 and X2

Figure 3.4: The estimates of squared coherence of X1 and X2

The estimates of squared coherence, ( ˜γ

0,n

12 )2,( ˜γ

2,n

12 )2, are in Figure 3.4. The estimates of squared

65

coherence are close to each other except on the frequency near zero and are close to the theoretical

squared coherence function which is close to c2
γ

from Theorem 3.3.12,iii).

= .5, since C C (cid:48) =(cid:169)(cid:173)(cid:173)(cid:171)10 5

5

5

(cid:170)(cid:174)(cid:174)(cid:172), as it was expected

3.4.3 Correlated sample paths with diagonizable H: Case III

Here, we only consider the case where C = I but H is not diagonal, but diagonizable matrix,

H = PΛP(cid:48), where PP(cid:48) = I and Λ =(cid:169)(cid:173)(cid:173)(cid:171)h1
Figure 3.5 shows the sample paths with h1 = .4, h2 = .8, C = I, P =(cid:169)(cid:173)(cid:173)(cid:171)cos θ

(cid:170)(cid:174)(cid:174)(cid:172) . The two sample paths are dependent, but their depen-
(cid:170)(cid:174)(cid:174)(cid:172) , θ = .2 ∗ π,

dency is quite different from Case II and shows different behaviour in sample paths and coherence.

sin−θ
cos θ

with t = i/1000, i = 0, 1, .., 1000.

sin θ

0

0

h2

Figure 3.5: sample paths of X1 and X2

The roughness of two sample paths are similar since both paths are governed by h1 = .4. From
Figure 3.6, it is seen that the estimates of squared coherence are very close to each other except on
frequency zero, and the theoretical coherence function is very close to 1 as it was expected from
Theorem 3.3.12. As h1, and h2 are getting close to each other, the theoretical coherence function
12 )2, are getting smaller. This is observed well in
˜γ12 and the squared of its estimates, ( ˜γ

12 )2,( ˜γ

2,n

0,n

66

Figure 3.7, 3.8, and this results were predicted from Theorem 3.3.12.

Figure 3.6: θ = .2 ∗ π, h1 = .4, h2 = .8, coherence of X1,X2

Figure 3.7: θ = .2 ∗ π, h1 = .65, h2 = .8,

67

Figure 3.8: θ = .2 ∗ π, h1 = .65, h2 = .8,

68

CHAPTER 4

ESTIMATING HURST INDICES IN OPERATOR SCALING GAUSSIAN RANDOM

FIELD

4.1 Estimation method

Let X be an operator scaling Gaussian random field that has stationary increment and operator
scaling property (1.7). Its semi-variogram is defined as (1.8) and it has the following property:
υE,H(cE h) = c2HυE,H(h). Also, note that if X is operator scaling for E and H > 0, it is also
for E/H and 1 and (E, H) is not uniquely defined. In the following sections, eigenvectors of E,
u1, u2, and the ratios of H to eigenvalues of E, H/λ1, H/λ2, are estimated where E = UΛU(cid:48), U =
(u1, u2), Λ = diag(λ1, λ2).
Eigenvectors u1, u2, are estimated first, and H/λ1, H/λ2, are computed by using the relation

{X(cλiuit)}t∈R

f .d

= {cH X(uit)}t∈R for any c > 0,

(4.1)

which means that {X(uit)}t∈R is self-similar with Hurst index H/λi, i = 1, 2.
Different sampling regimes are considered in this chapter.
In Section 4.1, it is assumed that
independent sample paths are obtained from OSGRF in a fixed domain or increasing domain. Later
in simulation part, Section 4.2, the performance of the estimators is investigated with X observed
on a grid in a fixed domain, {X(i/n, j/n), i, j = 0, 1, .., n}, in addition to the independent sample
paths case. Throughout this chapter, notation E is used either for scaling matrix E or expected
value, but with the context, there should be no confusion on what it indicates.

4.1.1 Estimation of eigenvectors of E

We make some assumptions on the matrix E and semi-variogram υE,H.

ASSUMPTION

69

(A1) υE,H(t) = τE(t)2H with τE in (1.8), H = 1 and a1 ≤ a2, a1, a2 ∈ (0, 1).
(A2) The matrix E in (1.7) is diagonizable with eigenvalues λ2 (cid:44) λ1.

By (A2), we have the following:

E =(cid:169)(cid:173)(cid:173)(cid:171) cos θ

− sin θ

(cid:170)(cid:174)(cid:174)(cid:172)

(cid:169)(cid:173)(cid:173)(cid:171)λ1

0

sin θ
cos θ

Since H = 1,

(cid:170)(cid:174)(cid:174)(cid:172) ,

θ ∈ [0, π).

0

(cid:169)(cid:173)(cid:173)(cid:171)cos θ − sin θ
(cid:170)(cid:174)(cid:174)(cid:172)
= { 2


cos θ

sin θ

λ2

f .d.

{X(t)}t∈R2

(4.2)
where u1 = (cos θ,− sin θ), u2 = (sin θ, cos θ) and Ba1, Ba2 are two independent fractional Brownian
motion with semi-variogram || · ||2a1, || · ||2a2, respectively.
Define for i = 2, .., n,

i=1

Bai((cid:104)t, ui(cid:105))}t∈R2,

(4.3)

(4.4)

(4.6)

(4.7)

(4.8)

Note that {∇1X(i/n), ∇2X(i/n), i = 2, .., n} are (covariance) stationary processes for fixed n, and
(4.5)

=

where ∇Baj(i) = Baj(i) − 2Baj(i − 1) + Baj(i − 2),
In an increasing domain, the followings are defined analogous to (4.5-4.6), which is used in later
Section 4.1.2. For i = 2, .., n,

0

0

0

(cid:170)(cid:174)(cid:174)(cid:172) − 2X(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171)(i − 2)/n
∇1X(i/n) = X(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:170)(cid:174)(cid:174)(cid:172) ,
(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171)
(cid:170)(cid:174)(cid:174)(cid:172) − 2X(cid:169)(cid:173)(cid:173)(cid:171)
∇2X(i/n) = X(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:170)(cid:174)(cid:174)(cid:172) .
(cid:12)(cid:12)(cid:12)(cid:12)cos θ
(cid:12)(cid:12)(cid:12)(cid:12)a2∇Ba2(i),
(cid:12)(cid:12)(cid:12)(cid:12)sin θ
(cid:12)(cid:12)(cid:12)(cid:12)a1∇Ba1(i) +
(cid:12)(cid:12)(cid:12)(cid:12)a2∇Ba2(i),
(cid:12)(cid:12)(cid:12)(cid:12)cos θ
(cid:12)(cid:12)(cid:12)(cid:12)a1∇Ba1(i) +
(cid:12)(cid:12)(cid:12)(cid:12)sin θ

∇1X(i/n) f .d.
∇2X(i/n) f .d.

(i − 2)/n

(i − 1)/n

i/n

0

0

n

n

n

=

n
j = 1, 2.

∇1X(i) = X(cid:169)(cid:173)(cid:173)(cid:171)i
∇2X(i) = X(cid:169)(cid:173)(cid:173)(cid:171)0

0

i

(cid:170)(cid:174)(cid:174)(cid:172) − 2X(cid:169)(cid:173)(cid:173)(cid:171)i − 1
(cid:170)(cid:174)(cid:174)(cid:172) − 2X(cid:169)(cid:173)(cid:173)(cid:171) 0

i − 1

0

(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171)i − 2
(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171) 0

i − 2

0

(cid:170)(cid:174)(cid:174)(cid:172) ,
(cid:170)(cid:174)(cid:174)(cid:172) .

70

0

i) υE,H

(cid:169)(cid:173)(cid:173)(cid:171)h
iii) limn
n

Proof. i) Since

ii) E(∇1X(i/n)2) = 8υE,H

the results are derived.
ii)

E(∇1X(i/n)2) = E

{∇1X(i), ∇2X(i), i = 2, ..., n} are stationary processes for fixed n, and

∇1X(i) f .d.
∇2X(i) f .d.

= |cos θ|a1∇Ba1(i) + |sin θ|a2∇Ba2(i),
= |sin θ|a1∇Ba1(i) + |cos θ|a2∇Ba2(i).

(4.10)
For the rest of the chapter, we restrict θ ∈ [0, π/2), since if θ ∈ [π/2, π), then we can just switch
{∇1X(i/n), ∇1X(i)} to {∇2X(i/n), ∇2X(i)}, and have the same results.
Lemma 4.1.1. Under the assumptions (A1, A2),

(4.9)

(cid:170)(cid:174)(cid:174)(cid:172) .

h

0

0

1/n

2/n

j = 1, 2.

(cid:169)(cid:173)(cid:173)(cid:171) 0

i=2 cov(∇j X(2)2, ∇j X(i)2) < ∞,

(cid:169)(cid:173)(cid:173)(cid:171)0
(cid:170)(cid:174)(cid:174)(cid:172) ∼ |h sin θ|2a1 for |h| → 0 and for fixed θ ∈ [0, π/2).
(cid:170)(cid:174)(cid:174)(cid:172) ∼ |h cos θ|2a1 and υE,H
(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:169)(cid:173)(cid:173)(cid:171)1/n
(cid:169)(cid:173)(cid:173)(cid:171)2/n
(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) , E(∇2X(i/n)2) = 8υE,H
(cid:169)(cid:173)(cid:173)(cid:171)h
(cid:170)(cid:174)(cid:174)(cid:172) =|h cos θ|2a1 + |h sin θ|2a2,
(cid:169)(cid:173)(cid:173)(cid:171)0
(cid:170)(cid:174)(cid:174)(cid:172) =|h sin θ|2a1 + |h cos θ|2a2,
(cid:18)
(cid:19)2
(cid:19)2
X(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
X(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:170)(cid:174)(cid:174)(cid:172) − X(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
(cid:170)(cid:174)(cid:174)(cid:172) − X(cid:169)(cid:173)(cid:173)(cid:171)(i − 2)/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:20)(cid:18)
(cid:19)(cid:18)
(cid:19)(cid:21)
X(cid:169)(cid:173)(cid:173)(cid:171)i/n
− X(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
(cid:170)(cid:174)(cid:174)(cid:172) − X(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171)(i − 2)/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:169)(cid:173)(cid:173)(cid:171)1/n
(cid:169)(cid:173)(cid:173)(cid:171)1/n
(cid:169)(cid:173)(cid:173)(cid:171)1/n
(cid:169)(cid:173)(cid:173)(cid:171)2/n
(cid:169)(cid:173)(cid:173)(cid:171)0
(cid:170)(cid:174)(cid:174)(cid:172) + 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) + 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) + 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) .

(cid:169)(cid:173)(cid:173)(cid:171)1/n

0

h

0

0

0

0

0

0

υE,H

υE,H

0

0

0

0

+ 2E

0

+ E

(cid:18)

0

= 2υE,H

0

0

71

iii) are derived by (4.7-4.10). Let j = 1. For j = 2, it is derived in the same way.

cov(∇j X(2)2, ∇j X(i)2) = 2 lim

cov(∇j X(2), ∇j X(i))2

n


i=2

n

(| cos θ|2a1cov(∇Ba1(2), ∇Ba1(i)) + | sin θ|2a2cov(∇Ba2(2), ∇Ba2(i)))2

n


i=2

lim
n

= 2 lim
n

n


i=2

and the fact that cov(∇Baj(2), ∇Baj(i)) ∼ ci2aj−4 for some c ∈ R, thus
cov(∇Baj(2), ∇Baj(i))2 < ∞, j = 1, 2.

n


lim
n

i=2

Define

Pn :=

∇1X(i/n)2

n − 1

,

Qn :=

Lemma 4.1.2. Assumptions (A1, A2) are satisfied. Then,

n


i=2

n


i=2

∇2X(i/n)2

n − 1

.

(cid:3)

(4.11)

i)

ii)

= Op(n−.5)

(cid:19)2a1

(cid:18)cos θ

− EPn
EQn
−

sin θ

Pn
Qn
EPn
EQn

= O(n2(a1−a2) sin2(a2−a1)

θ ∨ n2(a1−a2) cos2(a2+a1)

θ/sin4a1 θ)

σ2
P := lim
n

var(n.5+2a1Pn) = lim
n

2(ρa1(k) cos2a1 θ + na1−a2 ρa2(k) sin2a2 θ)2

Proof. i) Using (4.5-4.6), Lemma 4.1.1 iii),

n


k=−n

n


k=−n

n


n


≤ lim
n

cov(∇1X(2)2, ∇1X(i)2) < ∞,

i=2
var(n.5+2a1Qn) = lim
n

σ2
Q := lim
n

≤ lim
n

cov(∇2X(2)2, ∇2X(i)2) < ∞,

i=2
cov(n.5+2a1Pn,(n.5+2a1Qn) < ∞,

σPQ := lim
n

72

2(ρa1(k) sin2a1 θ + na1−a2 ρa2(k) cos2a2 θ)2

(4.12)

(4.13)

(4.14)

(4.15)

(4.16)

where ρaj(k) = cov(∇Baj(i), ∇Baj(i + k)), and ∇Baj(i) = Baj(i)− 2Baj(i − 1) + Baj(i − 2), Baj(i)
is a fractional Brownian motion with Hurst index aj, j = 1, 2. Therefore by Theorem 2 in [5],

n.5+2a1(cid:169)(cid:173)(cid:173)(cid:171) Pn − EPn

Qn − EQn

(cid:170)(cid:174)(cid:174)(cid:172) →d N(0, Σ),

where Σ11 = σ2
are of order n−2a1 by Lemma 4.1.1 ii).

P, Σ22 = σ2

Q, Σ12 = σPQ. i) follows using delta method and the fact that EPn, EQn

(4.17)

(cid:19)

(cid:169)(cid:173)(cid:173)(cid:171) 0

2/n

(cid:170)(cid:174)(cid:174)(cid:172)

.

(4.18)

(cid:3)

− EPn
n.5( Pn
Qn
EQn
+ (EPn)2
where σ2
(EQn)4 σ2
ii) By Lemma 4.1.1 ii), it is easily derived that

= limn n4a1(

(EQn)2 σ2

P/Q

P

1

) →d N(0, σ2

P/Q),

Q − 1
EQn
(cid:19)(cid:30)(cid:18)

(cid:18)

4υE,H

(cid:169)(cid:173)(cid:173)(cid:171)1/n

0

(cid:170)(cid:174)(cid:174)(cid:172) − υE,H

(cid:169)(cid:173)(cid:173)(cid:171)2/n

0

(cid:170)(cid:174)(cid:174)(cid:172)

EPn

(EQn)2 σPQ).

(cid:169)(cid:173)(cid:173)(cid:171) 0

1/n

(cid:170)(cid:174)(cid:174)(cid:172) − υE,H

4υE,H

cos2a1 θ + n2(a1−a2)(4 − 22a2)(4 − 22a1)−1 sin2a2 θ
sin2a1 θ + n2(a1−a2)(4 − 22a2)(4 − 22a1)−1 cos2a2 θ

EPn
EQn

=

=

Define for i = 4, ..., n,

1X(i/n) = X(cid:169)(cid:173)(cid:173)(cid:171)i/n

0

∇∗

Let

(cid:170)(cid:174)(cid:174)(cid:172) + X(cid:169)(cid:173)(cid:173)(cid:171)(i − 4)/n

0

(cid:170)(cid:174)(cid:174)(cid:172) .

(cid:170)(cid:174)(cid:174)(cid:172) − 2X(cid:169)(cid:173)(cid:173)(cid:171)(i − 2)/n
n


∇∗
1X(i/n)2
n − 3

0

i=4

.

P∗
n :=

73

Corollary 4.1.3. Under the assumptions (A1, A2),
i) P∗

− υE,H

υE,H

n
Qn

0

υE,H

(cid:30)
(cid:169)(cid:173)(cid:173)(cid:171)4/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:30)
(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:30)
(cid:169)(cid:173)(cid:173)(cid:171)4/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:30)
(cid:169)(cid:173)(cid:173)(cid:171)2/n

2/n

υE,H

0

0

υE,H

sin θ

2/n

= O(n2(a1−a2) sin2(a2−a1)

(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:170)(cid:174)(cid:174)(cid:172) = Op(n−.5)
(cid:18)2 cos θ
(cid:19)2a1
(cid:170)(cid:174)(cid:174)(cid:172) −
(cid:169)(cid:173)(cid:173)(cid:171)2/n
(cid:170)(cid:174)(cid:174)(cid:172) = Op(n−.5)
(cid:170)(cid:174)(cid:174)(cid:172) − 22a1 = O(n2(a1−a2) sin2a2 θ/cos2a1 θ)

0

0

(cid:169)(cid:173)(cid:173)(cid:171)4/n
(cid:169)(cid:173)(cid:173)(cid:171)4/n

0

(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172)

− υE,H

ii) υE,H

iii) P∗

n
Pn

iv) υE,H

θ ∨ n2(a1−a2) cos2a2 θ/sin2a1 θ)

The estimation method for θ is the following:
Qn and Pn
n
Qn .

Step 1) Estimate 2a1 by the ratios of P∗

(cid:19)

/log 2.

(cid:18) P∗

n
Pn

ˆ2a1 = log

Step 2) Estimate θ with ˆ2a1 and Pn
Qn .

cos ˆθ
sin ˆθ

=

Then the estimator of θ is

(cid:18) Pn
(cid:19)1/ ˆ2a1
(cid:18)(cid:18) Pn
(cid:19)1/ ˆ2a1(cid:19)

Qn

.

Qn

.

(4.19)

(4.20)

(4.21)

ˆθn := cot−1
Theorem 4.1.4. Under the assumptions (A1, A2)
n2(a1−a2)(sin2a2 θ
sin2a1 θ

(cid:18)

i)
ii)

ˆ2a1 − 2a1 = Op
ˆθn − θ = Op(n2(a1−a2)(sin2a2 θ ∨ sin2−4a1 θ cos2(a2+a1)

) ∨ n−.5

(cid:19)
(cid:1)2a1 and P∗

n
Qn

∨ cos2(a2+a1) θ
sin4a1 θ
−(cid:0) cos θ

sin θ

Proof. By (4.19), ˆ2a1 − 2a1 has the same order as Pn
Qn
i) follows from Lemma 4.1.2 and its corollary.

θ) ∨ n−.5 sin2 θ)

−(cid:0) 2 cos θ

sin θ

(4.22)

(cid:1)2a1, therefore

74

For ii), note that(cid:0) Pn

Qn

and

(cid:1)1/ ˆ2a1 − cos θ

sin θ has the same order of ˆ2a1 − 2a1 and Pn
Qn

−(cid:0) cos θ

sin θ

(cid:1)2a1. By (4.20)

dθ

d cot θ

= − sin2 θ,

the result follows since

ˆθn − θ = Op(cid:0)dθ/d cot θ(cid:0) Pn

(cid:1)(cid:1)
1/ ˆ2a1 − cos θ
sin θ
∨ cos2(a2+a1) θ
= Op(n2(a1−a2) sin2 θ(sin2a2 θ
sin4a1 θ
sin2a1 θ
= Op(n2(a1−a2)(sin2a2 θ ∨ sin2−4a1 θ cos2(a2+a1)

Qn

) ∨ n−.5 sin2 θ)
θ) ∨ n−.5 sin2 θ).

Remark 2. (a). Theorem 4.1.4 implies that the estimates of 2a1 and cos θ/sin θ , ˆ2a1 and(cid:0) Pn

Qn

(cid:3)

(cid:1)1/ ˆ2a1,

n and Pn directly, instead of the ratios of P∗

respectively, are not good when θ is close to zero. But the estimate of θ is uniformly good throughout
θ ∈ [0, π).
(b). In (4.19) (step 1), one can use the ratio of P∗
n/Qn
and Pn/Qn to estimate 2a1, and will get the same asymptotic result as in Theorem 4.1.4, which is
because of the previous corollary iii),iv).
(c). Theorem 4.1.4 implies consistency of 2 ˆa1, ˆθn in probability. Also note that the bias and variance
n2(a1−a2)(sin2a2 θ
, Op(n−1), respectively, following the
of 2 ˆa1 are of order Op
sin2a1 θ
P∗
Qn . Likewise, the bias and the variance of ˆθn are of order
order of the bias and variance of Pn
n
Qn,
n2(a1−a2)(sin2a2 θ ∨ sin2−4a1 θ cos2(a2+a1) θ), n−1 sin4 θ, respectively.
(d). By (4.12-4.15), it can be seen that the variances σ2
Q, σP/Q in (4.16-4.18) are dependent on
a1, thus the variance of 2 ˆa1, ˆθn are dependent on a1. More specifically, the smaller a1 is, the bigger
the variances σ2

Q, σP/Q are, as well as the variances of 2 ˆa1, ˆθn.

∨ cos2(a2+a1) θ
sin4a1 θ

P, σ2

P, σ2

(cid:18)

(cid:19)

)

75

4.1.2 Estimation of ratios of H(= 1) to eigenvalues of E

Define ˆu1 =(cid:169)(cid:173)(cid:173)(cid:171) cos ˆθ

− sin ˆθ

(cid:170)(cid:174)(cid:174)(cid:172) , and ˆu2 =(cid:169)(cid:173)(cid:173)(cid:171)sin ˆθ

cos ˆθ

(cid:170)(cid:174)(cid:174)(cid:172) . Using the relation,

{X(cuit)}t∈R

f .d.

= {c1/λi X(uit)}t∈R,

(4.23)
we estimate ai = 1/λi, i = 1, 2, which are Hurst indices along the directions ui, i = 1, 2, respectively.
Among the many well-known Hurst estimation methods, here, we employ discrete variation method.

(cid:18)i − 2m+1

(cid:19)

For a fixed domain, [0, 1] × [0, 1], and for m ∈ N ∩ {0}, 2m << n, define
ˆ∇m
j X(i/n) = X

j = 1, 2, i = 2m+1, ..., n.
j X(i/n)}i is a discretized sample path with direction ˆuj . If ˆθ = 0 then ˆ∇0

− 2X

+ X

ˆuj

ˆuj

ˆuj

n

n

n

(cid:18)i − 2m

(cid:19)

(cid:18) i

(cid:19)

(4.24)
j X(i/n) =

Note that { ˆ∇m
∇j X(i/n), otherwise ˆ∇0

ˆPm
n :=

j X(i/n) (cid:28) ∇j X(i/n). Similar to Pn, Qn in (4.11), define
n

ˆ∇m
2 X(i/n)2
n − 2m+1 + 1 .
n ; m = 1, 2, ..., (cid:96)n}, { ˆQm

ˆ∇m
1 X(i/n)2
n − 2m+1 + 1,

n


ˆQm
n :=

i=2m+1

i=2m+1

n ; m = 1, 2, ..., (cid:96)n} on

The estimates of a1 and a2 are log-regression of { ˆPm
{2m log 2, m = 1, 2, ..., (cid:96)n} respectively, i.e.
1
2

ˆa1 =

wm log2 ˆPm
n ,

(4.25)

(4.26)

(4.27)

(cid:96)n

(cid:96)n


m=1

m=1

ˆa2 =

1
2

(cid:96)n


m=1

(cid:96)n

(cid:96)n


m=1

m=1

aE
1 :=

aE
2 :=

1
2

1
2

wm log2 ˆQm
n ,

(cid:96)n


m=1

wm = 0,

mwm = 1.

n ),
wm log2 EX| ˆu1( ˆPm

wm log2 EX| ˆu2( ˆQm
n ),

where

Define

where the expectation is for X given ˆuj, j = 1, 2.

76

Theorem 4.1.5. In a fixed domain with the assumptions (A1, A2),
i)



1 − a1 =
aE

ii) If a2 − a1 ≤ .25 and a1 > .5,

n.5+2a1( ˆa1 − aE

OP(n−2(a2−a1)−a2)
OP(n−(2+4a2)(a2−a1))

1 ) →d N(0, σa1),
if a2 − a1 > .25,
if a2 − a1 ≤ .25.

Proof. For aE

1 , by (4.2),
EX| ˆu1( ˆ∇m

2 ) →d N(0, σa2),
n.5+2a2( ˆa2 − aE
2 − a2 = OP(n−(4a1−2)(a2−a1)).
aE
(cid:19)

(cid:18)22ma2

2a2
+ ca2

1
(1 + ca2/c∗
a1


n2a2
2a2
1

(cid:18)22ma1
(cid:18)22ma1

n2a1

(cid:19)
(cid:19)

n2a1

1 X(i/n)2) = c∗
a1
= c∗
a1

22m(a2−a1)n2a1−2a2),

= E(∇Ba1(i)2)(u(cid:48)

1 ˆu1)2a1, ca2 = E(∇Ba2(i)2), and 
1 = u(cid:48)

where c∗
a1
(4.22). Note that u(cid:48)
n−(2+4a2)(a2−a1) depending on whether a2 − a1 > .25 or a2 − a1 ≤ .25, the result follows.
For aE
2 ,

2 ˆu1 which has the order of
n2a1−2a2 is either less than n2(a1−a2)−a2 or

1 ˆu1 ∼ 1. Since the order of 


2a2
1

EX| ˆu2( ˆ∇m

(cid:18)22ma1
(cid:18)22ma2
(cid:19)
2 X(i/n)2) = ca1

= c∗
a2
= E(∇Ba2(i)2)(u(cid:48)

(cid:19)
+ c∗
a2
n2a1
(1 + ca1/c∗
a2


2a1
2

n2a2

(cid:18)22ma2

(cid:19)

n2a2
2a1
2

22m(a1−a2)n2a2−2a1),

(4.28)

(4.29)

(4.30)

(4.31)

where ca1 = E(∇Ba1(i)2), c∗
1 ˆu2 which has the order of
a2
n2a2−2a1 is of order n−(4a1−2)(a2−a1) when a2 − a1 ≤ .25, and a1 > .5 or
2a1
(4.22). Since 

2
divergent otherwise, the result follows. Asymptotic normality of ˆa1, ˆa2 follows from (4.17) with
(cid:3)
delta method.

2 ˆu2)2a2 and 
2 = u(cid:48)

77

Remark 3. (a). If the condition of ii) in Theorem 4.1.5 is not met, then ˆa2 goes to a1. This is
2a1
because 

2

log2(1 + ca1/c∗
2a1
a2

2

n2a2−2a1 in (4.31) diverges and therefore as n → ∞
22m(a1−a2)n2a2−2a1) (cid:39) log2(ca1/c∗
a2

(cid:96)n


and

22m(a1−a2)n2a2−2a1),

2a1
2

aE
2 =

1
2

m=1

wm log2 EX| ˆu2( ˆ∇m

2 X(i/n)2) ∼ a2 + (a1 − a2) = a1.

(b). Note that here it is assumed that independent sample paths are observed. However, it is
i − ai, i = 1, 2, even when ˆui and {X(t ˆui)}t∈R are
expected to have the same bias results for aE
dependent (e.g., they are in the same sample surface.) This is because ˆPm
n converge a.s.
to its expectation. For the same reason, the bias of the estimators in Theorem 4.1.6 and Theorem

n and ˆQm

4.1.10 remain true when {X(cid:169)(cid:173)(cid:173)(cid:171)t

0

(cid:170)(cid:174)(cid:174)(cid:172) , X(cid:169)(cid:173)(cid:173)(cid:171)0

t

(cid:170)(cid:174)(cid:174)(cid:172) , X(t ˆui)}t∈R is from one sample surface.

(c). Note that a1 is estimated two times in the whole estimating procedure: The first time is when θ
is estimated in section 4.1.1 (Theorem 4.1.4) and the other is in this section 4.1.2 (Theorem 4.1.5).
This is because ˆ2a1 was estimated to get the estimate of θ so that (a1, a2) are estimated. However,
ˆ2a1 estimated the first time has both larger bias and variance (Theorem 4.1.4) compared to ˆa1 in
this section (Theorem 4.1.5).

As the above lemma indicates that the estimator ˆa2 is not good when a1 is small, one can naturally
use the sample in a different domain (i.e., an increasing domain.) In an increasing domain, the
j X(i).
estimators are obtained in the same way in (4.25-4.27) except that ˆ∇m

j X(i/n) is changed to ˆ∇m

Theorem 4.1.6. For an increasing domain with the assumptions (A1, A2),
i)



1 − a1 =
aE

1 ) →d N(0, σa1),
if a2 − a1 > .25,
if a2 − a1 ≤ .25.

√

n( ˆa1 − aE

OP(n−a2)
OP(n−4a2(a2−a1))

78

ii)

Proof.



2 − a2 =
aE

√

n( ˆa2 − aE

OP(n−a1)
OP(n−4a1(a2−a1))

2 ) →d N(0, σa2),
if a2 − a1 > .25,
if a2 − a1 ≤ .25.

EX| ˆu1( ˆ∇m

1 X(i)2) = c∗
= c∗

(22ma2)
2a2
a122ma1 + ca2

1
a122ma1(1 + ca2/c∗
2a2
a1

1

22m(a2−a1)),

= E(∇Ba1(i)2)(u(cid:48)

1 ˆu1)2a1, ca2 = E(∇Ba2(i)2) and 
1 = u(cid:48)

where c∗
2 ˆu1 which has the order of
a1
(4.22). Since the order of 
2a2 is either less than n−a2 or n−4a2(a2−a1) depending on whether
a2 − a1 > .25 or a2 − a1 ≤ .25, the result follows. For ii),
22ma1 + c∗
a222ma2
a222ma2(1 + ca1/c∗
2a1
a2

2

2 X(i)2) = ca1

= c∗

22m(a1−a2)),

EX| ˆu2( ˆ∇m

2a1
2

= E(∇Ba2(i)2)(u(cid:48)

2 ˆu2)2a2 and 
2 = u(cid:48)

1 ˆu2 which has the order of
2 − a2. Asymptotic normality of ˆa1, ˆa2 follows from

where ca1 = E(∇Ba1(i)2), c∗
a2
(4.22). The result follows for the order of aE
(cid:18)
n

the fact that
(cid:18)
n


ˆ∇m
1 X(i/n)2
n − 2m+1 + 1
ˆ∇m
2 X(i/n)2
n − 2m+1 + 1

i=2m+1

√

√

n

n

i=2m+1

(cid:19)
(cid:19)

− E( ˆ∇m

1 X(i)2)

− E( ˆ∇m

2 X(i)2)

→d N(0, Σ1),

→d N(0, Σ2),

which can be derived in the same way as (4.17).

(cid:3)

4.1.3

For H < 1

From now on, the results for H = 1 has been shown. This constraint is relaxed in this section, and it
will be seen that the results remain the same except the estimators are for hi = H/λi = Hai, i = 1, 2.

79

Since H < 1, (4.5-4.6),(4.9-4.10) no longer hold. Nonetheless, most of the results hold in the
similar way.

ii) E(∇1X(i/n)2) = 8υE,H

ASSUMPTION
(A1(cid:48)) υE,H(x) = τE(x)2H with τE in (2), H < 1 and a1 ≤ a2, a1, a2 ∈ (0, 1).
Lemma 4.1.7. Under the assumptions (A1(cid:48), A2),

0

i) υE,H

(cid:169)(cid:173)(cid:173)(cid:171)h
iii) limn n4h1
n

Proof. i) Since

0

0

h

i=2 cov(∇j X(2/n)2, ∇j X(i/n)2) < ∞,

(cid:169)(cid:173)(cid:173)(cid:171)0
(cid:170)(cid:174)(cid:174)(cid:172) ∼ |h sin θ|2h1 for |h| → 0.
(cid:170)(cid:174)(cid:174)(cid:172) ∼ |h cos θ|2h1 and υE,H
(cid:169)(cid:173)(cid:173)(cid:171)1/n
(cid:169)(cid:173)(cid:173)(cid:171)2/n
(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H
(cid:170)(cid:174)(cid:174)(cid:172) , E(∇2X(i/n)2) = 8υE,H
(cid:169)(cid:173)(cid:173)(cid:171)h
(cid:170)(cid:174)(cid:174)(cid:172) =(|h cos θ|2a1 + |h sin θ|2a2)H,
(cid:169)(cid:173)(cid:173)(cid:171)0
(cid:170)(cid:174)(cid:174)(cid:172) =(|h sin θ|2a1 + |h cos θ|2a2)H,

j = 1, 2.

υE,H

υE,H

0

(cid:169)(cid:173)(cid:173)(cid:171) 0

1/n

(cid:170)(cid:174)(cid:174)(cid:172) − 2υE,H

(cid:169)(cid:173)(cid:173)(cid:171) 0

2/n

(cid:170)(cid:174)(cid:174)(cid:172) .

h

the results are derived. ii) It is proved in the same way as in Lemma 4.4.1.
iii) Let j = 1. For j = 2, the proof goes in the same way. Note that cov(∇j X(2/n)2, ∇j X(i/n)2) =
2cov(∇j X(2/n), ∇j X(i/n))2. For large enough n,

cov(∇1X(2/n), ∇1X(i/n)) ∼ c1υ

(4)
E,H

(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172) is the fourth derivative of υE,H

0

for some constant c1, where υ

(cid:169)(cid:173)(cid:173)(cid:171)i/n

0

(4)
E,H

(cid:169)(cid:173)(cid:173)(cid:171)i/n

0

(cid:170)(cid:174)(cid:174)(cid:172) with respect to i. This

80

cov(∇1X(2/n), ∇1X(i/n)) = −2

0

0

0

0

0

2

2 + E

(cid:170)(cid:174)(cid:174)(cid:172) X(cid:169)(cid:173)(cid:173)(cid:171)1
E(X(cid:169)(cid:173)(cid:173)(cid:171)i + 1
(cid:170)(cid:174)(cid:174)(cid:172)) = −2υE,H
(cid:26)

2(cid:19)(cid:30)
2(cid:19)(cid:30)
(cid:18)
(cid:18)
X(cid:169)(cid:173)(cid:173)(cid:171)1
X(cid:169)(cid:173)(cid:173)(cid:171)i + 1
(cid:169)(cid:173)(cid:173)(cid:171)i
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172) + E
(cid:169)(cid:173)(cid:173)(cid:171)(i − 1)/n
(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:169)(cid:173)(cid:173)(cid:171)(i + 1)/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172) + 6υE,H
(cid:170)(cid:174)(cid:174)(cid:172) − 4υE,H
(cid:27)
(cid:169)(cid:173)(cid:173)(cid:171)(i − 3)/n
(cid:169)(cid:173)(cid:173)(cid:171)(i − 2)/n
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:170)(cid:174)(cid:174)(cid:172) + υE,H
(cid:170)(cid:174)(cid:174)(cid:172) is the same as the fourth derivative of (|i cos θ|2a1 + n2a1−2a2|i sin θ|2a2)H,

(cid:169)(cid:173)(cid:173)(cid:171)i/n

0

− 4υE,H

υE,H

0

0

0

0

0

.

is because

and therefore

Since n2h1υ

(4)
E,H

which is


k≥(cid:96)≥m≥p≥0
k+(cid:96)+m+p=4

ck,(cid:96),m,p(|i cos θ|2a1 + n2(a1−a2)|i sin θ|2a2)H−k(|i cos θ|2a1−1 + n2(a1−a2)|i sin θ|2a2−1)k−(cid:96)

× (|i cos θ|2a1−2 + n2(a1−a2)|i sin θ|2a2−2)(cid:96)−m(|i cos θ|2a1−3 + n2(a1−a2)|i sin θ|2a2−3)m−p
× (|i cos θ|2a1−4 + n2(a1−a2)|i sin θ|2a2−4)p < c2i2a1H−4

for some constants ck,(cid:96),m,p, c2. Therefore n2h1cov(∇1X(2/n), ∇1X(i/n)) < c2i2a1H−4, and the
results follow since a1H < 1.

(cid:3)

Lemma 4.1.8. Assumptions (A1(cid:48), A2) are satisfied. Then,

(cid:170)(cid:174)(cid:174)(cid:172)

(cid:169)(cid:173)(cid:173)(cid:171)2/n
(cid:30)

0

υE,H

υE,H

(cid:30)
(cid:169)(cid:173)(cid:173)(cid:171) 0

2/n

2/n

(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:170)(cid:174)(cid:174)(cid:172) = Op(n−.5)
(cid:19)2h1
(cid:18)cos θ
(cid:170)(cid:174)(cid:174)(cid:172) −

sin θ

i) Pn
Qn

ii) υE,H

− υE,H

(cid:169)(cid:173)(cid:173)(cid:171)2/n

0

(cid:170)(cid:174)(cid:174)(cid:172)

= O(n2(a1−a2) sin2(a2−a1)

θ ∨ n2(a1−a2) cos2a2 θ/sin2a1 θ)

81

follows

Therefore

n.5+2h1(Pn − EPn) → N(0, σP),
n.5+2h1(Qn − EQn) → N(0, σQ).

− EPn
EQn
Q).
σ2
where σ2
For ii), it is easily derived from the fact that

= limn n4h1( 1
EQ2
n

n.5( Pn
Qn
+ EP2
n
EQ4
n

P/Q

σ2
P

) → N(0, σP/Q),

(cid:30)

(cid:169)(cid:173)(cid:173)(cid:171)2/n

0

(cid:170)(cid:174)(cid:174)(cid:172)

υE,H

υE,H

(cid:18)cos2a1 θ + 2n2a1−a2 sin2a2 θ

sin2a1 θ + 2n2a1−a2 cos2a2 θ

(4.32)

(4.33)

(cid:3)

(cid:19) H

.

Proof. i) By Theorem 2 in [5] and

σ2
P := lim
n

σ2
Q := lim
n

var(n.5+2h1Pn) < 2 lim
n
var(n.5+2h1Qn) < 2 lim
n

n4h1

n4h1

n

n


i=2

i=2

cov(∇1X(2/n)2, ∇1X(i/n)2) < ∞,

cov(∇1X(2/n)2, ∇1X(i/n)2) < ∞,

By the above results and the same way as in Theorem 4.1.4, the next theorem follows.

Theorem 4.1.9. Under the assumptions (A1(cid:48), A2),

i)
ii)

ˆ2h1 − 2h1 = Op
ˆθn − θ = Op(n2(a1−a2) ∨ n−.5)
(4.23) is now changed to the following equation.

n2(a1−a2)(sin2a2 θ
sin2a1 θ

(cid:19)

∨ cos2a2 θ
sin2a1 θ

) ∨ n−.5

{X(cuit)}t∈R

f .d.

= {cH/λi X(uit)}t∈R.

Using the same estimator as (4.26-4.27), now one can obtain the estimator of hi, hi = H/λi, i = 1, 2.

wm log2 ˆPm
n .

wm log2 ˆQm
n ,

(4.34)

(4.35)

(cid:170)(cid:174)(cid:174)(cid:172) =

2/n

(cid:169)(cid:173)(cid:173)(cid:171) 0
(cid:18)

(cid:96)n

(cid:96)n


m=1

m=1

ˆh1 =

ˆh2 =

1
2

1
2

82

n , ˆQm

n defined same as (4.25).

where ˆPm
Theorem 4.1.10. In a fixed domain with the assumptions (A1(cid:48), A2),
i)



1 − h1 =
hE

ii) If a2 − a1 ≤ .25 and a1 > .5,

n.5+2a1(ˆh1 − hE

OP(Hn−2(a2−a1)−a2)
OP(Hn−(2+4a2)(a2−a1))

1 ) →d N(0, σh1),
if a2 − a1 > .25,
if a2 − a1 ≤ .25.

Proof. For hE

1 , by Lemma 4.1.7 ii), it is derived that

EX| ˆu1( ˆ∇m

1 X(i/n)2) = 8

ba1

+ 


n2a1

(cid:19)

2a2
1

(cid:18)22ma2

(cid:19)(cid:19) H − 2

2 ) → N(0, σh2),
n.5+2a2(ˆh2 − hE
2 − h2 = OP(Hn−(4a1−2)(a2−a1)).
hE
(cid:19)
(cid:18)22(m+1)a1
(cid:18)
(cid:18)
(cid:18)22ma1
(cid:18)22mHa1
(cid:19)(cid:16)
ba1
n2a2
n2a1
a1 n2a1−2a2)H
22m(a2−a1)b−1
a1(1 + 

2a2
8bH
1
(cid:18)22mh1
(cid:19)(cid:16)
a122Ha1(1 + 22(a2−a1)
− 2bH
1 + ˜ca2/c∗
h1
n2h1
1 ˆu1)2h1, ˜ca2 = E(∇Ba2−(1−H)a1(i)2)bH−1
= E(∇Bh1(i)2)(u(cid:48)
a1

22m(a2−a1)b−1
22m(a2−a1)n2a1−2a2 + o(

2a2
1

2a2
H

1

n2Ha1

2a2
1

2a2
1




+ 


n2a2

(cid:19)(cid:19) H
(cid:18)22(m+1)a2
a1 n2a1−2a2)H(cid:17)
n2a1−2a2)(cid:17)

,

=

= c∗
h1
1 ˆu1)2a1, c∗
h1

where ba1 = (u(cid:48)
and 
1 =
u(cid:48)
2 ˆu1. The last equality follows by Taylor expansion and the fact that E(∇Bh(i)2) = 8 − 2 · 22h.
n2a1−2a2 is either less than n2(a1−a2)−a2 or n−(2+4a2)(a2−a1)
Since 
1 has the order of (4.22), 

depending on whether a2 − a1 > .25 or a2 − a1 ≤ .25, the result follows.

2a2
1

83

For hE

2 , in the same way as above, it is derived that

EX| ˆu2( ˆ∇m

2 X(i/n)2) = 8

n2a1




(cid:19)

2a1


2

(cid:18)22ma2

(cid:19)(cid:19) H − 2

(cid:18)
(cid:18)22(m+1)a1
n2a1
a2 n2a2−2a1)H

(cid:19)
(cid:18)
(cid:18)22ma1
(cid:19)(cid:16)
(cid:18)22mHa2
2a1
+ ba2
2
n2a2
22m(a1−a2)b−1
a2(1 + 

2a1
8bH
2
(cid:18)22mh2
(cid:19)(cid:16)
a222Ha2(1 + 22(a1−a2)
− 2bH
1 + ˜ca1/c∗
h2
n2h2
2 ˆu2)2h2, ˜ca1 = E(∇Ba1−(1−H)a2(i)2)bH−1
= E(∇Bh2(i)2)(u(cid:48)
a2

22m(a1−a2)b−1
22m(a1−a2)n2a2−2a1 + o(

2a1
2

n2a2

(cid:19)(cid:19) H
(cid:18)22(m+1)a2
a2 n2a2−2a1)H(cid:17)
n2a2−2a1)(cid:17)

2a2
H

2

2a2
2




+ ba2

=

n2Ha2

= c∗
h2
2 ˆu2)2a2, c∗
h2

,

where ba2 = (u(cid:48)
and 
2 =
u(cid:48)
1 ˆu2. The last equality follows by Taylor expansion and the fact that E(∇Bh(i)2) = 8 − 2 · 22h.
n2a2−2a1 is of order n−(4a1−2)(a2−a1) when a2− a1 ≤ .25, and
Since 
2 has the order of (4.22), 

2 − h2. Asymptotic normality
a1 > .5 or divergent otherwise, the results follow for the order of hE
of ˆh1, ˆh2 follows from (4.32-4.33).
(cid:3)

2a1
2

4.2 Simulation results

The simulation of operator scaling Gaussian random field with semi-variogram (1.8) and
diagonal matrix E was performed with the algorithm in [9](Figure 4.1, Figure 4.6). As we have
diagonizable matrix E in this chapter, the algorithm in [9] was modified. Two different sampling
regimes were employed in this section: 1) independent sample paths on exact directions 2) sample
surface on grid lines. In Section 4.2.1, independent sample paths were simulated on exact directions.
In Section 4.2.2, sample surfaces were simulated on grid lines.

4.2.1

Independent sample paths on exact directions

Here it is assumed that independent sample paths on exact directions are obtained from OSGRF
with different parameters (θ, a1, a2) when H = 1, (θ, h1, h2) when H < 1. More specifically, two

independent sample paths, {X(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:170)(cid:174)(cid:174)(cid:172) , n = 213, i = 1, .., 213}, are
n ˆu1(cid:1), n =
obtained to estimate θ, and with that directions ˆθ, two independent sample paths, {X(cid:0) i

(cid:170)(cid:174)(cid:174)(cid:172) , n = 213, i = 1, ..., 213}, {X(cid:169)(cid:173)(cid:173)(cid:171) 0

i/n

0

84

213, i = 1, ..., 213}, {X(cid:0) i

n ˆu2(cid:1), n = 213, i = 1, ..., 213}, are used to estimate (a1, a2) for H = 1 and

(h1, h2) for H < 1. In an increasing domain regime, {X(i ˆu1), i = 1, ..., 213}, {X(i ˆu2), i = 1, ..., 213}
are used to estimate (a1, a2) for H = 1.

4.2.1.1

For H = 1

Figure 4.1 shows the simulations of random field X whose semi-variogram is (1.8) with θ = 0 and
various a1, a2 in a fixed domain [0, 1]×[0, 1]. In the same model with different θ, the picture would
look just like Figure 4.1 in a different angle (i.e., the picture rotated by the angle θ from Figure 4.1.)
In Figure 4.2, cos θ/sin θ was estimated in OSGRF for {θ = i/20 ∗ π/2, i = 1, ..., 19}, for each

of {(a1, a2)} in Figure 4.1. In other words, (cid:18) Pn

(cid:19)1/ ˆαn

Qn

are drawn and its convergence rate is in (4.21). It is seen that, for θ close to zero, the estimation is
not good as it was expected by (4.21), but as seen in Figure 4.3, the estimate of θ, which one needs
ultimately, is good for the whole range of θ, and this result was expected from (4.22).

In Figure 4.3, the estimates of θ, θ ∈ [0, π/2], were drawn in OSGRF for {θ = i/20 ∗ π/2, i =
1, ..., 19} in each case of {(a1, a2)} as in Figure 4.1. For these estimates, the sample paths of

{X(cid:169)(cid:173)(cid:173)(cid:171)i/n

0

(cid:170)(cid:174)(cid:174)(cid:172) ; n = 213, i = 1, .., 213}, {X(cid:169)(cid:173)(cid:173)(cid:171) 0

(cid:170)(cid:174)(cid:174)(cid:172) , n = 213, i = 1, 2, .., 213} were used.

It is seen from
Figure 4.3 that ˆθ works well in any range of set (θ, a1, a2) as all other lines are close to the black
line.

i/n

In Figure 4.4, the estimates of a1, a2 were drawn for each case of (a1, a2) with varying θ. The
different colors represent a different set of (a1, a2), and in each color, the solid line is for ˆa1, and the
dotted line is for ˆa2. The estimates ˆa1, ˆa2 were obtained in a fixed domain [0, 1] × [0, 1] with m =

n ˆu1(cid:1), n = 213, i = 1, .., 213}, {X(cid:0) i

n ˆu2(cid:1), n = 213, i = 1, .., 213}

1, 2, 3, 4(= (cid:96)n), and sample paths {X(cid:0) i

for a1, a2, respectively, where ˆθ for { ˆuj( ˆθ), j = 1, 2} are from the previous step in Figure 4.3. As it
was expected in Theorem 4.1.5, Figure 4.4 reveals the fact that ˆa1 performs well for any range of
the set (θ, a1, a2), but ˆa2 performs well only when a1 is large and a2 − a1 is small– {(a1, a2); a1 >

85

.5 and a2 − a1 < .25}. In Figure 4.4, it is seen that for {(a1, a2) = (.1, .3),(.2, .3),(.3, .8)}, ˆa1 and
ˆa2 are both close to a1. For {(a1, a2) = (.6, .8),(.8, .9)}, ˆa1 is close to a1 and ˆa2 is close to a2.

In Figure 4.5, the same estimates of θ were used as in Figure 4.3 and Figure 4.4, but a1, a2
are estimated in an increasing domain [0, 213] × [0, 213]. More specifically, the estimates ˆa1, ˆa2
were obtained with m = 1, 2, 3, 4(= (cid:96)n), and the sample paths {X(i ˆu1), i = 1, .., 213}, {X(i ˆu2), i =
1, .., 213} for a1, a2, respectively, where ˆθ for { ˆuj( ˆθ), j = 1, 2} are the same as in Figure 4.3, and
Figure 4.4. As it was expected in Theorem 4.1.6, it is seen in Figure 4.5 that both ˆa1, ˆa2 work well
in any set of the range (θ, a1, a2), as the solid line which is for ˆa1 is close to a1 and the dotted line
which is for ˆa2 is close to a2 in any case of (θ, a1, a2).

For H < 1

4.2.1.2
Figure 4.6 shows sample surface of OSGRF with various (a1, a2) when H < 1, and θ = 0. The
sample surfaces look different than those in Figure 4.1 as there is no longer clear grid lines in
Figure 4.6. Figure 4.7 shows the results for the estimator ˆθ, and Figure 4.8 shows the results for
the estimator ˆh1, ˆh2 with m = 1, 2, 3, 4(= (cid:96)n). The estimators ˆθ, ˆh1 ˆh2 behave similarly as in the case
when H = 1. In Figure 4.7, ˆθ was drawn for {θ = i/20 ∗ π/2, i = 1, ..., 19}, for each (h1, h2) with
various H in each of the graphs (a)–(f). It is noticeable that the bias of ˆθ increases as h2 − h1 and
therefore a2−a1 decreases, and the standard error of ˆθ increases as a1, a2 are getting smaller. These
results were predicted in Theorem 4.1.9. In Figure 4.8, it is observed that for graphs (a), (b), (c),
the estimators ˆh1, ˆh2 work well for any range of θ, whereas for graphs (d), (e), only ˆh1 works well.
For graph (f), ˆh1 works well, and ˆh2 performs better than that of graphs (d), (e) but not as much
as graph (a-c). These results coincide with Theorem 4.1.10, as for (a-c), (a1 = h1/H, a2 = h2/H)
falls into the region {(a1, a2); a1 > .5 and a2 − a1 < .25}, whereas in (d), (e), (a1, a2) is far from
the region, and in F, (a1, a2) does not fall into the region but pretty close to the boundary.

86

(cid:27)
(cid:27)
(cid:170)(cid:174)(cid:174)(cid:172) , X(cid:169)(cid:173)(cid:173)(cid:171)(cid:98)(i − n/2) tan θ(cid:99)/n + 1/2
(cid:170)(cid:174)(cid:174)(cid:172) , i = 1, ..., n
(cid:170)(cid:174)(cid:174)(cid:172) , X(cid:169)(cid:173)(cid:173)(cid:171)(cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2
(cid:170)(cid:174)(cid:174)(cid:172) , i = 1, ..., n

i/n

i/n

(cid:27)

(cid:98)(i − n/2) tan θ(cid:99)/n + 1/2

were chosen to estimate θ, and with ˆθ,

S2 =

i/n

(cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2

S1 =

(cid:26)

X(cid:169)(cid:173)(cid:173)(cid:171)

4.2.2

sample surface on a grid

(cid:26)

X(cid:169)(cid:173)(cid:173)(cid:171)i/n

j/n

Here it is assume that sample surface from OSGRF is available on a grid. More specifically,
are simulated and used to estimate (θ, a1, a2) for H = 1 and (θ, h1, h2) for
H < 1. Since OSGRF with the parameter (θ (cid:44) 0, a1, a2, H) is the same as OSGRF (0, a1, a2, H)
were generated from OSGRF (0, a1, a2, H) and the

rotated by the angle θ,

(cid:170)(cid:174)(cid:174)(cid:172) , i, j = 1, ..., n

(cid:27)
(cid:170)(cid:174)(cid:174)(cid:172) , i, j = 1, ..., n
(cid:26)
X(cid:169)(cid:173)(cid:173)(cid:171)i/n
(cid:26)

j/n

i/n

X(cid:169)(cid:173)(cid:173)(cid:171)

sample points

were chosen to estimate (a1, a2) for H = 1, (h1, h2) for H < 1. The estimation method for θ is the
same as before except that ∇j X(i/n), j = 1, 2, in (4.3-4.4) are computed with sample points in S1.
Likewise, the estimation method for aj, hj, j = 1, 2, are the same as in the previous section except
that ˆ∇m

j X(i/n), j = 1, 2, in (4.24) are computed with every 2m-th points in S2.

4.2.2.1 For H = 1
For each θ from {θ = i/20∗π/2, i = 1, ..., 19}, the set S1 with n = 214 was chosen from the simulated
sample surface to estimate θ. The results are shown in Figure 4.9. With ˆθ, and n = 214, the set
j X(i/n) similarly to (4.24) with every 2m-th sample point from S2 to
S2 was chosen to compute ˆ∇m
estimate (a1, a2). Figure 4.10 and Figure 4.11 show the results for ˆa1, ˆa2 with m = 1, 2, 3, 4(= (cid:96)n),
and m = 5, 6, 7, 8(= (cid:96)n), respectively.

Figure 4.10 looks very different than the one obtained in Section 4.2.1.1, especially for the
dotted lines which are for ˆa2. This is because samples are on grid lines so that for small ˆθ − θ, (cid:98)(i −
n/2) tan( ˆθ − θ)(cid:99)/n + 1/2 will be the same for many adjacent i values, which results in sample points

87

on the same vertical or horizontal lines in S2. This affects the estimator of a2 in a different way for
2 X(i/n) is computed
a different value of a2. This is explained well with (4.30-4.31). For ˆa2, as ˆ∇m
with S2,

ˆ∇m

(i − 2m)/n

(i − 2m+1)/n

2 X(i/n) =X(cid:169)(cid:173)(cid:173)(cid:171)(cid:98)(i − 2m+1 − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2
(cid:170)(cid:174)(cid:174)(cid:172)
− 2X(cid:169)(cid:173)(cid:173)(cid:171)(cid:98)(i − 2m − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2
(cid:170)(cid:174)(cid:174)(cid:172)
+ X(cid:169)(cid:173)(cid:173)(cid:171)(cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2
(cid:170)(cid:174)(cid:174)(cid:172)
(cid:16)(Ba1((cid:98)(i − 2m+1 − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2)
(cid:18)22ma2
(cid:19)

− 2Ba1((cid:98)(i − 2m − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2) + Ba1((cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2))2(cid:17)
(cid:16)(Ba2((i − 2m+1)/n) − 2Ba2((i − 2m)/n) + Ba2(i/n))2(cid:17)
(cid:18) 1
(cid:19)

2 X(i/n)2) = E

i = 2m+1, ..., n,

EX| ˆu2( ˆ∇m

i/n

+ E

= ˜c

+ ca2

n2a1

,

n2a2

the expectation of its squared (4.30) is changed to

(cid:26)

where ˜c is either zero or E(Ba1(i) − Ba1(i − 1))2 = 2, and ca2 = E(∇Ba2(i)2). This is because for
ˆθ − θ being small, (cid:98)2m tan( ˆθ − θ)(cid:99) (cid:28) 1, and therefore
(cid:98)(i − 2m+1 − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2, (cid:98)(i −
2m − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2, (cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2
are all the same or at most one of
them has an increment of 1/n for about 2m+1(cid:98)n tan( ˆθ − θ)(cid:99) among {i = 2m+1, ..., n}. This leads the
change in (4.31), replacing ca1/c∗
a2


22m(a1−a2) by ˜c/ca22−2ma2,

2a1
2

(cid:27)

n2a2

2 X(i/n)2) = log2 ca2 + log2

+ log2(1 + ˜c/ca22−2ma2n2a2−2a1).

log2 EX| ˆu2( ˆ∇m

(cid:96)n
Since ˜c is either 0 or 2, the last term is either 0 or close to log2( ˜c/ca22−2ma2n2a2−2a1). Therefore,
m=1 wm log2 E( ˆ∇m
1
2 becomes the weighted average
2
between 0 and a2 with weights around 2m+1 tan( ˆθ − θ), 1 − 2m+1 tan( ˆθ − θ), respectively. Note that
the larger a2 is, the bigger the bias of ˆa2 is, as ˆa2 becomes much smaller than a2. Therefore when

2 X(i/n)2) is either a2 or close to 0, and aE

(cid:18)22ma2

(cid:19)

88

a2 is small, the estimator ˆa2 actually performs better with samples on grid lines than with samples
on the exact direction. However, there is not much difference for ˆa1, as it always performs well
whether samples are on grid or exact directions. This can be explained by (4.28-4.29). Since the
term n2a1−2a2 in (4.29) goes to zero as the sample size grows, it always makes the estimator ˆa2
work well.

If one increases the scale m that is used for the distance between sample points in estimators
ˆa1, ˆa2 (4.26-4.27), then the estimators behave similarly to the previous case where sample points
are laid exactly along the direction. This can be seen in Figure 4.11. Figure 4.11 looks similar
to Figure 4.4 in a way that ˆa2 works well only when parameters are in the region {(a1, a2); a1 >
.5 and a2 − a1 < .25}, and ˆa1 works well in general. It is because by choosing larger scale m, the
effect of sample points lying on the grid lines is diminished. In other words, the distance between
adjacent i, which are used for ˆ∇m
j X(i/n) in (4.24) and in S2, is 2m. Therefore, with the higher m,
(cid:98)(i − n/2) tan( ˆθ − θ)(cid:99)/n + 1/2 would not overlap for these i (i.e. every 2m-th sample in S2 lies on
the different vertical and horizontal grid lines.)

4.2.2.2 For H < 1

(cid:26)

(cid:27)

X(cid:169)(cid:173)(cid:173)(cid:171)i/n

(cid:170)(cid:174)(cid:174)(cid:172) , n = 212, i, j = 1, ..., n

j/n

. The
OSGRF with θ = 0 and various a1, a2, H were simulated at
same method was used to estimate θ, h1, h2 as in section 4.2.2.1. Here, {m = 1, 2, 3, 4(= (cid:96)n)} and
{m = 3, 4, 5, 6(= (cid:96)n)} was used for both h1 = .6, h2 = .7 and h1 = .2, h2 = .4 with various H < 1.
The results are shown in Figure 4.12-4.17. The estimators behave similarly to the case when H = 1
in Section 4.2.2.1. For h1 = .6, h2 = .7, the estimators with larger scale m = 3, 4, 5, 6(= (cid:96)n) perform
well as ˆh1, ˆh2 are close to h1, h2, respectively, whereas the estimator ˆh2 with m = 1, 2, 3, 4(= (cid:96)n) is
poor since both ˆh1, ˆh2 are close to h1. For h1 = .2, h2 = .4, the estimators with m = 1, 2, 3, 4(= (cid:96)n)
performs better as ˆh1, ˆh2 are close to h1, h2 in a wider range of θ than when m = 3, 4, 5, 6(= (cid:96)n).

89

(a)

(c)

(b)

(d)

Figure 4.1: (a): a1 = .1, a2 = .3, (b): a1 = .2, a2 = .3, (c):a1 = .3, a2 = .8, (d):a1 = .6, a2 = .8,
(e): a1 = .8, a2 = .9

(e)

90

Figure 4.2: Estimates of cos θ/sin θ

Figure 4.3: Estimates of θ

91

Figure 4.4: Estimates of a1 and a2 in a fixed domain

Figure 4.5: Estimates of a1 and a2 in an increasing domain

92

(a)

(c)

(b)

(d)

Figure 4.6: (a): h1 = .6, h2 = .7, H = 75, (b): h1 = .6, h2 = .7, H = .85, (c):h1 = .1, h2 = .3, H =
.4, (d):h1 = .1, h2 = .3, H = .8 (e): h1 = .2, h2 = .3, H = .4

(e)

93

Figure 4.7: Estimates of θ

Figure 4.7 (a) h1 = .6, h2 = .7

Figure 4.7 (b) h1 = .6, h2 = .8

94

H =.75(red)

.85(green)
.95(blue)

H =.85(red)
.9(green)

Figure 4.7 (cont’d)

Figure 4.7 (c) h1 = .8, h2 = .9

Figure 4.7 (d) h1 = .3, h2 = .8

95

H =.93(red)

.98(green)

H =.85(red)

.95(green)

Figure 4.7 (cont’d)

Figure 4.7 (e) h1 = .1, h2 = .3

Figure 4.7 (f) h1 = .6, h2 = .9

96

H =.4(red)

.6(green)
.8(blue)

H =.95(red)

.99(green)

Figure 4.8: Estimates of h1, h2

Figure 4.8 (a) h1 = .6, h2 = .7

Figure 4.8 (b) h1 = .6, h2 = .8

97

H =.75(red)

.85(green)
.95(blue)

H =.85(red)
.9(green)

Figure 4.8 (cont’d)

Figure 4.8 (c) h1 = .8, h2 = .9

Figure 4.8 (d) h1 = .3, h2 = .8

98

H =.93(red)

.98(green)

H =.85(red)

.95(green)

Figure 4.8 (cont’d)

Figure 4.8 (e) h1 = .1, h2 = .3

Figure 4.8 (f) h1 = .6, h2 = .9

99

H =.4(red)

.6(green)
.8(blue)

H =.95(red)

.99(green)

Figure 4.9: Estimates of θ in a sample surface

Figure 4.10: Estimates of a1 and a2 in a sample surface with m = 1, 2, 3, 4

100

Figure 4.11: Estimates of a1 and a2 in a sample surface with m = 5, 6, 7, 8

Figure 4.12: Estimates of θ

101

Figure 4.13: Estimates of h1 and h2 in an fixed domain with m = 1, 2, 3, 4

Figure 4.14: Estimates of h1 and h2 in an fixed domain with m = 3, 4, 5, 6

102

Figure 4.15: Estimates of θ

Figure 4.16: Estimates of h1 and h2 in an fixed domain with m = 1, 2, 3, 4

103

Figure 4.17: Estimates of h1 and h2 in an fixed domain with m = 3, 4, 5, 6

104

CHAPTER 5

CONCLUSION AND DISCUSSION

In this dissertation, statistical inference was made on self-similar processes/random fields with
stationary increments. First, Hurst parameter, which captures self-similarity, was estimated in
multivariate self-similar stochastic processes and random field, OFBM and OSGRF, respectively.
Second, the dependency within multivariate random fields with stationary increments was measured
in the spectral domain.

Hurst estimation method with wavelet transform and eigen decompositon was examined thor-
oughly in a OFBM in continuous sample path, discrete sample path, and discrete noisy data. It was
revealed both theoretically and empirically that there is an interplay of Hurst parameters, h1, h2,
the covariance matrix Γ(1, 1) of XH(1), and the choice of the scale parameter j of wavelet function
on the performance of the estimator. If continuous sample paths are observed, then the bigger j is,
the better the estimators of h2, h1 perform. If discrete sample paths are observed, then it is better
to choose slightly smaller j than the maximum j = log2 n. Especially if hi is small e.g. hi < .5,
or determinant of Γ(1, 1) is small, then the estimators of h1, h2 have large bias for j close to its
maximum. However, the larger j is, the smaller the standard error of the estimators is. Considering
both bias and standard error, it is best to choose a set of j that is slightly smaller than log2 n or wide
range of j except the very large or small j when discrete sample paths are given. If noise is present
in the processes, then it is better to choose much smaller j than in the previous cases, but with the
noise term, the estimator does not work well unless both h1, h2 are small. For the estimator of θ, it
is always the best to choose the largest j = log2 n, since it has the both smallest bias and standard
error among all possible j in any set of (θ, h1, h2).

Hurst estimator was developed in OSGRF, and its performance was analyzed for different
sampling regimes. The performance of the estimator is affected by the several factors such as values
of the parameters, whether it is sampled in grid lines or it is sampled along the exact direction, and
whether samples are on a fixed domain or an increasing domain. It was observed that, in a fixed

105

domain setting, the estimator of a1(or h1) always performs well, whereas the estimator of a2(or h2)
performs well if the parameters are in the region of {(a1, a2); a1 > .5, a2 − a2 < .25} provided that
samples are obtained on the exact directions. In a fixed domain, if samples are on grid lines, then
the different choice of scale parameter m used in discrete variation method affects the performance
of the estimator ˆa2(or ˆh2) differently for a different set of (a1, a2). If the parameters fall in the
region {(a1, a2); a1 > .5, a2 − a2 < .25}, higher values of m works better for ˆa2(or ˆh2), whereas
lower values of m work better for ˆa2(or ˆh2) if the parameter a2 is less than .5. In an increasing
domain with H = 1, the estimators ˆa1, ˆa2 perform well in any range of (a1, a2).

The concept of coherence was extended and defined in multivariate random fields with stationary
increments, and its estimator was developed. The estimator was applied to OFBM with sample
paths observed in a fixed domain, and its behavior was observed both theoretically and empirically.
It was revealed that OFBM in a fixed domain has different dependence structure for different
matrices D(H), C in (3.23). Particularly, if the scaling matrix D(H) is diagonizable, not diagonal,
then the squared coherence function converges to constant function of 1 as the sample size grows in
a fixed domain, which implies that there exists strong correlation between two stochastic processes.
Moreover, in that case, the squared coherence converges to 1 faster if the two Hurst parameters
have a bigger difference (i.e., the bigger α2 − α1(h1 − h2) is, the stronger correlation exists in two
series.) If both C and D are diagonal matrix, then the squared coherence is zero function which
means that the two series are independent, having no correlation. If D is diagonal matrix, but not
C , then the squared coherence is close to constant function cγ, whose value is in between 0 and 1,
cγ =

C12√
C11C22
Hurst estimation in multivariate stochastic processes and random field is not only an interesting
topic but also has practical value as self-similarity and, more generally, fractality are seen in many
objects in nature. For future work, a better Hurst estimation method needs to be developed in
OFBM when noise is present, since the current method in this dissertation works well only for a
small range of parameters when noise is present. Hurst estimation in OFBM with various form of
Hurst matrix H is also of interest, especially OFBM with Jordan form H. A better Hurst estimation

, C(cid:96),k = (C C (cid:48))(cid:96),k, (cid:96), k = 1, 2.

106

method for h2 in OSGRF is needed as the method in this dissertation works well only for some
range of parameters. The coherence function defined in this dissertation is a useful measure in
capturing dependence structure in two random fields that are increment stationary. In future work,
further applications of the coherence in various random fields will be of interest as it will reveal
more of its practical value.

107

APPENDIX

108

APPENDIX

A.1 Appendix

Theorem 11 (Complex-valued Gaussian random measure) Let µ be a symmetric measure on
(RN, B(RN)) and let G = {A ∈ B(RN) : µ(A) < ∞}. Then there exists a complex-valued

Gaussian process indexed by G, i.e. (cid:101)Wµ = {(cid:101)Wµ(A) : A ∈ G}, written as(cid:101)Wµ(A) = W1(A) + iW2(A)

that satisfies the following conditions,

(cid:16)(cid:101)Wµ(A)(cid:101)Wµ(B)(cid:17)
• for any A ∈ G,(cid:101)Wµ(−A) =a.s. (cid:101)Wµ(A) ⇐⇒ W1(A) =a.s. W1(−A) and W2(A) =a.s. −W2(−A).
In general, we denote the Gaussian process(cid:101)Wµ by(cid:101)W.

• E

= µ(A ∩ B);

Proof. We only give the sketch of the proof. Put

E+ = {x ∈ RN : x1 ≥ 0, xi ∈ R, i = 2, . . . , N},
E− = {x ∈ RN : x1 ≤ 0, xi ∈ R, i = 2, . . . , N},

then (−E−) = E+ and RN = E+ ∪ E−. Without loss of the generality, we assume µ(E+ ∩ E−) = 0.
As in the real case, we can take two independently scattered Gaussian random measures W1 and
W2 on (E+, E+, 1

4 µ). Define, for A ∈ G,

(cid:101)W1(A) := W1(A ∩ E+) + W1(−(A ∩ E−)),
(cid:101)W2(A) := W2(A ∩ E+) + W2(−(A ∩ E−)),

and

(cid:101)W(A) :=(cid:101)W1(A) + i(cid:101)W2(A).

Then {(cid:101)W(A), A ∈ G} is a complex-valued Gaussian process with zero mean. We now verify the two

conditions. Notice that the distinction of the symbols between the conditions and constructions,

109

then the second condition is obvious. It is easy to see that

E((cid:101)W(A)(cid:101)W(B)) =E((cid:101)W1(A)(cid:101)W1(B) +(cid:101)W2(A)(cid:101)W2(B))
+ i E((cid:101)W1(B)(cid:101)W2(A) −(cid:101)W1(A)(cid:101)W2(B))
1
2 µ(AB) +

1
2 µ(AB) = µ(AB),

=

which yields the first condition.

(cid:3)

110

BIBLIOGRAPHY

111

BIBLIOGRAPHY

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

Achard, P. and Coeurjolly, J.-F. (2010) Discrete variations of the fractional Brownian motion
in the presence of outliers and an additive noise. Statistics Surveys, 4, 117-147.
Amblard, P.-O. and Coeurjolly, J.-F. (2011) Identification of the multivariate fractional
Brownian motion. IEEE Transactions on Signal Processing, 59, 5152-5168.
Amblard, P.-O., Coeurjolly, J.-F. and Philippe, A. (2012) Basic properties of the multivariate
fractional Brownian motion. arXiv:1007.0828
Arby, P. and Didier, G. (2018) Wavelet estimation for operator fractional Brownian motion.
Bernoulli, 24, 895-928.
Arcones, M. (1994) Limit theorems for nonlinear functionals of a stationary Gaussian
sequence of vectors. The Annals of Probability, 22, 2242-2274.
Bardet, J.-M. and Surgailis, D. (2013) Nonparametric estimation of the local Hurst function
of multifractional Gaussian processes. Stochastic Processes and their Applications, 123,
1004-1045.
Benson, D.A, Meerschaert, M.M, Baeumer, B and Scheffler, H.-P. (2006) Aquifer operator-
scaling and the effect on solute mixing and dispersion. WATER RESOURCES RESEARCH,
42, W01415, doi:10.1029/2004WR003755.
Biermé, H, Meerschaert, M.M, and Scheffler, H.-P. (2007) Operator scaling stable random
fields. Stochastic Processes and their Applications, 117, 312-332.
Biermé, H. and Lacaux, C. (2018) Fast and exact synthesis of some operator scaling Gaussian
random fields. Applied and Computational Harmonic Analysis.

[10] Breton, J.-C. and Coeurjolly, J.-F. (2012) Confidence intervals for the Hurst parameter of
a fractional Brownian motion based on finite sample size. Stat Inference Stoch Pocess, 15,
1-26.

[11] Brockwell, P, and Davis, R. (1991) Time series : Theory and Methods. Springer Series in

Statistics.

[12] Chan, G., and Wood, A.T.A (1999) Simulation of stationary Gaussian vector fields. Statistics

and Computing, 9, 254-268.

[13] Coeurjolly, J.-F., Amblard, P.-O. and Achard, S. (2013) Wavelet analysis of the multivariate
fractional Brownian motion . ESAIM: Probability and Statistics, EDP sciences, 17, 592-604.
[14] Coeurjolly J.-F. (2008) Hurst exponent estimation of locally self-similar Gaussian processes

using sample quantiles. The Annals of Statistics, 36, 1404-1434.

112

[15] Coeurjolly J.-F. (2001) Estimating the parameters of a fractional Brownian motion by discrete

variations of its sample paths. Statistical Inference for Stochastic processes, 4, 199-227.

[16] Cramér, H. (1940) On the theory of stationary random processes. Ann. Math. 41, 215–230.
[17] Delgado, R. (2007) A reflected fBm limit for fluid models with ON/OFF sources under heavy

traffic . Stochastic Processes and their Applications, 117, 188-201.

[18] Delbeke, L. and Arby, P. (2000) Stochastic integral representation and properties of the
wavelet coefficients of linear fractional stable motion Stochastic Processes and their Appli-
cations, 86, 177-182.

[19] Didier, G. and Pipiras, V. (2011) Integral representations and properties of operator frac-

tional Brownian motions. Bernoulli, 17, 1-33.

[20] Düker, M.-C. (2018) Limit theorems for Hilbert space-valued linear processes under long

range dependence. Stoch. Proc. Appl. 128 , 1439–1465.

[21] Düker, M.-C. Limit theorems for multivariate long-range dependent processes. ArXiv:

1704.08609.

[22] Fay, G., Moulines, E., Roueff, F. and Taqqu, M. (2009) Estimators of long-memory: Fourier

versus wavelets. Journal of Econometrics, 151, 159-177.

[23] Flandrin, P. (1989) On the spectrum of fractional Brownian motions. IEEE Transactions on

Information Theory, 35, 197-199.

[24] Flandrin, P. (1992) Wavelet analysis and synthesis of fractional Brownian motion. IEEE

Transactions on Information Theory, 38, 910-917.

[25] Fox, R. and Taqqu, M. (1986) Large-Sample Properties of Parameter Estimates for Strongly

Dependent Stationary Gaussian Time Series. The Annals of Statistics, 14, 517-532.

[26] Gloter, A. and Hoffmann, M. (2007) Estimation of the hurst parameter from discrete noisy

data. The Annals of Statistics, 5, 1947-1974.

[27] Gikhman, I. I. and Skorokhod, A. V. (2004) The Theory of Stochastic Processes. I. Springer-

Verlag, Berlin.

[28] Helgason, H., Pipiras, V. and Arby, P. (2011) Fast and exact synthesis of stationary multi-
variate Gaussian time series using circulant embedding. Signal Processing, 91, 1123-1133.
Itô, K. (1954) Stationary random distributions. Mem. Coll. Sci. Univ. Kyoto. Ser. A. Math.
28 , 209–223.

[29]

[30] Kleiber, W. (2017) Coherence for multivariate random fields. Statist. Sinica 27, 1675–1697.
[31] Kleiber, W. and Nychka, D.W.(2015) Equivalent kriging. Spat. Stat. 12 , 31–49.
[32] Lim. C, and Stein, M.L. (2008) Properties of spatial cross-periodograms using fixed-domain

asymptotics. J. Multivar. Anal. 99, 1962–1984.

113

[33] Lim, C.Y. and Meerschaert , M.M, and , Scheffler, H.-P. (2014) Parameter estimation for

operator scaling random fields. Journal of Multivariate Analysis, 123, 172–183.

[34] Maejima, M. and Mason, J.D. (1994) Operator-self-similar stable processes. Stoch. Process.

Appl. 54, 139–163.

[35] Majewski, K. (2005) Fractional Brownian heavy traffic approximations of multiclass feed-

forward queueing networks. Queueing systems, 50, 199-230.

[36] Mandelbrot, B. and Van Ness, J. (1968) Fractional Brownian motions, fractional noises and

applications. SIAM Review, 10, 422-437.

[37] Marinucci, D. and Robinson, P. M. (2000) Weak convergence of multivariate fractional

processes. Stoch. Process. Appl. 86 , 103–120.

[38] Mason, J.D. and Jurek, Z.J. (1993) Operator-limit distributions in probability theory. Wiley

Series in Probability and Statistics.

[39] Mason, J.D. and Xiao, Y. (2002) Sample path properties of operator-self-similar Gaussian

random fields. Theory Probab. Appl. 46 , 58–78.

[40] Masry, E. (1993) The wavelet transform of stochastic processes with stationary increments
and its application to fractional Brownian motion. IEEE Transactions on Information Theory,
39, 260-264.

[41] Meerschaert, M.M. and Scheffler, H.-P. (2001) Limit Distributions for Sums of Independent
Random Vectors. Heavy tails in theory and practice. John Wiley & Sons, Inc., New York.
[42] Moulines, E., Roueff, F. and Taqqu, M. (2007) Central limit theorem for the log-regression
wavelet estimation of the memory parameter in the Gaussian semi-parameter context. Frac-
tals, 15, 301-313.

[43] Moulines, E., Roueff, F. and Taqqu, M. (2007) On the spectral density of the wavelet
coefficients of long-memory time series with application to the log-regression estimation of
the memory parameter. Journal of Time series Analysis, 28, 155-187.

[44] Pitt, L.D. (1978) Scaling limits of Gaussian vector fields. J. Multivar. Anal. 8 , 45–54.
[45] Račkauskas, A. and Suquet, C. (2011) Operator fractional brownian motion as limit of

polygonal lines processes in Hilbert space. Stoch. Dyn. 11 , 49–70.

[46] Robinson, P. M. (1995) Gaussian semiparametric estimation of long range dependence. The

Annals of Statistics, 23, 1630-1661.

[47] Robinson, P. M. (1995) Log-periodogram regression of time series with long range depen-

dence. The Annals of Statistics, 23, 1048-1072.

[48] Roueff, F. and Taqqu, M. (2009) Asymptotic normality of wavelet estimators of the memory

parameter for linear processes. Journal of Time series Analysis, 30, 534-558.

114

[49] Tewifik, A.H. and Kim , M.(1992) Correlation Structure of the Discrete Wavelet Coefficients
of Fractional Brownian Motion. IEEE Transactions on Information Theory, 38, 904 - 909.
[50] Veitch, D. and Arby, P. (1999) A wavelet-Based joint estimator of the parameters of long-

range dependence. IEEE Transactions on Information Theory, 45, 878-897.

[51] Wendt, H., Didier, G., Combrexelle, S. and Abry, P. (2017) Multivariate Hadamard self-

similarity: testing fractal connectivity. Physica D: Nonlinear Phenomena 356/357, 1–36.

[52] Yaglom, A.M. (1957) Certain types of random fields in n-dimensional space similar to

stationary stochastic processes. (Russian) Teor. Veroyatnost. i Primenen 2 , 292–338.

115