PARAMETER ESTIMATION FOR GAUSSIAN RANDOM FIELDS AND MULTIVARIATE GAUSSIAN RANDOM PROCESSES UNDER FIXED-DOMAIN ASYMPTOTICS By Haoxiang Feng A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Statistics—Doctor of Philosophy 2024 ABSTRACT This dissertation explores parameter estimation for Gaussian random fields and multivari- ate Gaussian random processes under fixed-domain asymptotics, a crucial framework for modeling spatial and temporal data. Unlike increasing-domain asymptotics, fixed-domain asymptotics involve a growing number of observations within a fixed, bounded region, lead- ing to denser data. This scenario is common in applications such as image processing, where the spatial domain is constrained by the finite size of the sensor array. First, we study the parameter estimation for a Gaussian field with a multiplicative co- variance function, which is particularly relevant in computer experiments. We propose an increment-based estimator for estimating variance and scale parameters, and subsequent analysis shows that the estimator is both strongly consistent and asymptotically normal. Next, we extend the analysis to the bivariate Ornstein-Uhlenbeck process, constructing an explicit estimator that is strongly consistent and asymptotically normal. This estimator, requiring no prior parameter information, is shown to have the same asymptotic covariance matrix as that of the maximum likelihood estimator (MLE). Finally, we investigate asymptotic properties of MLE for the isotropic powered exponen- tial field. Unlike the Matérn model, the spectral density of the powered exponential model poses analytical challenges. We also establish conditions for the equivalence of Gaussian measures, providing a contrast to the orthogonality conditions found in earlier studies. Copyright by HAOXIANG FENG 2024 This dissertation is dedicated to my parents. iv ACKNOWLEDGEMENTS The research in this dissertation was partially supported by the NSF grant DMS-2153846. I would like to express my deepest gratitude to my advisor, Dr. Yimin Xiao, for his guidance, support, and patience throughout my PhD research. My sincere thanks also go to Dr. Haolei Weng, Dr. Chih-Li Sung, and Dr. Antoine Ayache for serving on my guidance committee and for their invaluable feedback and encouragement. I would also like to thank my friends for their support. I am grateful to the Department of Statistics and Probability and Michigan State Uni- versity for their support during my time in East Lansing. Lastly, I extend my heartfelt thanks to my beloved Yaqian for her enduring love and understanding. v TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 2 A CLASS OF ORNSTEIN-UHLENBECK FIELDS . . . . . . . . CHAPTER 3 THE MULTIVARIATE ORNSTEIN-UHLENBECK PROCESS . . CHAPTER 4 THE POWERED EXPONENTIAL FIELD . . . . . . . . . . . . . BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 6 49 58 80 vi CHAPTER 1 INTRODUCTION Multivariate Gaussian random fields defined on Rd are widely used to model spatial and temporal data. When fitting a random field to the data, it is common to assume that its covariance function belongs to a family of functions parameterized by a set of parameters which are then estimated based on the data. This turns the modeling problem into a pa- rameter estimation one (see, e.g., Chilès and Delfiner 2012, for an introduction to spatial statistical techniques). Three prevalent asymptotic frameworks exist in spatial statistics: increasing-domain asymptotics, fixed-domain asymptotics (see, e.g., M. L. Stein 1999a, chap. 3.3) and mixed- domain asymptotics (see, e.g., Lahiri 2003). In increasing-domain asymptotics, the distance between neighboring observation points remains above a positive threshold, causing the sam- pling region to expand as the number of observations increases. In fixed-domain asymptotics which is also called infill asymptotics in Cressie (1993), the number of observation points rises within a bounded sampling domain, resulting in increasingly dense observation points. In mixed-domain asymptotics, the sampling region expands, and at the same time, increasingly dense observation points fill in any given subregion of the sampling region. Fixed-domain asymptotics can occur in the process of spatial data. For example, we can model the analog signal a camera receives as Xa(s) = f (s) + X(s), s ∈ T, where T ⊂ R2 is a fixed and bounded set because the charge-coupled device (CCD) of a camera has a finite surface. We use f (s) to represent the original signal, and X(s) is the intrinsic physical noise which is modeled by a random field. Subsequently, an analog-to- digital converter converts Xa(s) to a digital signal, which can be modeled as Xd(n) = ⟨Xa(s), ϕn(s)⟩, 0 ≤ n ≤ N, 1 where {ϕn(s)} are the sensor responses. Here the two-dimensional vector n represents the pixel and N denotes the image size. If we plan to study the performance of a denoising algorithm based on measurements {Xd(n)} as N increases, it seems reasonable to adopt the fixed-domain asymptotic framework (see, e.g., Mallat 2008, chap. 11 for denoising). In the theory of fixed-domain asymptotics for multivariate Gaussian fields, there is a key concept called the equivalence of Gaussian measures. Specifically, for two probability measures P0 and P1 on a measurable space (Ω, F), say that P0 is absolutely continuous with respect to P1 if for all A ∈ F, P1(A) = 0 implies P0(A) = 0. Define P0 and P1 to be equivalent, written P0 ≡ P1, if they are mutually absolutely continuous. Define P0 and P1 to be orthogonal, written P0 ⊥ P1, if there exists A ∈ F such that P1(A) = 0 implies P0(A) = 1. Feldman (1958) proved that two Gaussian measures are either equivalent or orthogonal. The concept plays an important role in both prediction and estimation of multivariate Gaussian fields under fixed-domain asymptotics (see, e.g., M. L. Stein 1988; M. Stein 1990; M. L. Stein 1993; M. L. Stein 1999b; M. L. Stein 2004). Since this work focuses mainly on parameter estimation problems of multivariate Gaus- sian processes and univariate Gaussian fields under fixed-domain asymptotics, we first review recent results and relevant methods in the literature. On the one hand, Mason and Xiao (2002) studied sample path properties of operator fractional Brownian motions. Amblard and Coeurjolly (2011) studied parameter estimation of multivariate fractional Brownian mo- tion with an increment-based method and proved the strong consistency and asymptotic normality under increasing-domain asymptotics. Didier and Pipiras (2011) provided the integral representations of operator fractional Brownian motions in the spectral and time domains, respectively. Subsequently, Abry and Didier (2018) constructed estimators for op- erator fractional Brownian motion with wavelets, and showed the estimators are consistent and asymptotically normal under increasing-domain asymptotics. It is of interest to note that the corresponding high-pass filter of the wavelet used by Abry and Didier (2018) is a type of increment. On the other hand, Zhou and Xiao (2018) estimated the fractal indices of bi- 2 variate stationary Gaussian processes with the increment-based method under fixed-domain asymptotics and proved the consistency and asymptotic normality. Gneiting, Kleiber, and Schlather (2010) introduced a flexible parametric family of matrix-valued covariance func- tions for multivariate stationary Gaussian random fields, where each constituent component is a Matérn field. Zhang and Cai (2015) gave a sufficient condition for two probability measures corresponding to matrix-valued Matérn covariance functions to be equivalent and displayed an explicit example where cokriging is identical to kriging (best linear unbiased prediction). Velandia et al. (2017) showed that the MLE for a bivariate Ornstein-Uhlenbeck process is strongly consistent and asymptotically normal under fixed-domain asymptotics given some prior information of the parameter. Bachoc et al. (2022) provided conditions for equivalence of measures associated with multivariate Gaussian random fields and studied misspecified cokriging prediction for multivariate Matérn and generalized Wendland fields under fixed-domain asymptotics. Compared with MLE, the increment-based estimator can achieve linear computational complexity, and it is easily computed with no maximization required in many cases. Unsur- prisingly, there has been a long history of using increment-based methods to study properties of random fields. Actually, it can be traced back to the quadratic variation theorem in Lévy (1940). Since then, there has been a growing number of papers focused on parameter estima- tion for univariate random fields with increments under fixed-domain asymptotics. Kent and Wood (1997) estimated the fractal index of Gaussian processes using increments. Chan and Wood (2000) consistently estimated the fractal dimension with asymptotic normality using increments of the random field observed on a regular grid in R2. Anderes (2010) proved both variance and scale parameters from the Matérn model can be consistently estimated when d > 4 based on the increment-based method, but asymptotic distributions were not given. The increment-based method was then generalized by Loh, Sun, and Wen (2021) and Loh (2015) to estimate the smoothness parameter of the Matérn field irregularly sampled on Rd. However, both papers did not study the asymptotic distribution of the estimators. 3 In Chapter 2, we propose a method based on findings from Chan and Wood (2000) and Loh (2015) for estimating the variance and scale parameters in the model considered in Ying (1993) and study both its strong consistency and asymptotic normality with irregu- larly spaced data under fixed-domain asymptotics. Models with a multiplicative covariance function like the one considered in Ying (1993) are popular in computer experiments; (see, e.g., Sacks, Welch, et al. 1989; Sacks, Schiller, and Welch 1989; Paulo 2005; Bayarr et al. 2009; Peng and Wu 2014). In contrast to the increasing-domain asymptotics where the MLE of all (identifiable) parameters is consistent and asymptotically normal under some mild regularity conditions (Mardia and Marshall 1984), there is no general result for the asymptotic properties of MLE under fixed-domain asymptotics. However, there is quite a bit of literature on fixed- domain asymptotics of MLE when assuming that the covariance belongs to a parametric family. For the univariate Ornstein-Uhlenbeck process, Ying (1991) showed the product of its variance and scale parameters can be consistently estimated by MLE and the estimator is asymptotically normal. Regarding the isotropic Matérn model, with known smoothness parameter ν and free variance σ2 and scale parameter α, Zhang (2004) showed the MLE of σ2α2ν is strongly consistent for dimension d ≤ 3, and Du, Zhang, and Mandrekar (2009) showed the asymptotic normality of the estimator when d = 1. Wang and Loh (2011) extended the asymptotic result to dimension d ≤ 3. Bevilacqua et al. (2019) studied the MLE for the generalized Wendland model and derived similar results to those for the Matérn model. In Chapter 3, we construct an estimator for the bivariate Ornstein-Uhlenbeck process con- sidered in Velandia et al. (2017) and study its strong consistency and asymptotic normality. The construction is inspired by the results in Ying (1991) and Zhang (2004). Compared to the MLE from Velandia et al. (2017), our estimator has an explicit form and does not require any prior information of the parameter. Meanwhile, it turns out that the estimator has the same asymptotic covariance matrix as that of the MLE, echoing the fact that cokriging is 4 identical to kriging for this model as shown in Zhang and Cai (2015). In Chapter 4, we analyze the asymptotic properties of MLE for the isotropic powered ex- ponential field using tools from M. L. Stein (2004), Zhang (2004), and Wang and Loh (2011). Compared with the Matérn model, the spectral density of the powered exponential model cannot be analytically expressed except in some special cases, which brings new challenges to the analysis. Furthermore, we also establish the parameter condition for the equivalence of Gaussian measures, which contrasts with the orthogonality condition in Theorem 5 from Anderes (2010). 5 CHAPTER 2 A CLASS OF ORNSTEIN-UHLENBECK FIELDS 2.1 Introduction Stationary Gaussian fields with multiplicative covariance functions have been success- fully applied to the modeling of computer experiments. Sacks, Welch, et al. (1989) and Sacks, Schiller, and Welch (1989) proposed the use of a Gaussian field with the multiplica- tive powered exponential covariance function in their modeling of computer experiments. Observing the undesirable properties of the corresponding Gaussian field for modeling com- puter experiments, M. L. Stein (1989) proposed using a stationary Gaussian field model, X(u), u ∈ [0, 1]d, with mean 0 and the multiplicative Matérn covariance function, cov (X (u) , X (v)) =σ2 d (cid:89) i=1 21−ν Γ(ν) (θi|u[i] − v[i]|)νKν(θi|u[i] − v[i]|), ∀ u, v ∈ [0, 1]d, (2.1) where θ1, · · · , θd and σ2 are strictly positive parameters, and Kν is the modified Bessel function of the second kind with order ν > 0. By Equation (32) of M. L. Stein (1999a), the spectral density corresponding to Eq. (2.1) is f (ω) = σ2 d (cid:89) i=1 Γ(ν + 1 2)θ2ν i Γ(ν)π1/2 1 i + (ω[i])2)ν+(1/2) . (θ2 (2.2) The parameter ν controls the smoothness of the random field X. Specifically, X is m times mean square differentiable if and only if ν > m. If ν = n+ 1 2, Eq. (2.1) reduces to the product of an exponential function and a polynomial, in that cov (X (u) , X (v)) =σ2 d (cid:89) i=1 exp (−θi|u[i] − v[i]|) n (cid:88) k=0 (n + k)! (2n)! 6  (2θi|u[i] − v[i]|)n−k,  (2.3)    n k for n = 0, 1, · · · by Eq.(8.468) of Gradshteyn et al. (1981). Therefore, when ν = 0.5, Eq. (2.1) reduces to cov (X (u) , X (v)) = σ2 d (cid:89) exp (−θi|u[i] − v[i]|) . (2.4) i=1 In this regard, Ying (1993) showed both the strong consistency and the asymptotic normality of the maximum likelihood estimator (MLE) of (σ2, θ1, · · · , θd) when d ≥ 2. Subsequently, van der Vaart (1996) proved the MLE is also asymptotically efficient for dimension d = 2. When ν = 1.5, Eq. (2.1) reduces to cov (X (u) , X (v)) = σ2 d (cid:89) i=1 (1 + θi|u[i] − v[i]|) exp (−θi|u[i] − v[i]|) , (2.5) and Loh (2005) constructed a consistent estimator of (σ2, θ1, · · · , θd) based on the maximum likelihood method when d ≥ 3. In this chapter, we use the increment-based method from Chan and Wood (2000) and Loh (2015) to estimate the variance and scale parameters in Eq. (2.4) based on irregularly sampled data and study both its consistency and asymptotic normality under infill asymp- totics. For simplicity, we first study the dimension d = 2 in detail; therefore, there are only three parameters (λ, µ, σ2)′ being defined in our model below. Then we generalize the method to an arbitrary d in Section 2.5. Our main motivation to study this particular case is that, on the one hand, properties of the MLE in this case have been well studied so that a comprehensive comparison can be made between the MLE and the increment-based method in terms of the consistency and the asymptotic normality. On the other hand, one obvious advantage of the increment-based method is that the construction of the estimator does not involve the inverse of the covariance matrix compared to the MLE, which greatly reduces the computational complexity of the estimation. Actually, in our case, the computational complexity is O(n) where n is the sample size. However, to our best knowledge, the appli- cability of the increment-based method on a Gaussian field with a multiplicative covariance function under an irregular sampling scheme had not yet been studied. The proofs in this case lay the foundation for the asymptotic analysis of more general cases. 7 The rest of this chapter is organized as follows. In Section 2.2, we study the Ornstein- Uhlenbeck (OU) process with a covariance function Γµ(t) = σ2e−µ|t| by the increment-based method, which is the first step towards the asymptotic analysis of the OU field with the covariance function as Eq. (2.4). In Section 2.3, we establish the strong consistency and the asymptotic normality of the estimator of (λ, µ, σ2)′. In Section 2.4, we present a simulation study on the efficiency of the estimator in finite-sample cases. In Section 2.5, we generalize the increment-based method to any dimension d and derive the corresponding strong consistency and asymptotic normality results. We end the introduction with some notation. For any real-valued sequences an, bn, an ∼ bn means limn→∞ bn/an = 1. We write j = (j[1], · · · , j[d]), with brackets used to denote components of j, and 0 ≤ j ≤ n is equivalent to 0 ≤ j[ℓ] ≤ n[ℓ] for ℓ = 1, · · · , d. If j ∈ Zd is a multi-index, we write |j| = (cid:80)d ℓ=1 |j[ℓ]|. 2.2 Sampling over R In this section, we explore the performance of the quadratic variation built on the incre- ment introduced in Section 2 in Loh (2015). For the Ornstein-Uhlenbeck process X with mean 0 and the covariance function Γµ(t) = σ2e−µ|t|, we aim to estimate σ2µ with observa- tions X(tn,1), · · · , X(tn,n), where 0 = tn,1 < tn,2 < · · · < tn,n−1 < tn,n = 1. For brevity, we write tn,i = ti and X(ti) = Xi, i = 1, · · · , n. For ℓ ∈ N+, and i = 1, · · · , n − ℓ, we plan to design the increment such that ℓ (cid:88) k=0 ai,ktq (i+k) =    0, ℓ! ∀ q = 0, · · · , (ℓ − 1), if q = ℓ, (2.6) where we use the convention 00 = 1. The system of equations can be expressed in the 8 following matrix form          t0 (i+0) t0 (i+1) t1 (i+0) ... t1 (i+1) ... · · · · · · . . . t0 (i+ℓ) t1 (i+ℓ) ... tℓ (i+0) tℓ (i+1) · · · tℓ (i+ℓ)                            ai,0 ai,1 ... ai,ℓ                   0 0 ... ℓ! . = (2.7) The matrix on the left hand side is a Vandermonde matrix and the last column of the inverse of the Vandermonde matrix is (cid:32) 1 0≤s≤ℓ,s̸=0(ti+0 − ti+s) , (cid:81) 1 0≤s≤ℓ,s̸=1(ti+1 − ti+s) (cid:81) , · · · , (cid:81) 1 0≤s≤ℓ,s̸=ℓ(ti+ℓ − ti+s) (cid:33)′ . (2.8) Therefore, ai,k = (cid:81) ℓ! 0≤s≤ℓ,s̸=k(ti+k − ti+s) , ∀ k = 0, · · · , ℓ. (2.9) Over here, we only consider the case ℓ = 1; therefore, the coefficients of the increment can be simplified as ai,0 = −1/∆i, ai,1 = 1/∆i, (2.10) where ∆i = ti+1 − ti; (see Kent and Wood 1997, for a full exposition on the effect of ℓ ). Then define the quadratic variation as where Vn = n−1 (cid:88) i=1 (∇Xi)2, ∇Xi = ai,0Xi + ai,1Xi+1. (2.11) (2.12) To study the limiting moments of Vn, we impose a regularity condition on sampling points {ti}n i=1 as follows. Condition 1. For n ≥ 2, define ti = φ((i − 1)/(n − 1)), i = 1, · · · n, where φ : R → R is a twice continuously differentiable function satisfying φ(0) = 0, φ(1) = 1 and min 0≤s≤1 φ(1)(s) > 0. 9 It follows from Condition 1 that there exist positive constants C0 and C1 such that 0 < C0/n ≤ min (ti+1 − ti) ≤ max (ti+1 − ti) ≤ C1/n. 1≤i≤n−1 1≤i≤n−1 With the Taylor expansion, we have So ex = 1 + x + exθx 2! x2, 0 < θx < 1. Γµ(t) = σ2 − σ2µ|t| + e−µ|t|θ 2! σ2µ2|t|2, where 0 < θ < 1 depends on (−µ|t|). Then E((∇Xi)2) = 1 (cid:88) 1 (cid:88) ai,k1ai,k2Γµ(t(i+k2) − t(i+k1)) k1=0 = −σ2µ k2=0 1 (cid:88) k1=0 1 (cid:88) 1 (cid:88) + k1=0 k2=0 1 (cid:88) k2=0 ai,k1ai,k2|t(i+k2) − t(i+k1)| ai,k1ai,k2 e−µ|t(i+k2)−t(i+k1)|θ 2! σ2µ2|t(i+k2) − t(i+k1)|2 (2.13) (2.14) (2.15) (2.16) = −2(σ2µ)ai,0ai,1(|ti+1 − ti|) + O(1). And the O(1) in the above equation comes from the fact that the last term in the second last equation can be uniformly bounded over i = 1, · · · , (n − 1), namely, 1 (cid:88) 1 (cid:88) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) k1=0 1 (cid:88) k2=0 (cid:12) 1 (cid:88) (cid:12) (cid:12) (cid:12) k1=0 1 (cid:88) k2=0 1 (cid:88) k1=0 k2=0 ≤ ≤ ai,k1ai,k2 e−µ|t(i+k2)−t(i+k1)|θ 2! ai,k1ai,k2 e−µ|t(i+k2)−t(i+k1)|θ 2! σ2µ2|t(i+k2) − t(i+k1)|2 σ2µ2|t(i+k2) − t(i+k1)|2 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ai,k1ai,k2 1 2 σ2µ2|t(i+k2) − t(i+k1)|2 (cid:12) (cid:12) (cid:12) (cid:12) (2.17) =(cid:12) (cid:12)ai,k0ai,k1∆2 i (cid:12) (cid:12)σ2µ2 =σ2µ2. 10 Meanwhile, by the Taylor expansion again, we have − 2(σ2µ) n−1 (cid:88) i=1 (ai,0ai,1|ti+1 − ti|) =2(σ2µ) =2(σ2µ) =2(σ2µ) n−1 (cid:88) i=1 n−1 (cid:88) i=1 n−1 (cid:88) i=1 (cid:18) 1/ φ( i n − 1 ) − φ( (cid:19) ) i − 1 n − 1 i − 1 n − 1  (cid:32) (cid:18) φ(1)( 1/ (n − 1) φ(1)( i−1 n−1) 1 − 1 + αi ) 1 (n − 1) + φ(2)( i − 1 n − 1 + θi n − 1 ) 1 2(n − 1)2 (cid:19) (2.18) n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) (cid:33)−2 n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1)   =2(σ2µ)n2 (cid:90) 1 0 {φ(1)(s)}−1ds + O(n), as n → ∞, where θi, αi ∈ (0, 1) depend on i. Notice that the second last equation above holds when n is big enough since n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) > −1, uniformly over i = 1, · · · , (n − 1) under this case. As for the O(n) term in Eq. (2.18), we only need to notice that for n sufficiently big, (cid:32) 1 + αi n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) (cid:33)−2 < 2, uniformly over i = 1, · · · , (n − 1). Therefore, (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (n − 1) φ(1)( i−1 n−1) (cid:32) 1 + αi (cid:33)−2 n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ maxx∈[0,1] |φ(2)(x)| minx∈[0,1] |φ(1)(x)|2 , (2.19) uniformly over i = 1, · · · , (n − 1) when n big enough. Therefore, with Eqs. (2.16) and (2.18), we have E(Vn) = n−1 (cid:88) i=1 E((∇Xi)2) = −2(σ2µ) = 2(σ2µ)n2 n−1 (cid:88) i=1 (cid:90) 1 0 ai,0ai,1(|ti+1 − ti|) + O(n) (2.20) {φ(1)(s)}−1ds + O(n), 11 as n → ∞. From Eq. (2.20), we expect (cid:18) Vn/ 2n2 {φ(1)(s)}−1ds (cid:19) (cid:90) 1 0 a.s.−−→ σ2µ, as n → ∞. So we first study Vn/ E(Vn) and then replace E(Vn) by 2n2 (cid:82) 1 estimator of σ2µ. 0 {φ(1)(s)}−1ds to get the desired Lemma 2.1. If the sampling function φ meets Condition 1, then lim n→∞ n var(Vn/ E(Vn)) = 2 (cid:82) 1 0 {φ(1)(s)}−2ds (cid:16)(cid:82) 1 0 {φ(1)(s)}−1ds (cid:17)2 ≜ ϕ0. Proof. By Eq. (2.11), we have var(Vn) = var (cid:33) (∇Xi)2 (cid:32)n−1 (cid:88) i=1 n−1 (cid:88) n−1 (cid:88) i=1 j=1 n−1 (cid:88) n−1 (cid:88) i=1 j=1 n−1 (cid:88) n−1 (cid:88) = = = cov((∇Xi)2, (∇Xj)2) 2[E(∇Xi∇Xj)]2 (cid:34) 1 (cid:88) 1 (cid:88) 2 ai,k1aj,k2Γµ(ti+k1 − tj+k2) . (cid:35)2 If i = j, similar to the derivation of E(Vn), we have i=1 j=1 k1=0 k2=0 (2.21) (2.22) Pn ≜ n−1 (cid:88) 2 i=1 (cid:34) 1 (cid:88) 1 (cid:88) ai,k1ai,k2Γµ(ti+k1 − ti+k2) (cid:35)2 (cid:34) n−1 (cid:88) =2 i=1  k1=0 k2=0 2(σ2µ)(n − 1) φ(1)( i−1 n−1) (cid:32) × 1 − 1 + αi n−1 (cid:88) =2 (cid:34) 2(σ2µ) i=1 (cid:20) (2σ2µn)2n =2 (n − 1) φ(1)( i−1 n−1) (cid:90) 1 n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1) (cid:35)2 + O(1) (cid:21) {φ(1)(s)}−2ds + O(n2) . 0 12 (cid:33)−2 n−1 + θi φ(2)( i−1 n−1) 2(n − 1)φ(1)( i−1 n−1)  (cid:35)2  + O(1) (2.23) If i ̸= j, because of the symmetry, we only need to study the case j > i, then tj+k2 −ti+k1 ≥ 0 for all k1, k2 = 0, 1. By the Taylor expansion (see, e.g., Loh 2015, Lemma 4), Γµ(tj+k2 − ti+k1) = 1 (cid:88) s=0 Γ(s) µ (tj − ti) s! [(tj+k2 − ti+k1) − (tj − ti)]s+ (cid:90) tj+k2 −ti+k1 tj −ti [(tj+k2 − ti+k1) − t]Γ(2) µ (t)dt. Meanwhile, regarding Γµ(t) = σ2e−µ|t|, Γ(2) µ (t) = σ2µ2e−µt for t ≥ 0. So ai,k1aj,k2Γµ(tj+k2 − ti+k1) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ai,k1aj,k2 (cid:32)(cid:90) tj+k2 −ti+k1 tj −ti [(tj+k2 − ti+k1) − t] Γ(2) µ (t)dt (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (2.24) |ai,k1aj,k2| (cid:18)(cid:90) I |(tj+k2 − ti+k1) − t| dt (cid:19) (cid:12) (cid:12)Γ(2) µ (t)(cid:12) (cid:12) max t∈I (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) 1 (cid:88) k1=0 1 (cid:88) k2=0 1 (cid:88) k1=0 1 (cid:88) k2=0 1 (cid:88) = ≤ k1=0 C 2 1 C 2 0 ≤4 k2=0 σ2µ2, where I ≜ ((tj − ti) ∧ (tj+k2 − ti+k1), (tj − ti) ∨ (tj+k2 − ti+k1)) and the last inequality holds by noticing that max k1,k2 |(tj+k2 − ti+k1) − (tj − ti)| ≤ max{∆i, ∆j} ≤ C1 n . (2.25) Therefore, Qn ≜ (cid:88) 2 (cid:34) 1 (cid:88) 1 (cid:88) ai,k1aj,k2Γµ(ti+k1 − tj+k2) = O(n2). (2.26) (cid:35)2 i̸=j Combining Eqs. (2.22), (2.23) and (2.26), we have k2=0 k1=0 lim n→∞ n var(Vn/ E(Vn)) = lim n→∞ n n−1 (cid:88) n−1 (cid:88) 2 (cid:34) 1 (cid:88) 1 (cid:88) i=1 j=1 k1=0 k2=0 ai,k1aj,k2Γµ(ti+k1 − tj+k2) /(E(Vn))2 (cid:35)2 (2.27) n(Pn + Qn)/f 2 n = lim n→∞ (cid:82) 1 0 (φ(1)(s))−2ds (cid:104)(cid:82) 1 0 (φ(1)(s))−1ds = 2 (cid:105)2 . 13 Remark 2.1. By Jensen’s inequality, (cid:82) 1 0 (φ(1)(s))−2ds [(cid:82) 1 0 (φ(1)(s))−1ds]2 ≥ 1, and the equality holds if and only if φ(1)(s) = const on [0, 1]. This implies the estimator has the same asymptotic variance as that of MLE only when the sampling scheme is regular. If we consider a family of sampling functions (cid:40) φi(x) = (cid:18) (x + i i + 2 1 i )2 − 1 i2 (cid:19)2 (cid:41) : i ∈ N+ , then (cid:82) 1 0 {φ(1) (cid:104)(cid:82) 1 0 {φ(1) i (s)}−2ds i (s)}−1ds (cid:105)2 → +∞, as i → ∞. Loh (2015) showed that Vn/ E(Vn) a.s.−−→ 1. In the following, we derive a weaker upper bound for the convergence rate based on Lemma 2.1, which will be used in the subsequent multidimensional analysis. Lemma 2.2. ∀ ξ > 0, there exist a finite constant C5 > 0 such that when n is sufficiently large, P (cid:18)(cid:12) (cid:12) (cid:12) (cid:12) Vn E(Vn) (cid:12) (cid:12) − 1 (cid:12) (cid:12) (cid:19) ≥ ξ ≤ 2e−C5ξ √ √ n/ ϕ0. (2.28) Proof. The proof follows the argument of Lemma 1 given in Zhou and Xiao (2018). Writing (cid:32) Yn = ∇X1 (cid:112)E(Vn) , · · · , ∇Xn−1 (cid:112)E(Vn) (cid:33)′ , Σn = (Σi,j)(n−1)×(n−1) = cov(Yn). (2.29) Let Λn = diag(λ1, · · · , λn−1) be the diagonal matrix whose diagonal entries are the eigenval- ues of Σn, and let U = (U1, · · · , Un−1)′ ∼ N (0, In−1). Then, we have Vn/ E(Vn) = Y′ nYn d= U′ΛnU. 14 We apply the Hanson and Wright inequality (Hanson and Wright 1971) to bound the tail probability of the quadratic forms and obtain P(|Vn/ E(Vn) − E(Vn/ E(Vn))| > ξ) = P(|U′ΛnU − tr(Λn)| > ξ) (cid:40) (cid:32) ≤2 exp − min C3 ξ ∥Λn∥2 , C4 ξ2 ∥Λn∥2 F (cid:33)(cid:41) , (2.30) where ∥Λn∥2 are positive constants independent of Λn and ξ. Note that and ∥Λn∥F are the ℓ2 norm and Frobenius norm of Λn, respectively; C3 and C4 ϕn ≜ var(Vn/ E(Vn)) = var(U′ΛnU) = 2tr(Λ2 n) = 2∥Λn∥2 F . By Lemma 2.1, we have ∥Λn∥2 F = 1 2ϕn ∼ 1 2nϕ0. Combining the above with the fact that ∥Λn∥2 ≤ ∥Λn∥F , we see that ≥ 1 ∥Λn∥2 1 ∥Λn∥2 F (cid:114) n ϕ0 , ∼ 2n ϕ0 as n → ∞. Hence, Eq. (2.30) decays n when n → ∞. Consequently, when n is sufficiently large, P(|Vn/ E(Vn) − 1| > ξ) ≤ 2e−C5ξ √ √ n/ ϕ0, when n is sufficiently large. Meanwhile exponentially with rate √ where C5 ≜ C3. Theorem 2.1. Let Vn be as in Eq. (2.11), and suppose Condition 1 holds. Then Vn/ E(Vn) a.s.−−→ 1, as n → ∞. (2.31) Proof. Combining Lemma 2.2 and the Borel-Cantelli lemma, we have the conclusion. 2.3 Sampling over R2 Let X(u), u ∈ [0, 1]2 be a zero-mean univariate Gaussian field with covariance function cov (X (u) , X (v)) = σ2 exp (−λ|u[1] − v[1]| − µ|u[2] − v[2]|) , (2.32) 15 where (σ2, λ, µ)′ ∈ R3 + . By Ying (1993), we know that σ2, λ and µ can be consistently estimated. Let φh and φv satisfy Condition 1. We consider the observation locations over R2 which are of the form u(i) ≜ (cid:18) φh (cid:18) i[1] − 1 m − 1 (cid:19) , φv (cid:18) i[2] − 1 n − 1 (cid:19)(cid:19)′ , where 1 ≤ i[1] ≤ m, 1 ≤ i[2] ≤ n, and n/m → ρ ∈ (0, ∞). For simplicity, we write Xi = X(u(i)). Since the observation region is in R2, we will use letters v and h to emphasize operations in the sense of the vertical and horizontal directions respectively. To simplify the notation, let φ(h,j) = φh (cid:18) j − 1 m − 1 (cid:19) , φ(v,l) = φv (cid:18) l − 1 n − 1 (cid:19) , for j = 1, · · · , m and l = 1, · · · , n. Let ∆h,j = φ(h,j+1) − φ(h,j), ∆v,k = φ(v,k+1) − φ(v,k), for j = 1, · · · , (m − 1) and k = 1, · · · , (n − 1). Define the index set I(m,n) = {i : (1, 1)′ ≤ i ≤ (m − 1, n − 1)′}, (2.33) If j ∈ I(m,n) is a vector, ∆h,j and ∆v,j are defined as ∆h,j = φ(h,j[1]+1) − φ(h,j[1]), ∆v,j = φ(v,j[2]+1) − φ(v,j[2]). Similar to the increment in Eq. (2.10), we define 4 sets of increments {ah;i,0, ah;i,1}, {av;i,0, av;i,1}, {ah;i,(0,0), ah;i,(1,0)}, {av;i,(0,0), av;i,(0,1)}, 16 as ah;i,0 = ah;i,(0,0) = −1/∆h,i, ah;i,1 = ah;i,(1,0) = 1/∆h,i, av;i,0 = av;i,(0,0) = −1/∆v,i, av;i,1 = av;i,(0,1) = 1/∆v,i. (2.34) (2.35) Based on Condition 1, we also have positive constants Ch,0 and Ch,1 such that 0 < Ch,0/m ≤ min 1≤i≤(m−1) (cid:18) φh (cid:18) i (cid:19) m − 1 − φh (cid:19)(cid:19) (cid:18) i − 1 m − 1 (cid:18) ≤ max 1≤i≤(m−1) φh (cid:18) i m − 1 (cid:19) − φh (cid:19)(cid:19) (cid:18) i − 1 m − 1 ≤ Ch,1/m. Cv,0 and Cv,1 can be similarly defined. 2.3.1 The vertical and horizontal increments Similar to the definition of Vn in Eq. (2.11) for the univariate case, for each fixed j = 1, · · · , (m − 1), let where Vv,j = n−1 (cid:88) k=1 (cid:0)∇vX(j,k) (cid:1)2, ∇vX(j,k) = av;(j,k),(0,0)X(j,k) + av;(j,k),(0,1)X(j,k+1). (2.36) (2.37) From the structure of the covariance function Eq. (2.32), the Gaussian process X(i,·) has the same distribution across 1 ≤ i ≤ m. Therefore, to estimate σ2µ, we first construct the weighted quantity ¯Zv,(m,n) = m−1 (cid:88) j=1 ∆h,jVv,j/ E(Vv,j) = m−1 (cid:88) j=1 ∆h,jVv,j/ E(Vv,1). (2.38) Switching the direction from vertical to horizontal, for each fixed j = 1, · · · , (n − 1), we can also define ∇hX(k,j) = ah;(k,j),(0,0)X(k,j) + ah;(k,j),(1,0)X(k+1,j), Vh,j = m−1 (cid:88) k=1 (cid:0)∇hX(k,j) (cid:1)2. 17 Therefore, to estimate σ2λ, we construct the weighted quantity ¯Zh,(m,n) = n−1 (cid:88) j=1 ∆v,jVh,j/ E(Vh,j) = n−1 (cid:88) j=1 ∆v,jVh,j/ E(Vh,1). (2.39) From the definitions of ¯Zv,(m,n) and ¯Zh,(m,n), it is reasonable to just study statistical properties of ¯Zv,(m,n). To analyze the variance of ¯Zv,(m,n), we first note that for any j, k = 1, · · · , (m−1), cov(Vv,j, Vv,k) = e−2λ|φ(h,j)−φ(h,k)| var(Vv,1). P (λ) =    1 e−λ|φ(h,j)−φ(h,k)| e−λ|φ(h,j)−φ(h,k)| 1   , Actually, denote and B(µ) = (cid:0)e−µ|φ(v,p)−φ(v,q)|(cid:1) n×n , where p, q = 1, · · · , n, then the 2n-vector Y2n ≜ (cid:0)(cid:0)X(j,l), l ∈ {1, · · · , n}(cid:1), (cid:0)X(k,l), l ∈ {1, · · · , n}(cid:1)(cid:1)′ has the distribution (2.40) (2.41) (2.42) Y2n ∼ N (0, σ2P (λ) ⊗ B(µ)), (2.43) where ⊗ is the Kronecker product. Let    1 0 0 0   , R1 =    0 0 0 1   . R2 = Since Vv,j is a quadratic form, it can be expressed as Vv,j = Y ′ 2n(R1 ⊗ F )Y2n, where F is a n × n semi-positive definite matrix. Similarly, Vv,k = Y ′ 2n(R2 ⊗ F )Y2n. 18 (2.44) (2.45) Therefore, cov(Vv,j, Vv,k) = cov (Y ′ 2n(R1 ⊗ F )Y2n, Y ′ 2n(R2 ⊗ F )Y2n) =2 tr(cid:0)(R1 ⊗ F )σ2(P (λ) ⊗ B(µ))(R2 ⊗ F )σ2(P (λ) ⊗ B(µ))(cid:1) (2.46) =2σ4 tr((R1P (λ)R2P (λ)) ⊗ (F B(µ)F B(µ))) =e−2λ|φ(h,j)−φ(h,k)|2σ4 tr((F B(µ)F B(µ))) =e−2λ|φ(h,j)−φ(h,k)| var(Vv,1). Based on Eqs. (2.40) and (2.46), we have var (cid:0) ¯Zv,(m,n) (cid:1) = var (cid:33) ∆h,jVv,j/ E(Vv,1) (cid:32)m−1 (cid:88) j=1 = var (cid:18) Vv,1 (cid:19) m−1 (cid:88) m−1 (cid:88) E(Vv,1) j=1 k=1 e−2λ|φ(h,j)−φ(h,k)|(∆h,j∆h,k). (2.47) Based on the definition of {∆h,j}m−1 j=1 and Condition 1, we have m−1 (cid:88) m−1 (cid:88) j=1 k=1 e−2λ|φ(h,j)−φ(h,k)|(∆h,j∆h,k) m→∞−−−→ (cid:90) (cid:90) [0,1]2 e−2λ|x−y| dx dy. (2.48) A direct calculation shows that (cid:90) (cid:90) Cλ ≜ [0,1]2 e−2λ|x−y| dx dy = 1 λ − 1 − e−2λ 2λ2 . (2.49) We are ready to prove the following theorem on asymptotic properties of ¯Zv,(m,n). Theorem 2.2. Let ¯Zv,(m,n) be as in Eq. (2.38). Suppose the two sampling functions φh(·) and φv(·) satisfy Condition 1. Then n var (cid:0) ¯Zv,(m,n) (cid:1) = 2Cλ lim n→∞ and (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) v (s)}−2ds v (s)}−1ds (cid:17)2 , ¯Zv,(m,n) a.s.−−−→ n→∞ 1. 19 (2.50) (2.51) Proof. Combining Lemma 2.1 and Eqs. (2.47) to (2.49), we have Eq. (2.50). As for the almost sure convergence, let ξ > 0, then by Lemma 2.2, P (cid:0)| ¯Zv,(m,n) − E( ¯Zv,(m,n))| > ξ(cid:1) (cid:33) (cid:32)m−1 (cid:88) ≤ P ∆h,i Vv,i E(Vv,i) (cid:12) (cid:12) − 1 (cid:12) (cid:12) > ξ (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ m−1 (cid:88) i=1 ≤m P i=1 (cid:18) P ∆h,i (cid:18)(cid:12) (cid:12) (cid:12) (cid:12) Vv,1 E(Vv,1) √ Vv,i E(Vv,i) (cid:12) (cid:12) − 1 (cid:12) (cid:12) √ ϕ0, ≤2me−C5 nξ/Ch,1/ (cid:19) > ξ/m (cid:12) (cid:12) − 1 (cid:12) (cid:12) (cid:19) (2.52) > ξ/Ch,1 when n is big enough. So applying the Borel-Cantelli lemma, we have the almost sure convergence. 2.3.2 The square increment In this part, we focus on estimating σ2µλ. Define the increment {bk}1 k=0 as b(0,0) = b(1,1) = 1, b(1,0) = b(0,1) = −1, and we call it the square increment since its support is the four vertices of the square; (see Chan and Wood 2000, for a fuller description of it) . In the following, we reserve the symbol s to emphasize the square increment. Based on the structure of the covariance function, for each fixed 1 ≤ j ≤ m, {X(j,l)}n l=1 contains information of σ2µ. And for each fixed 1 ≤ l ≤ n, {X(j,l)}m j=1 contains information of σ2λ. So naturally we define the quantity as Define the quadratic variation and its transformation ∇sXi = (cid:80)1 k=0 bkXi+k ∆h,i∆v,i . Vs,(m,n) = (cid:88) (∇sXi)2, i∈I(m,n) ¯Zs,(m,n) = Vs,(m,n)/ E (cid:0)Vs,(m,n) (cid:1), 20 (2.53) (2.54) (2.55) then for fixed index i, j ∈ I(m,n), E (∇sXi∇sXj) = 1 ∆h,i∆v,i∆h,j∆v,j = σ2 ∆h,i∆v,i∆h,j∆v,j 1 (cid:88) k,l=0 1 (cid:88) k,l=0 bkbl E (Xi+kXj+l) bkble−λ|φ(h,(i+k)[1])−φ(h,(j+l)[1])|−µ|φ(v,(i+k)[2])−φ(v,(j+l)[2])| (cid:32) 1 (cid:88) bk1bl1e−λ|φ(h,(i[1]+k1))−φ(h,(j[1]+l1))| (cid:33) bk2bl2e−µ|φ(v,(i[2]+k2))−φ(v,(j[2]+l2))| (cid:33) = σ2 ∆h,i∆h,j 1 ∆v,i∆v,j (cid:32) 1 (cid:88) × =σ2 l1,k1=0 (cid:32) 1 (cid:88) l2,k2=0 ah;i,kah;j,le−λ|φ(h,(i[1]+k))−φ(h,(j[1]+l))| (cid:33) (2.56) l,k=0 (cid:32) 1 (cid:88) l,k=0 × av;i,kav;j,le−µ|φ(v,(i[2]+k))−φ(v,(j[2]+l))| (cid:33) , where b0 = −1, b1 = 1. So we have decomposed the covariance between ∇sXi and ∇sXj into the product of the covariance in the univariate case. Based on the derivation of Eqs. (2.18) 21 and (2.20), we have E (cid:0)Vs,(m,n) (cid:1) = (cid:88) E (cid:2)(∇sXi)2(cid:3) i∈I(m,n) = (cid:88) σ2 i∈I(m,n) (cid:32) 1 (cid:88) × = σ2 ah;i,kah;i,le−λ|φ(h,(i[1]+k))−φ(h,(i[1]+l))| (cid:33) (cid:32) 1 (cid:88) l,k=0 av;i,kav;i,le−µ|φ(v,(i[2]+k))−φ(v,(i[2]+l))| (cid:33) l,k=0 (cid:32)m−1 (cid:88) (cid:32) 1 (cid:88) ah;i1,kah;i1,le−λ|φ(h,(i1+k))−φ(h,(i1+l))| (cid:33)(cid:33) l,k=0 i1=1 (cid:32) n−1 (cid:88) (cid:32) 1 (cid:88) av;i2,kav;i2,le−µ|φ(v,(i2+k))−φ(v,(i2+l))| (cid:33)(cid:33) × (2.57) i2=1 (cid:18) l,k=0 (cid:90) 1 = σ2 2λm2 0 (cid:90) 1 {φ(1) h (s)}−1ds + O(m) (cid:19) {φ(1) v (s)}−1ds + O(n) (cid:18) × 2µn2 (cid:19) 0 = 4σ2λµ(mn)2 as n → ∞. (cid:90) 1 0 {φ(1) h (s)}−1ds (cid:90) 1 0 {φ(1) v (s)}−1ds + O(n3), Theorem 2.3. Let ¯Zs,(m,n) be as in Eq. (2.55). Suppose the two sampling functions φh(·) and φv(·) satisfy Condition 1. Then mn var (cid:0) ¯Zs,(m,n) (cid:1) = 2 lim n→∞ and (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) v (s)}−2ds v (s)}−1ds (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) (cid:17)2 h (s)}−2ds h (s)}−1ds (cid:17)2 , Proof. As with the univariate case, the variance of Vs,(m,n) can be written as ¯Zs,(m,n) a.s.−−−→ n→∞ 1. var (cid:0)Vs,(m,n) (cid:1) = var   (cid:88) (∇sXi)2   2[E (∇sXi∇sXj)]2. i∈I(m,n) (cid:88) (cid:88) = i∈I(m,n) j∈I(m,n) 22 (2.58) (2.59) (2.60) As with the derivation of Eq. (2.23), based on Eqs. (2.17) and (2.19), we have that 1 (cid:88) l,k=0 av;i,kav;i,le−µ|φ(v,(i[2]+k))−φ(v,(i[2]+l))| = 2µ (n − 1) (cid:16) i[2]−1 n−1 φ(1) v (cid:17) + O(1), (2.61) as n → ∞ uniformly over i ∈ I(m,n). We first deal with the case i = j, P(m,n) ≜ (cid:88) 2(cid:8)E (cid:2)(∇sXi)2(cid:3)(cid:9)2 i∈I(m,n) (cid:88) = (cid:40) 2 σ2 (cid:32) 1 (cid:88) ah;i,kah;i,le−λ|φ(h,(i[1]+k))−φ(h,(i[1]+l))| (cid:33) i∈I(m,n) (cid:32) 1 × (cid:88) l,k=0 l,k=0 av;i,kav;i,le−µ|φ(v,(i[2]+k))−φ(v,(i[2]+l))| (cid:33)(cid:41)2 (cid:34) 2λ (m − 1) φ(1) h ( i[1]−1 m−1 ) (cid:35)2 (cid:34) + O(1) 2µ (cid:35)2 + O(1) (n − 1) v ( i[2]−1 n−1 ) φ(1) (cid:35)2  + O(1) (2.62) = 2σ4 (cid:88) i∈I(m,n)  (cid:34)  m−1 (cid:88) = 2σ4  i1=1 (cid:34) n−1 (cid:88)   × 2µ  i2=1 (cid:20) = 2σ4 (2λm)2m (cid:20) × (2µn)2n (cid:90) 1 0 2λ (m − 1) φ(1) h ( i1−1 m−1) (n − 1) φ(1) v ( i2−1 n−1 ) (cid:90) 1  (cid:35)2   + O(1) (cid:21) h (s)}−2ds + O(m2) {φ(1) 0 (cid:21) v (s)}−2ds + O(n2) {φ(1) = 2(4σ2λµ)2(mn)3 (cid:90) 1 0 {φ(1) h (s)}−2ds (cid:90) 1 0 {φ(1) v (s)}−2ds + O(n5), as n → ∞. Secondly, we consider the case where i and j are different but in the same row or column. Because of the similarity between the row and the column case, we only show the case 23 j[1] = i[1] and j[2] ̸= i[2]. Based on Eq. (2.24), we have Q(m,n) ≜ (cid:88) 2{E [∇sXi∇sXj]}2 i∈I(m,n) j[1]=i[1] j[2]̸=i[2] (cid:88) = (cid:40) 2 σ2 (cid:32) 1 (cid:88) ah;i,kah;i,le−λ|φ(h,(i[1]+k))−φ(h,(i[1]+l))| (cid:33) i∈I(m,n) j[1]=i[1] j[2]̸=i[2] (cid:32) 1 (cid:88) l,k=0 av;i,kav;j,le−µ|φ(v,(i[2]+k))−φ(v,(j[2]+l))| (cid:33)(cid:41)2 × ≤ n l,k=0 (cid:88) (cid:40) 2 σ2 (cid:32) 1 (cid:88) i∈I(m,n) l,k=0 ah;i,kah;i,le−λ|φ(h,(i[1]+k))−φ(h,(i[1]+l))| (cid:33) (2.63) (cid:21)(cid:19)(cid:41)2 (cid:18)(cid:20) 4 × C 2 v,1 C 2 v,0 µ2 = 2σ4n (cid:88) (cid:34) 2λ i∈I(m,n)  (cid:34)  m−1 (cid:88) = 2σ4n  i1=1 = 8σ4Cqλ2m3n2 2λ (m − 1) φ(1) h ( i1−1 m−1) (cid:90) 1 (m − 1) φ(1) h ( i[1]−1 m−1 ) (cid:35)2 (cid:20) + O(1) (cid:21)2 µ2 4 C 2 v,1 C 2 v,0 (cid:35)2  (cid:40) n−1 (cid:88) + O(1)  i2=1 (cid:21)2(cid:41) (cid:20) 4 C 2 v,1 C 2 v,0 µ2 {φ(1) h (s)}−2ds + O(n4), as n → ∞, where 0 C 2 v,1 C 2 v,0 Finally, we consider the case j[1] ̸= i[1] and j[2] ̸= i[2]. Based on Eq. (2.24), it should be Cq ≜ 4 µ2. easy to see that as n → ∞. R(m,n) ≜ (cid:88) 2{E [∇sXi∇sXj]}2 = O(n4), (2.64) i∈I(m,n) j[1]̸=i[1] j[2]̸=i[2] Combining Eqs. (2.62) to (2.64), we see that var (cid:0)Vs,(m,n) (cid:1) = 2(4σ2λµ)2(mn)3 (cid:90) 1 0 {φ(1) h (s)}−2ds (cid:90) 1 0 {φ(1) v (s)}−2ds + O(n5), (2.65) 24 as n → ∞. Consequently, it follows from Eqs. (2.57) and (2.65) that Eq. (2.58) holds. With Eq. (2.58) on hand, following the proof of Theorem 2.1, we have This proves Theorem 2.3. ¯Zs,(m,n) a.s.−−−→ n→∞ 1. 2.3.3 The asymptotic distribution of the estimator In summary, we have those three quantities ¯Zv,(m,n) = ¯Zh,(m,n) = m−1 (cid:88) i=1 n−1 (cid:88) i=1 ∆h,iVv,i/ E(Vv,i), ∆v,iVh,i/ E(Vh,i), ¯Zs,(m,n) = Vs,(m,n)/ E (cid:0)Vs,(m,n) (cid:1). (2.66) (2.67) (2.68) (2.69) Denote ¯Z(m,n) ≜ (cid:0) ¯Zv,(m,n), ¯Zh,(m,n), ¯Zs,(m,n) (cid:1)′, and Φ(m,n) ≜ cov (cid:0) ¯Z(m,n) (cid:1). Lemma 2.3. Suppose the two sampling functions φh(·), φv(·) satisfy Condition 1 and lim m→∞ n m = ρ. Then nΦ(m,n) →  2Cλ        (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) v (s)}−2ds v (s)}−1ds (cid:17)2 2ρCµ 0 0 0 0 (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) h (s)}−2ds h (s)}−1ds (cid:17)2 0 0 0         ≜ Φ0. (2.70) Proof. Elements on the main diagonal of Φ0 are implied by Theorems 2.2 and 2.3. Moreover, by the Cauchy-Schwarz inequality (cid:12) (cid:12)cov (cid:0) ¯Zv,(m,n), ¯Zs,(m,n) (cid:1)(cid:12) (cid:12) ≤ (cid:113)(cid:0)var (cid:0) ¯Zv,(m,n) (cid:1)(cid:1)(cid:0)var (cid:0) ¯Zs,(m,n) (cid:1)(cid:1) = O(n −3 2 ). So Φ0(1, 3) = 0. It remains to prove Φ0(1, 2) = 0. Let fλ(t) = e−λ|t|, then f (1) λ (t) =    −λe−λt, if t > 0, λeλt, if t < 0. 25 (2.71) The key step in our proof is to show that for any i, j ∈ I(m,n), (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) E (cid:0)Xi+(0,1) − Xi (cid:1)(cid:0)Xj+(1,0) − Xj (cid:12) (cid:12) (cid:1) (cid:12) (cid:12) (cid:12) ≤ σ2λµ∆h,j∆v,i. (2.72) Actually, (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) E (cid:0)Xi+(0,1) − Xi (cid:1)(cid:0)Xj+(1,0) − Xj (cid:12) (cid:12) (cid:1) (cid:12) (cid:12) (cid:12) =σ2 e−λ|φ(h,i[1])−φ(h,j[1])|−µ|φ(v,i[2])−φ(v,j[2])| (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) + e−λ|φ(h,i[1])−φ(h,j[1]+1)|−µ|φ(v,i[2]+1)−φ(v,j[2])| − e−λ|φ(h,i[1])−φ(h,j[1]+1)|−µ|φ(v,i[2])−φ(v,j[2])| − e−λ|φ(h,i[1])−φ(h,j[1])|−µ|φ(v,i[2]+1)−φ(v,j[2])| (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) =σ2 (e−λ|φ(h,i[1])−φ(h,j[1])| − e−λ|φ(h,i[1])−φ(h,j[1]+1)|) (cid:12) (cid:12) (e−µ|φ(v,i[2])−φ(v,j[2])| − e−µ|φ(v,i[2]+1)−φ(v,j[2])|) (cid:12) (cid:12) (cid:12) . Since (cid:0)φ(h,i[1]) − φ(h,j[1]) (cid:1)(cid:0)φ(h,i[1]) − φ(h,j[1]+1) (cid:1) ≥ 0, it follows from Eq. (2.71) and the mean value theorem that (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (e−λ|φ(h,i[1])−φ(h,j[1])| − e−λ|φ(h,i[1])−φ(h,j[1]+1)|) (cid:12) (cid:12) (cid:12) ≤ λ∆h,j. The same method shows (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (e−µ|φ(v,i[2])−φ(v,j[2])| − e−µ|φ(v,i[2]+1)−φ(v,j[2])|) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ µ∆v,i. (2.73) (2.74) Therefore, Eq. (2.72) holds. Based on the definition of ∇vXi in Eq. (2.37), we have that (cid:12) (cid:12)E (cid:0)Xi+(0,1) − Xi (cid:1)(cid:0)Xj+(1,0) − Xj (cid:1)(cid:12) (cid:12) ≤ σ2λµ, (2.75) |E (∇vXi∇hXj)| = ∆h,j∆v,i 26 which implies that for any (i, j)′ ∈ I(m,n), cov (Vv,i, Vh,j) = = n−1 (cid:88) m−1 (cid:88) k=1 n−1 (cid:88) s=1 m−1 (cid:88) k=1 s=1 cov (cid:0)(∇vX(i,k))2, (∇hX(s,j))2(cid:1) 2(cid:0)E ∇vX(i,k)∇hX(s,j) (cid:1)2 (2.76) Finally, ≤ 2mn(σ2λµ)2. ΦN (1, 2) = cov (cid:32)m−1 (cid:88) i=1 ∆h,iVv,i/ E(Vv,i), (cid:33) ∆v,jVh,j/ E(Vh,j) n−1 (cid:88) j=1 m−1 (cid:88) n−1 (cid:88) = i=1 j=1 ∆h,i∆v,j cov(Vv,i, Vh,j)/ E(Vv,i)/ E(Vh,j) (2.77) ≤ 2mn(σ2λµ)2/ E(Vv,1)/ E(Vh,1) = O((mn)−1). This proves Φ0(1, 2) = 0. Theorem 2.4. √ n( ¯Z(m,n) − 1) → N (0, Φ0). (2.78) Proof. The following argument imitates the method of Chan and Wood. Let L = (m−1)(n− 1). Fix an 3-vector f ∈ R3 + and define the 3L×3L diagonal matrix FL by FL = diag{f ′, · · · , f ′} for L ≥ 1. Define the 3L-vector WL = (Y′ L(j), j ∈ I(m,n))′, where (cid:32) YL(j) = (cid:112)∆h,j ∇vXj (cid:112)E(Vv,1) , (cid:112)∆v,j ∇hXj (cid:112)E(Vh,1) , ∇sXj (cid:112)E(Vs,(m,n)) (cid:33)′ . By construction, we have √ SL ≜ nf ′(cid:0) ¯Z(m,n) − 1(cid:1) LFLWL − E(W′ Let VL denote the covariance matrix of WL. Note that each entry of VL is of the form LFLWL)). n(W′ = √ (2.79) L (i, j) with a, b ∈ {v, h, s}. For example, σab (cid:32) σvh L (i, j) = cov (cid:112)∆h,i ∇vXi (cid:112)E(Vv,1) , (cid:112)∆v,j ∇hXj (cid:112)E(Vh,1) (cid:33) . (2.80) 27 For convenience, σab L (·) will be denoted as σab(·) below. V 1 2 L is defined as the symmetric positive definite square root of VL. Denote by ΛL = diag(λ1,L, · · · , λ3L,L) the diagonal matrix whose diagonal entries are eigenvalues of 2n 1 2 V 1 2 1 2 L FLV L . Then for ϵL ∼ N (0, I3L×3L), we have n 1 2 W′ LFLWL d= n 1 2 ϵ′ LV L FLV 1 2 1 2 L ϵL d= 1 2 ϵ′ LΛLϵL. (2.81) Therefore, for all |θ| < max(λ1,L,··· ,λ3L,L), the cumulant generating function of SL is given by 1 κL(θ) ≜ ln(cid:0)E(eθSL)(cid:1) = − 1 2 3L (cid:88) q=1 {ln(1 − θλq,L) + θλq,L}, (2.82) (see Khuri 2009, chap. 5). To obtain the limit of κL(θ) as n → ∞, we first prove tr(cid:0)Λ4 L (cid:1) = 3L (cid:88) q=1 λ4 q,L → 0, as n → ∞. (2.83) Direct calculation shows that tr(cid:0)Λ4 L (cid:1) = 16n2 tr (cid:16) (V 1 2 L FLV 1 2 L )4(cid:17) = 16n2 tr (cid:0)(VLFL)4(cid:1) = 16n2 (cid:88) · · · (cid:88) u1∈{v,h,s} u4∈{v,h,s} fu1 · · · fu4∆L(u1, · · · , u4), where ∆L(u1, · · · , u4) (cid:88) = · · · (cid:88) i1∈I(m,n) i4∈I(m,n) σu1u2(i1, i2)σu2u3(i2, i3)σu3u4(i3, i4)σu4u1(i1, i4). (2.84) (2.85) For a, b ∈ {v, h, s} and i, j ∈ I(m,n), if we are able to find an upper bound ˆσab(i − j) for each (cid:12)σab(i, j)(cid:12) (cid:12) 2000, (7.15)). (cid:12), the proof can be completed by imitating the stationary case; (see Chan and Wood 28 • For the vertical case, we have (cid:112)∆h,i cov |σvv(i, j)| (cid:12) (cid:32) (cid:12) (cid:12) = (cid:12) (cid:12) (cid:112)∆h,i∆h,j E(Vv,1) = ≤ Ch,1 E(Vv,1)m ≤ Ch,1 E(Vv,1)m ∇vXi (cid:112)E(Vv,1) , (cid:112)∆h,j ∇vXj (cid:112)E(Vv,1) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) |cov (∇vXi, ∇vXj)| e−λ|φ(h,i[1])−φ(h,j[1])| (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) s,k=0 (cid:12) (cid:12) (cid:12) av;i,kav;j,sΓµ(φ(v,(i[2]+k)) − φ(v,(j[2]+s))) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) s,k=0 (cid:12) (cid:12) (cid:12) av;i,kav;j,sΓµ(φ(v,(i[2]+k)) − φ(v,(j[2]+s))) (cid:12) (cid:12) . Then imitating the proof of Lemma 2.1, we get that there exist some constants C1, C2 > 0 such that |σvv(i, j)| ≤   C1 mn ,  C2 mn2 , if i[2] − j[2] = 0 , if i[2] − j[2] ̸= 0 , (2.86) uniformly over i, j ∈ I(m,n) when n is big enough. Define ˆσvv(x) = C1 mn 10(x[2]) + C2 mn2 1R\0(x[2]), then when n is big enough. |σvv(i, j)| ≤ ˆσvv(i − j), • For the square case, noticing Eq. (2.56), we have that there exist some constants C3, C4, C5, C6 > 0, such that |σss(i, j)| =|E (∇sXi∇sXj)|/ E(Vs,(m,n)) ≤ (cid:0)C3m10((i − j)[1]) + C41R\0((i − j)[1])(cid:1) 1 (mn)2 × (cid:0)C5n10((i − j)[2]) + C61R\0((i − j)[2])(cid:1), (2.87) 29 uniformly over i, j ∈ I(m,n) when n is big enough. So we define ˆσss(·) as ˆσss(x) = 1 (mn)2 (cid:0)C3m10(x[1]) + C41R\0(x[1])(cid:1) × (cid:0)C5n10(x[2]) + C61R\0(x[2])(cid:1). • For the cross term between the vertical and the horizontal, (cid:12)σvh(i, j)(cid:12) (cid:12) (cid:12) (cid:12) (cid:32) (cid:12) (cid:12) = (cid:12) (cid:12) cov (cid:112)∆h,i ∇vXi (cid:112)E(Vv,1) , (cid:112)∆v,j ∇hXj (cid:112)E(Vh,1) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (2.88) = (cid:112)∆h,i (cid:112)∆v,j (cid:112)E(Vv,1)(cid:112)E(Vh,1) |E (∇vXi∇hXj)|. Then by Eq. (2.75), there exists some constant C7 such that (cid:12)σvh(i, j)(cid:12) (cid:12) (cid:12) ≤ C7 (mn) 3 2 , uniformly over i, j ∈ I(m,n) when n is big enough. In this case, ˆσvh(x) = C7 (mn) 3 2 is a constant function with respect to x. • Lastly, we study the cross term between the vertical and the square, ∇vXi (cid:112)E(Vv,1) , ∇sXj (cid:112)E(Vs,(m,n)) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)cov (cid:0)∇vXi, (∇vXj+(1,0) − ∇vXj)/∆h,j (cid:1)(cid:12) (cid:12) |σvs(i, j)| (cid:12) (cid:32) (cid:12) (cid:12) = (cid:12) (cid:12) cov (cid:112)∆h,i = = (cid:112)∆h,i (cid:112)E(Vv,1) E(Vs,(m,n)) (cid:112)∆h,i (cid:112)E(Vv,1) E(Vs,(m,n)) (cid:12) (cid:12) (cid:12)e−λ|φ(h,i[1])−φ(h,j[1]+1)| − e−λ|φ(h,i[1])−φ(h,j[1])|(cid:12) (cid:12) (cid:12)/∆h,j (2.89) × (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) s,k=0 (cid:12) (cid:12) (cid:12) av;i,kav;j,sΓµ(φ(v,(i[2]+k)) − φ(v,(j[2]+s))) (cid:12) (cid:12) ≤λ (cid:112)∆h,i (cid:112)E(Vv,1) E(Vs,(m,n)) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) s,k=0 av;i,kav;j,sΓµ(φ(v,(i[2]+k)) − φ(v,(j[2]+s))) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) , 30 where the last inequality is based on Eq. (2.73). Therefore, we conclude that there exist some constants C8, C9 > 0 such that |σvs(i, j)| ≤   C8 m1.5n,  C9 m1.5n2 , if i[2] − j[2] = 0 , if i[2] − j[2] ̸= 0 , (2.90) uniformly over i, j ∈ I(m,n) when n is big enough. Define ˆσvs(x) = C8 m1.5n 10(x[2]) + C9 m1.5n2 1R\0(x[2]), then |σvs(i, j)| ≤ ˆσvs(i − j), when n is big enough. So combining all the cases above, we have |∆L(u1, · · · , u4)| (cid:88) · · · (cid:88) i1∈I(m,n) i4∈I(m,n) (cid:88) · · · (cid:88) ≤ ≤ i1∈I(m,n) when n is big enough. i4∈I(m,n) Define the index set |σu1u2(i1, i2)σu2u3(i2, i3)σu3u4(i3, i4)σu4u1(i1, i4)| (2.91) ˆσu1u2(i1 − i2)ˆσu2u3(i2 − i3)ˆσu3u4(i3 − i4)ˆσu4u1(i1 − i4), D(m,n) ≜ {i − j : i, j ∈ I(m,n)}. (2.92) For each triple (h1, h2, h3) which satisfies ha ∈ D(m,n), 1 ≤ a ≤ 3, the cardinality of the set {(i1, i2, i3, i4) : ia ∈ In, a = 1, · · · , 4; ha = ia − ia+1, 1 ≤ a ≤ 3} is bounded by L. It follows that when n is big enough |∆L(u1, · · · , u4)| ≤L (cid:88) · · · (cid:88) ˆσu1u2(h1)ˆσu2u3(h2)ˆσu3u4(h3)ˆσu4u1(h1 + h2 + h3) h1∈D(m,n)  3 (cid:89) h3∈D(m,n)  ≤C10 (cid:88)  ˆσuaua+1(h)  . a=1 ha∈D(m,n) 31 (2.93) The last inequality holds since there exists some constant C10 > 0 such that ˆσab(h) ≤ C10 mn , for all possible a, b ∈ {v, h, s} and h ∈ D(m,n) when n is big enough. Direct calculation shows that (cid:88) ˆσvv(h) = O( 1 n ), h∈D(m,n) (cid:88) h∈D(m,n) ˆσvh(h) = O( √ 1 mn ), (cid:88) ˆσss(h) = O( 1 mn ), h∈D(m,n) (cid:88) h∈D(m,n) ˆσvs(h) = O( √ 1 mn ). Combining Eqs. (2.84) and (2.93), we conclude that tr(cid:0)Λ4 L (cid:1) = n2O(n−3) = O(n−1), (2.94) as n → ∞. Namely, Eq. (2.83) holds. Meanwhile, note that Eq. (2.83) implies that as n → ∞, max 1≤q≤3L {λq,L} ≤ (cid:32) (cid:88) 1≤q≤3L (cid:33) 1 4 λ4 q,L → 0. (2.95) Expanding Eq. (2.82) about θ = 0 using Taylor’s theorem, we obtain κL(θ) = 1 2 3L (cid:88) q=1 (cid:26) 1 2 (θλq,L)2 + 1 3 (θλq,L)3 + 1 4 (θλq,L)4(1 − θ∗ q,Lλq,L)−4 (cid:27) , (2.96) for some θ∗ which satisfies 0 ≤ (cid:12) (cid:12)θ∗ q,L q,L Let us first consider the term (cid:80)3L q=1 (cid:12) (cid:12) ≤ |θ|. 1 2(θλq,L)2, 3L (cid:88) q=1 1 2 λ2 q,L = 1 2 tr(cid:0)Λ2 L (cid:1) = 2n tr(cid:0)(VLFL)2(cid:1) √ = var( nW ′ LFLWL) = nf ′Φ(m,n)f . It follows from Lemma 2.3 that, as n → ∞, 3L (cid:88) q=1 1 2 q,L = nf ′Φ(m,n)f → f ′Φ0f . λ2 32 (2.97) (2.98) Second, Eq. (2.83) implies that, as n → ∞, (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 3L (cid:88) q=1 λ3 q,L (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ max 1≤q≤3L {λq,L} 3L (cid:88) q=1 λ2 q,L → 0. (2.99) Third, note that δ = sup L≥1 max 1≤q≤3L {λq,L} is positive and finite. If we restrict attention to |θ| ≤ 1 2δ , then (1 − θ∗ q,Lλq,L) > 1/2, (2.100) for all q and L. Hence, combining Eq. (2.98), Eq. (2.99) and Eq. (2.100), we have that for θ ∈ (− 1 2δ , 1 2δ ), κL(θ) → θ2f ′Φ0f /2, as n → ∞, (2.101) which leads to the conclusion that SL → N (0, f ′Φ0f ) in distribution. This completes the proof. Let Rv ≜ E(Vv,1)/ (cid:18) (cid:90) 1 2n2 Rh ≜ E(Vh,1)/ 2m2 (cid:18) (cid:19) , (cid:19) , {φ(1) v (s)}−1ds 0 (cid:90) 1 0 (cid:18) {φ(1) h (s)}−1ds (cid:90) 1 (2.102) (2.103) . (2.104) (cid:19) {φ(1) v (s)}−1ds (cid:90) 1 0 Rs ≜ E(Vs,(m,n))/ 4(mn)2 {φ(1) h (s)}−1ds 0 Based on Eqs. (2.18), (2.20) and (2.57), we have that rn ≜ (Rv, Rh, Rs) = (cid:0)σ2µ + O(n−1), σ2λ + O(n−1), σ2λµ + O(n−1)(cid:1), (2.105) as n → ∞. Define the estimator of p ≜ (cid:0)σ2µ, σ2λ, σ2λµ(cid:1)′ as ˆZ(m,n) = ( ¯Zv,(m,n)Rv, ¯Zh,(m,n)Rh, ¯Zs,(m,n)Rs)′, 33 and σ2µ        P ≜ σ2λ  .       σ2λµ Consequently, the estimator of (λ, µ, σ2) is defined as ˆλ = ˆZ(m,n)[3]/ ˆZ(m,n)[1], ˆµ = ˆZ(m,n)[3]/ ˆZ(m,n)[2], ˆσ2 = ˆZ(m,n)[1] ˆZ(m,n)[2]/ ˆZ(m,n)[3]. We end up with the corollary below: Corollary 2.1. Suppose the two sampling functions φh(·) and φv(·) satisfy Condition 1 and √   n                              −                              λ µ σ2 ˆλ ˆµ ˆσ2 0 limm→∞ n m = ρ, then where Σ is defined as (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ (1) v (s)}−2ds (1) v (s)}−1ds (cid:17)2 λ2 2Cλ      and d−−−→ n→∞ N (0,Σ), (2.106) 2ρCµ (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ (1) h (s)}−2ds (1) h (s)}−1ds (cid:17)2 µ2 (cid:32) 2σ4 ρCµ (1) v (s)}−2ds (1) v (s)}−1ds h (s)}−2ds (1) h (s)}−1ds (1) (cid:17)2 σ2λ (cid:17)2 σ2µ −2Cλ −2ρCµ (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ h (s)}−2ds (1) h (s)}−1ds (1) (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ (cid:17)2 + Cλ (cid:82) 1 0 {φ (cid:16)(cid:82) 1 0 {φ (1) v (s)}−2ds (1) v (s)}−1ds (cid:17)2      (cid:33) ,             ˆλ ˆµ ˆσ2 a.s.−−−→ n→∞       λ µ σ2       . (2.107) Proof. Let Rn ≜ Rv        Rh        , Rs 34 then we have Since we have √ n (cid:16) ˆZ(m,n) − p (cid:17) = √ n(cid:2)Rn (cid:0) ¯Z(m,n) − 1(cid:1) + (rn − p)(cid:3). (2.108) Rn n→∞−−−→ P, √ nRn (cid:0) ¯Z(m,n) − 1(cid:1) d−→ N (0, PΦ0P′), based on Theorem 2.4. Moreover, as n → ∞, √ n(rn − p) = (cid:16) Therefore, we have O(n− 1 2 ), O(n− 1 2 ), O(n− 1 2 ) (cid:17)′ . √ (cid:16) ˆZ(m,n) − p n (cid:17) d−→ N (0, PΦ0P′). (2.109) Applying the delta method, we get Eq. (2.106). The almost sure convergence is the direct consequence of Theorems 2.2 and 2.3. 2.4 Simulations Consider the 2-dimensional OU field with λ = 0.5, µ = 10, σ2 = 4, n m = 2 and two sampling functions as φh(x) = (cid:32)(cid:18) 100 102 x + (cid:19)2 − 1 100 φv(x) = (cid:32)(cid:18) 20 22 x + (cid:19)2 − 1 20 1 202 (cid:33) , 1 1002 (cid:33) . With m = 100, 150, 200, · · · , 500, we run 1000 realizations for each case and estimate λ, µ and σ2. 35 Figure 2.1 Simulated estimation for the 2-dimensional OU field . Plots in the first row present the averaged absolute value of bias for each sample size and each parameter; plots in the second row present the distribution of normalized bias when n = 1000, where the red curve is the density of N (0, 1). 2.5 Extensions to higher-dimensional spaces In this section, we generalize the results for dimension d = 2 to higher-dimensional spatial processes. Let X(u), u ∈ [0, 1]d, denote a spatial Gaussian process with the covariance function Eq. (2.4). Like the case d = 2, for ℓ = 1, · · · , d, φℓ satisfies Condition 1. We consider the N = (cid:81)d ℓ=1 nℓ observation locations over Rd which are of the form u(i) ≜ (cid:18) φ1 (cid:18) i[1] − 1 n1 − 1 (cid:19) , φ2 (cid:18) i[2] − 1 n2 − 1 (cid:19) , · · · , φd (cid:18) i[d] − 1 nd − 1 (cid:19)(cid:19)′ , where 1 ≤ i[ℓ] ≤ nℓ and N 1/d/nℓ → ρℓ ∈ (0, ∞), ℓ = 1, · · · , d, and we write Xi = X(u(i)). To simplify the notation, let φ(ℓ) j = φℓ (cid:18) j − 1 nℓ − 1 (cid:19) , for j = 1, · · · , nℓ and ℓ = 1, · · · , d. Furthermore, define ∆(ℓ) k = φ(ℓ) k[ℓ]+1 − φ(ℓ) k[ℓ], 36 for k ∈ IN ≜ {i : 1 ≤ i ≤ ((n1 − 1), · · · , (nd − 1))′} and ℓ = 1, · · · , d. If k = 1, · · · , (nℓ − 1) is a scalar, then ∆(ℓ) k = φ(ℓ) k+1 − φ(ℓ) k . To estimate σ2 (cid:81)d ℓ=1 θℓ, we first introduce the increment {bk}1 k=0 Accordingly, define bk = (−1)(d−|k|). ∇sXi = Vs,N = (cid:80)1 , k=0 bkXi+k (cid:81)d ℓ=1 ∆(ℓ) (∇sXi)2, i (cid:88) i∈IN ¯Zs,N = Vs,N / E(Vs,N ). Similarly to Eq. (2.57), it is not hard to see which is defined as (2.110) (2.111) (2.112) (2.113) E (Vs,N ) = σ2 d (cid:89) ℓ=1 (cid:18) 2θℓn2 ℓ (cid:90) 1 0 (cid:19) {φ(1) ℓ (s)}−1ds + O(N 2−1/d). (2.114) As for σ2 (cid:81) ℓ̸=j θℓ, j = 1, · · · , d, we first define {b(j) k } as follows: b(j) k = (−1)(d−1−|k|), where k = 0, · · · , (1 − ej) and {ej} is the standard basis of Rd; moreover, let ∇jXi = (cid:80)1−ej k Xi+k k=0 b(j) (cid:81) ℓ̸=j ∆(ℓ) i , Vj,m = (cid:88) (∇jXi)2, i∈Ij,m,N ¯Zj,N = nj −1 (cid:88) m=1 ∆(j) m Vj,m/ E (Vj,m) = nj −1 (cid:88) m=1 ∆(j) m Vj,m/ E (Vj,1), (2.115) (2.116) (2.117) (2.118) where Ij,m,N ≜ {i ∈ IN : i[j] = m}. Consequently, E (Vj,m) = σ2 (cid:89) ℓ̸=j (cid:18) 2θℓn2 ℓ (cid:90) 1 0 (cid:19) {φ(1) ℓ (s)}−1ds + O(N 2−3/d). (2.119) Finally, define the (d + 1)-vector ¯ZN ≜ ( ¯Z1,N , · · · , ¯Zd,N , ¯Zs,N )′, and ΦN ≜ cov (cid:0) ¯ZN study the limiting covariance of ¯ZN in the lemma below. (cid:1). We 37 Lemma 2.4. and N var (cid:0) ¯Zs,N (cid:1) → 2 d (cid:89) ℓ=1 (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 , N 1−1/dΦN →       Σd 0 0′ 0 ≜ Φ0, where Σd is a diagonal matrix with the j-th diagonal element as 2Cθj ρℓ (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 . (2.120) (2.121) (2.122) Proof. We follow the proof of Theorem 2.3 to establish the former part of Lemma 2.4. To begin with, var (Vs,N ) = var (cid:32) (cid:88) (cid:33) (∇sXi)2 i∈IN (cid:88) 2[E (∇sXi∇sXj)]2. (cid:88) = i∈IN j∈IN Imitating Eq. (2.62), we claim that PN ≜ (cid:88) 2(cid:8)E (cid:2)(∇sXi)2(cid:3)(cid:9)2 i∈IN = 2σ4 (cid:20) d (cid:89) ℓ=1 (2θℓnℓ)2nℓ (cid:21) ℓ (s)}−2ds + O(n2 ℓ ) {φ(1) . (cid:90) 1 0 (2.123) (2.124) As for the cross terms, without loss of generality, we assume i and j are the same only in terms of the first k coordinates where k = 0, · · · , (d − 1). Imitating Eq. (2.63), we have QN,k ≜ (cid:88) 2{E [∇sXi∇sXj]}2 i,j (cid:32) = O 2σ4 k (cid:89) ℓ=1 (cid:26) (2θℓnℓ)2nℓ (cid:90) 1 0 {φ(1) ℓ (s)}−2ds (cid:27) d (cid:89) (cid:33) . n2 j j=k+1 (2.125) By Eqs. (2.124) and (2.125), we see that var (Vs,N ) = 2σ4 d (cid:89) (cid:20) (2θℓnℓ)2nℓ (cid:90) 1 ℓ=1 0 38 (cid:21) ℓ (s)}−2ds {φ(1) + O(N 3−1/d). Combining Eq. (2.114) with the above equation, we establish Eq. (2.120). Next, we deal with Σd. Based on Theorem 2.2, we claim that nℓ var (cid:0) ¯Zj,N (cid:1) → 2Cθj (cid:89) ℓ̸=j (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 , (2.126) is a coefficient defined in Eq. (2.49). Moreover, since N 1/d/nℓ → ρℓ, the diagonal where Cθj elements of Σd have the form Eq. (2.122). In order to show that the off-diagonal elements of Σd are 0, for example, (Σd)1,2 = 0, we define {b(12) k } as follows: b(12) k = (−1)(d−2−|k|), (2.127) where k = 0, · · · , (1 − e1 − e2). Correspondingly, define ∇12Xi = (cid:80) (cid:81) k b(12) k Xi+k ℓ /∈{1,2} ∆(ℓ) i . (2.128) Similar to the derivation of Eq. (2.56), we claim that E (∇12Xi∇12Xj) =σ2e−θ1 (cid:12) (cid:12)φ(1) (cid:12) i[1]−φ(1) j[1] (cid:12) (cid:12) (cid:12)−θ2 (cid:12) (cid:12)φ(1) (cid:12) i[2]−φ(1) j[2] d (cid:89) (cid:12) (cid:12) (cid:12) (cid:32) 1 (cid:88) ℓ=3 l,k=0 i,ka(ℓ) a(ℓ) j,l e−θℓ (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:33) (cid:12) (cid:12) (cid:12) , (2.129) where a(ℓ) i,0 = −1/∆(ℓ) i and a(ℓ) i,1 = 1/∆(ℓ) i . Meanwhile, notice that ∇1Xi = (∇12Xi+e2 − ∇12Xi)/∆(2) i , ∇2Xj = (∇12Xj+e1 − ∇12Xj)/∆(1) j , same as the deduction of Eq. (2.75), we get |E (∇1Xi∇2Xj)| ≤ σ2θ1θ2 d (cid:89) ℓ=3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) l,k=0 i,ka(ℓ) a(ℓ) j,l e−θℓ 39 (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) . (2.130) And cov (V1,m, V2,n) (cid:88) (cid:88) 2|E (∇1Xi∇2Xj)|2 i∈I1,m,N j∈I2,n,N (cid:88) (cid:88) 2(cid:0)σ2θ1θ2 (cid:1)2 d (cid:89) (cid:32) 1 (cid:88) i,ka(ℓ) a(ℓ) j,l e−θℓ i∈I1,m,N j∈I2,n,N ℓ=3 l,k=0 = ≤ =n1n2 d (cid:89) ℓ=3 O(n3 ℓ ); (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:33)2 (cid:12) (cid:12) (cid:12) (2.131) see Lemma 2.1 for the derivation of O(n3 ℓ ) in the last equation. By Eq. (2.119) and the definitions of ¯Z1,N as well as ¯Z2,N , we see that (cid:32) (ΦN )1,2 = O n1n2 ℓ̸=1 n2 ℓ (cid:81)d ℓ=3 n3 ℓ (cid:81) ℓ̸=2 n2 ℓ (cid:81) (cid:33) = O(N −1); therefore, (Σd)1,2 = 0. Similarly, we can prove (Σd)i,j = 0, where i ̸= j and i, j = 1, · · · , d. Lastly, by the Cauchy-Schwarz inequality, we have that (ΦN )d+1,i = O(N 1/(2d)−1), (2.132) which leads to (Φ0)d+1,i = 0 for i = 1, · · · , d. Therefore, we have proven Lemma 2.4. Theorem 2.5. √ N 1−1/d( ¯ZN − 1) → N (0, Φ0). (2.133) Proof. The following is just a reiteration of the proof of Theorem 2.4 for higher-dimensional spaces. Let L = (cid:81)d ℓ=1(nℓ−1). Fix an (d+1)-vector f ∈ R(d+1) and define the (d+1)L×(d+1)L + diagonal matrix FL by FL = diag{f ′, · · · , f ′} for L ≥ 1. Define the (d + 1)L-vector WL = (Y′ L(j), j ∈ IN )′, where YL(j) = (cid:32)(cid:113) ∆(1) j ∇1Xj (cid:112)E(V1,1) , · · · , (cid:113) ∆(d) j ∇dXj (cid:112)E(Vd,1) , ∇sXj (cid:112)E(Vs,N ) (cid:33)′ . By construction, we have SL ≜ √ N 1−1/df ′( ¯ZN − 1) √ = N 1−1/d(W′ LFLWL − E(W′ LFLWL)). 40 (2.134) Let VL denote the covariance matrix of WL. Note that each entry of VL is of the form L (i, j) with a, b = 1, · · · , (d + 1). For example, σab and σ12 L (i, j) = cov (cid:32)(cid:113) ∆(1) j ∇1Xj (cid:112)E(V1,1) , (cid:113) ∆(2) j ∇2Xj (cid:112)E(V2,1) (cid:33) , σ1(d+1) L (i, j) = cov (cid:32)(cid:113) ∆(1) j ∇1Xj (cid:112)E(V1,1) , ∇sXj (cid:112)E(Vs,N ) (cid:33) . (2.135) (2.136) For convenience, σab L (·) will be denoted as σab(·) below. V 1 2 L is defined as the symmetric positive definite square root of VL. Denote by ΛL = diag(λ1,L, · · · , λ(d+1)L,L), the diagonal matrix whose diagonal entries are eigenvalues of 2 √ N 1−1/dV L FLV 1 2 1 2 L . Then for ϵL ∼ N (0, I(d+1)L×(d+1)L), we have √ N 1−1/dW′ LFLWL √ d= N 1−1/dϵ′ LV L FLV 1 2 1 2 L ϵL d= 1 2 ϵ′ LΛLϵL. (2.137) Therefore, for all |θ| < max(λ1,L,··· ,λ(d+1)L,L) , the cumulant generating function of SL is given 1 by κL(θ) ≜ ln(cid:0)E(eθSL)(cid:1) = − 1 2 (d+1)L (cid:88) {ln(1 − θλq,L) + θλq,L}, q=1 (2.138) (see Khuri 2009, chap. 5). To obtain the limit of κL(θ) as N → ∞, we first prove tr(cid:0)Λ4 L (cid:1) = (d+1)L (cid:88) q=1 λ4 q,L → 0, as N → ∞. (2.139) Direct calculation shows that tr(cid:0)Λ4 L (cid:1) = 16N 2−2/d tr (cid:16) 1 2 (V L FLV 1 2 L )4(cid:17) = 16N 2−2/d tr (cid:0)(VLFL)4(cid:1) (2.140) = 16N 2−2/d d+1 (cid:88) · · · d+1 (cid:88) u1=1 u4=1 fu1 · · · fu4∆L(u1, · · · , u4), 41 where ∆L(u1, · · · , u4) (cid:88) = · · · (cid:88) i1∈IN i4∈IN σu1u2(i1, i2)σu2u3(i2, i3)σu3u4(i3, i4)σu4u1(i1, i4). (2.141) For a, b = 1, · · · , (d + 1) and i, j ∈ IN , if we are able to find an upper bound ˆσab(i − j) for each (cid:12) (cid:12)σab(i, j)(cid:12) (cid:12), the proof can be completed by imitating the stationary case; (see Chan and Wood 2000, (7.15)). Without loss of generality, we only consider a, b ∈ {1, 2, (d + 1)} : • For the 1-1 case, we have (cid:12)σ11(i, j)(cid:12) (cid:12) (cid:12) (cid:12) (cid:32)(cid:113) (cid:12) (cid:12) = (cid:12) (cid:12) (cid:113) cov ∆(1) i ∆(1) = ≤ j E(V1,1) σ2C1,1 E(V1,1)n1 ∆(1) i ∇1Xi (cid:112)E(V1,1) , (cid:113) ∆(1) j ∇1Xj (cid:112)E(V1,1) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) |cov (∇1Xi, ∇1Xj)| (cid:12) (cid:12)φ(1) (cid:12) i[1]−φ(1) j[1] e−θ1 d (cid:89) (cid:12) (cid:12) (cid:12) ℓ=2 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) l,k=0 i,ka(ℓ) a(ℓ) j,l e−θℓ (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ σ2C1,1 E(V1,1)n1 d (cid:89) ℓ=2 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) l,k=0 i,ka(ℓ) a(ℓ) j,l e−θℓ (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) , where C1,1 is defined based on Eq. (2.13). Define ˆσ11(x) = C5 (cid:81)d ℓ=2 n2 ℓ n1 d (cid:89) (cid:16) ℓ=2 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , where C5, C (ℓ) 1 , and C (ℓ) 2 are positive constants. we claim that (cid:12)σ11(i, j)(cid:12) (cid:12) (cid:12) ≤ ˆσ11(i − j), uniformly over i, j ∈ IN when N is big enough. • For the square case, we have that |σss(i, j)| =|E (∇sXi∇sXj)|/ E(Vs,N ) (2.142) ≤ C6 N 2 d (cid:89) (cid:16) ℓ=1 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , 42 uniformly over i, j ∈ IN when N is big enough. So we define ˆσss(·) as ˆσss(x) = C6 N 2 d (cid:89) (cid:16) ℓ=1 where C6 is a positive constant. • For the 1-2 cross term, 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , (cid:12)σ12(i, j)(cid:12) (cid:12) (cid:12) (cid:12) (cid:32)(cid:113) (cid:12) (cid:12) = (cid:12) (cid:12) cov ∆(1) i ∇1Xi (cid:112)E(V1,1) , (cid:113) ∆(2) j ∇2Xj (cid:112)E(V2,1) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (2.143) (cid:113) ∆(1) i (cid:113) ∆(2) j (cid:112)E(V1,1)(cid:112)E(V2,1) According to Eq. (2.130), there exists some constant C7 such that |E (∇1Xi∇2Xj)|. = (cid:12)σ12(i, j)(cid:12) (cid:12) (cid:12) ≤ C7 2 (cid:81) (n1n2) 3 d (cid:89) (cid:16) ℓ=3 n2 ℓ ℓ=3 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , uniformly over i, j ∈ IN when n is big enough. In this case, ˆσ12(x) = C7 2 (cid:81)d (n1n2) 3 ℓ=3 n2 ℓ d (cid:89) (cid:16) ℓ=3 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , where C7 is a positive constant. • Lastly, we study the 1-s cross term , (cid:12)σ1s(i, j)(cid:12) (cid:12) (cid:12) (cid:12) (cid:32)(cid:113) (cid:12) (cid:12) = (cid:12) (cid:12) cov ∆(1) i ∇1Xi (cid:112)E(V1,1) , ∇sXj (cid:112)E(Vs,N ) (cid:33)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:113) ∆(1) i (cid:112)E(V1,1) E(Vs,N ) (cid:113) σ2 ∆(1) i d (cid:89) (cid:112)E(V1,1) E(Vs,N ) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:113) 1 (cid:88) l,k=0 ℓ=2 × θ1σ2 ∆(1) i (cid:112)E(V1,1) E(Vs,N ) = = ≤ (cid:16) (cid:12) (cid:12) (cid:12)cov ∇1Xi, (∇1Xj+(1,0) − ∇1Xj)/∆(1) j (cid:17)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12)φ(1) (cid:12) i[1]−φ(1) j[1]+1 (cid:12) (cid:12) (cid:12) (cid:12) e−θ1 (cid:12) (cid:12) (cid:12) − e−θ1 (cid:12) (cid:12)φ(1) (cid:12) i[1]−φ(1) j[1] (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) /∆(1) j (2.144) i,ka(ℓ) a(ℓ) j,l e−θℓ (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) d (cid:89) ℓ=2 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 1 (cid:88) l,k=0 i,ka(ℓ) a(ℓ) j,l e−θℓ 43 (cid:12) (cid:12)φ(ℓ) (cid:12) (i[ℓ]+k)−φ(ℓ) (j[ℓ]+l) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) , where the last inequality is based on Eq. (2.73). Define ˆσ1s(x) = C8 2 (cid:81)d (n1) 3 ℓ=2 n2 ℓ d (cid:89) (cid:16) ℓ=2 1 nℓ10(x[ℓ]) + C (ℓ) C (ℓ) 2 1R\0(x[ℓ]) (cid:17) , where C8 is a positive constant, then (cid:12)σ1s(i, j)(cid:12) (cid:12) (cid:12) ≤ ˆσ1s(i − j), when N is big enough. So combining all the cases above, we have |∆L(u1, · · · , u4)| |σu1u2(i1, i2)σu2u3(i2, i3)σu3u4(i3, i4)σu4u1(i1, i4)| (2.145) ˆσu1u2(i1 − i2)ˆσu2u3(i2 − i3)ˆσu3u4(i3 − i4)ˆσu4u1(i1 − i4), ≤ ≤ (cid:88) · · · (cid:88) i1∈IN (cid:88) i4∈IN (cid:88) · · · i1∈IN i4∈IN when N is big enough. Define the index set DN ≜ {i − j : i, j ∈ IN }. (2.146) For each triple (h1, h2, h3) which satisfies ha ∈ DN , 1 ≤ a ≤ 3, the cardinality of the set {(i1, i2, i3, i4) : ia ∈ IN , a = 1, · · · , 4; ha = ia − ia+1, 1 ≤ a ≤ 3} is bounded by L. It follows that when N is big enough |∆L(u1, · · · , u4)| ≤L (cid:88) · · · (cid:88) ˆσu1u2(h1)ˆσu2u3(h2)ˆσu3u4(h3)ˆσu4u1(h1 + h2 + h3) h1∈DN (cid:32) 3 (cid:89) ≤C9 h3∈DN (cid:88) (cid:33) ˆσuaua+1(h) . a=1 ha∈DN The last inequality holds since there exists some constant C9 > 0 such that (2.147) ˆσab(h) ≤ C9 N , 44 for all possible a, b ∈ {1, 2, (d + 1)} and h ∈ DN when N is big enough. Since the integrand ˆσab(h) is separable with respect to the counting measure, direct calculation shows that ˆσ11(h) = O ˆσ12(h) = O (cid:88) h∈DN (cid:88) h∈DN (cid:32) (cid:32) (cid:33) , 1 ℓ=2 nℓ (cid:81)d √ n1n2 1 (cid:81)d ℓ=3 nℓ (cid:33) , (cid:88) h∈DN (cid:88) h∈DN ˆσss(h) = O (cid:19) , (cid:18) 1 N ˆσ1s(h) = O (cid:32) 1 (cid:81)d ℓ=2 nℓ √ n1 (cid:33) . Combining Eqs. (2.140) and (2.147), we conclude that tr(cid:0)Λ4 L (cid:1) = N 2(1−1/d)O(N −3(1−1/d)) = O(N 1/d−1), (2.148) as N → ∞. Namely, Eq. (2.139) holds. Meanwhile, note that Eq. (2.139) implies that as N → ∞, max 1≤q≤(d+1)L {λq,L} ≤    1 4 (cid:88) λ4 q,L  → 0. 1≤q≤(d+1)L (2.149) Expanding Eq. (2.138) about θ = 0 using Taylor’s theorem, we obtain κL(θ) = 1 2 (d+1)L (cid:88) { q=1 1 2 (θλq,L)2 + 1 3 (θλq,L)3 + 1 4 (θλq,L)4(1 − θ∗ q,Lλq,L)−4}, (2.150) for some θ∗ q,L which satisfies 0 ≤ |θ∗ q,L| ≤ |θ|. Let us first consider the term (cid:80)(d+1)L q=1 1 2(θλq,L)2 : (d+1)L (cid:88) q=1 1 2 √ λ2 q,L = 1 2 tr(cid:0)Λ2 L (cid:1) = 2N 1−1/d tr(cid:0)(VLFL)2(cid:1) (2.151) = var( N 1−1/dW ′ LFLWL) = N 1−1/df ′ΦN f . It follows from Lemma 2.4 that, as N → ∞, (d+1)L (cid:88) q=1 1 2 q,L = N 1−1/df ′ΦN f → f ′Φ0f . λ2 (2.152) 45 Second, Eq. (2.139) implies that, as N → ∞, (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (d+1)L (cid:88) q=1 λ3 q,L (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ≤ max 1≤q≤(d+1)L {λq,L} (d+1)L (cid:88) q=1 λ2 q,L → 0. (2.153) Third, note that δ = supL≥1 max1≤q≤(d+1)L{λq,L} is positive and finite. If we restrict attention to |θ| ≤ 1 2δ , then (1 − θ∗ q,Lλq,L) > 1/2, (2.154) for all q and L. Hence, combining Eq. (2.152), Eq. (2.153) and Eq. (2.154), we have that for θ ∈ (− 1 2δ , 1 2δ ), κL(θ) → θ2f ′Φ0f /2, as n → ∞, (2.155) which leads to the conclusion that SL → N (0, f ′Φ0f ) in distribution. This completes the proof. In this final part, we propose the estimator of (θ1, · · · , θd, σ2)′ which is ˆθ1 ≜ ... ˆθd ≜ ˆσ2 ≜ (cid:81) ℓ̸=1 (cid:81)d ℓ=1 (cid:81) ℓ̸=d (cid:81)d ℓ=1 (cid:16) 2n2 ℓ (cid:16) 2n2 ℓ (cid:16) 2n2 ℓ (cid:16) 2n2 ℓ (cid:16)(cid:81)d ℓ=1 (cid:16) (cid:16)(cid:81) (cid:81)d ℓ=1 (cid:82) 1 0 {φ(1) (cid:82) 1 0 {φ(1) (cid:17) ℓ (s)}−1ds ℓ (s)}−1ds E (Vs,N ) (cid:17) E (V1,1) ¯Zs,N ¯Z1,N , (cid:17) ℓ (s)}−1ds (cid:17) ℓ (s)}−1ds (cid:82) 1 0 {φ(1) (cid:82) 1 0 {φ(1) (cid:82) 1 0 {φ(1) ℓ (s)}−1ds (cid:82) 1 0 {φ(1) 2n2 j 2n2 ℓ (cid:16) j̸=ℓ E (Vs,N ) ¯Zs,N ¯Zd,N E (Vd,1) (cid:17)(cid:17)(d−1) (cid:81)d ℓ=1 E (Vℓ,1) (E (Vs,N ))(d−1) (cid:17)(cid:17) j (s)}−1ds In a matrix form, Eq. (2.156) can be represented as                   ˆθ1 ... ˆθd ˆσ2 k1          ≜ . . .                   k(d+1) kd 46 r1 ... rd r(d+1) , (2.156) ℓ=1 (cid:81)d (cid:0) ¯Zs,N ¯Zℓ,N (cid:1)(d−1) .          , (2.157) where (cid:0)r1, · · · , rd, r(d+1) (cid:81)d (cid:0) ¯Zs,N (cid:1) is the coefficients of (cid:0)r1, · · · , rd, r(d+1) ¯Zs,N ¯Zd,N (cid:32) ¯Zs,N ¯Z1,N , · · · , (cid:1) ≜ , ℓ=1 ¯Zℓ,N (cid:1)(d−1) (cid:33) , (2.158) (cid:1) in Eq. (2.156). Based on and (cid:0)k1, · · · , kd, k(d+1) Eqs. (2.114) and (2.119), it is easy to see that kℓ = θℓ + O(N −1/d), k(d+1) = σ2 + O(N −1/d), (2.159) where ℓ = 1, · · · , d. In the following corollary, We state all the results in high-dimensional spaces : Corollary 2.2.           a.s.−−→          θ1 ... θd σ2          , (2.160) ˆθ1 ... ˆθd ˆσ2         and   N (d−1)/(2d)                                   −                   k1 ... kd k(d+1)   d−−−→ N →∞ N  0,         , ˜Σd b b′ c (2.161) ˆθ1 ... ˆθd ˆσ2 where ˜Σd is a diagonal matrix with the diagonal elements as (cid:17) (cid:16) ˜Σd j,j = 2Cθj ρℓ (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 θ2 j , and bj = −2Cθj ρℓ (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 θjσ2,  c = 2   d (cid:88) (cid:89) ρℓ Cθj j=1 ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds  σ4.  (cid:17)2 47 Proof. The almost sure convergence is implied by Eq. (2.159) and Lemma 2.4 and the Hanson-Wright inequality; see, e.g., Theorem 2.1. Applying the delta method to Theorem 2.5, we get that N (d−1)/(2d)                            −                    1  ...      1   1 r1 ... rd r(d+1)   d−−−→ N →∞ N  0,         , Γd e e′ f (2.162) where Γd is a diagonal matrix with the diagonal elements as (Γd)j,j = 2Cθj ρℓ (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 , and moreover,   ej = −2Cθj ρℓ (cid:89) ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds (cid:17)2 ,  f = 2   d (cid:88) (cid:89) ρℓ Cθj j=1 ℓ̸=j (cid:82) 1 0 {φ(1) (cid:16)(cid:82) 1 0 {φ(1) ℓ (s)}−2ds ℓ (s)}−1ds   ; (cid:17)2                   − k1 ... kd k(d+1)                   k1          = . . . kd                                     −                    1  ...      1   1 . r1 ... rd r(d+1) k(d+1) (2.163) ˆθ1 ... ˆθd ˆσ2                 Direct calculation shows the form of the right hand side of Eq. (2.161). This proves Corol- lary 2.2. 48 CHAPTER 3 THE MULTIVARIATE ORNSTEIN-UHLENBECK PROCESS 3.1 Introduction In this section, we review the definition and some basic properties of the multivariate Ornstein-Uhlenbeck (OU) process; (see Gardiner 1985, for a fuller exposition). The mul- tivariate Ornstein-Uhlenbeck (OU) process Z(t) ∈ Rn, t ∈ R is defined by the stochastic differential equation (SDE) dZ(t) + AZ(t)dt = BdW (t), (3.1) where A ∈ Rn×n is a matrix having eigenvalues with strictly positive real part, B ∈ Rn×n, and W (t) is a n-dimensional Wiener process. Note that eAtdZ(t) + AeAtZ(t)dt = d(cid:0)eAtZ(t)(cid:1), and lim t→−∞ eAtZ(t) = 0, so we have that is, eAtZ(t) = (cid:90) t −∞ eAuBdW (u), Z(t) = (cid:90) t −∞ e−A(t−u)BdW (u). Based on the form of Z(t), we see that E(Z(t)) = 0, cov(Z(t), Z(s)) = (cid:90) min(t,s) −∞ e−A(t−u)BB′e−A′(s−u)du. Let us define the stationary covariance matrix Σ as Σ ≜ cov(Z(t), Z(t)). 49 (3.2) (3.3) (3.4) (3.5) (3.6) (3.7) Then AΣ + ΣA′ = (cid:90) t Ae−A(t−u)BB′e−A′(t−u)du −∞ (cid:90) t + e−A(t−u)BB′e−A′(t−u)A′du (cid:90) t = −∞ d du −∞ e−A(t−u)BB′e−A′(t−u)du From Eq. (3.6), we see that if t > s, = BB′. cov(Z(t), Z(s)) = e−A(t−s) (cid:90) s −∞ e−A(s−u)BB′e−A′(s−u)du = e−A(t−s)Σ, t > s, (3.8) (3.9) and similarly, cov(Z(t), Z(s)) = Σe−A′(s−t), t < s. (3.10) Velandia et al. (2017) considered the parameter estimation of the multivariate OU under fixed-domain asymptotics when  A =   θ 0 0 θ   , BB′ = 2θ     σ2 1 σ1σ2ρ σ1σ2ρ σ2 2    . Combining Eqs. (3.8) to (3.10), we see that cov(Z(t), Z(s)) = e−θ|t−s|    σ2 1 σ1σ2ρ σ1σ2ρ σ2 2    , (3.11) (3.12) where σ2 1, σ2 2 are marginal variance parameters, θ > 0 is called a correlation decay parame- ter, and the quantity ρ with |ρ| < 1 is called the colocated correlation parameter (Gneiting, Kleiber, and Schlather 2010). For simplicity, let ψ = (σ2 2, ρ, θ)′. They showed the maxi- 2, ρ)′ is strongly consistent and asymptotically normal given some prior information of ψ in Velandia et al. (2017). However, the MLE in mum likelihood estimator (MLE) of (θσ2 1, θσ2 1, σ2 this case must be found numerically. Therefore, we plan to construct an explicit estimator of (θσ2 1, θσ2 2, ρ)′, which is computationally efficient in practice and asymptotically identical 50 to the MLE under fixed-domain asymptotics. The asymptotic results are formally stated in Theorem 3.3. 3.2 Asymptotic theory for estimation Under fixed-domain asymptotics, the random process Z(t) ≜ (Z1(t), Z2(t))′ is sampled increasingly densely over a fixed and bounded set T. Without loss of generality, we consider T = [0, 1], and the sampling design Tn = {t1, · · · , tn} with 0 ≤ t1 < t2 < · · · < tn ≤ 1. The corresponding observation is denoted as Zn = (Z ′ 1,n, Z ′ 2,n)′, where Zi,n = (Zi(t1), · · · , Zi(tn))′ for i = 1, 2. Let Σ(ψ) =    σ2 1 σ1σ2ρ σ1σ2ρ σ2 2    ⊗ R(θ), where ⊗ is the Kronecker product, and R(θ) = (cid:0)e−θ|ti−tj |(cid:1) 1≤i,j≤n (3.13) (3.14) (3.15) is a matrix-valued function. By Eq. (3.12), we see that Zn ∼ N (0, Σ(ψ)). The associated likelihood function is given by Ln(ψ) = (2π)−n|Σ(ψ)|−1/2 exp − (cid:26) Z ′ nΣ(ψ)−1Zn (cid:27) . 1 2 (3.16) In the rest of this section, we will denote by θ0, σ2 01, σ2 02 and ρ0 the true but unknown pa- rameter. The following theorem defines the asymptotic behavior of the MLE with respect to Ln(ψ). Theorem 3.1 (Velandia et al. (2017)). Let Tn be dense in [0, 1] as n goes to infinity. Let J = (aθ, bθ) × (aσ1, bσ1) × (aσ2, bσ2) × (aρ, bρ), with 0 < aθ ≤ θ0 ≤ bθ < ∞, 0 < aσ1 ≤ σ2 bσ1 < ∞, 0 < aσ2 ≤ σ2 02 ≤ bσ2 < ∞, and −1 < aρ ≤ ρ0 ≤ bρ < 1. Define ˆψ = (ˆσ2 1, ˆσ2 01 ≤ 2, ˆρ, ˆθ)′ as ˆψ = argmax Ln(ψ). (3.17) ψ∈J 51 Then and √ n                               ˆθˆσ2 1 ˆθˆσ2 2 ˆρ                                 θ0σ2 01 θ0σ2 02 ρ0                d →N               −                              0, ˆθˆσ2 1 ˆθˆσ2 2 a.s.−−−→ n→∞ a.s.−−−→ n→∞ ˆρ a.s.−−−→ n→∞ θ0σ2 01, θ0σ2 02, ρ0, 2θ2 0σ4 01 2θ2 0σ2 01σ2 02ρ2 0 θ0ρ0σ2 01(1 − ρ2 0) 2θ2 0σ2 01σ2 02ρ2 0 2θ2 0σ4 02 θ0ρ0σ2 02(1 − ρ2 0) θ0ρ0σ2 01(1 − ρ2 0) θ0ρ0σ2 02(1 − ρ2 0) (1 − ρ2 0)2 (3.18) (3.19) (3.20) (3.21)                               . Our estimator for (θ0σ2 02, ρ0)′ is inspired by the results in Ying (1991) where the parameter estimation for univariate OU under fixed-domain asymptotics were studied. With 01, θ0σ2 the observation Z1,n = (Z1(t1), · · · , Z1(tn))′, we have the likelihood function L1,n(σ2 1, θ) = (2π)−n/2(cid:12) (cid:12)σ2 1R(θ)(cid:12) (cid:12) −1/2 exp (cid:26) − Z ′ 1,nR(θ)−1Z1,n (cid:27) . 1 2σ2 1 (3.22) The following result can be found in Theorem 1 of Ying (1991). Theorem 3.2 (Ying (1991)). For any fixed θf > 0, let Then ˆσ2 1 = argmax σ2 1∈(0,∞) L1,n(σ2 1, θf ). ˆσ2 1θf a.s.−−−→ n→∞ σ2 01θ0, under the sampling design Eq. (3.13). Note that for the fixed θf > 0, ˆσ2 1 has an explicit form, that is, ˆσ2 1 = 1 n Z ′ 1,nR−1(θf )Z1,n. Therefore, we define the estimator of σ2 01θ0 as 1θ ≜ θf (cid:100)σ2 n Z ′ 1,nR−1(θf )Z1,n. 52 (3.23) (3.24) (3.25) (3.26) Similarly, the estimator of σ2 02θ0 is defined as 2θ ≜ θf (cid:100)σ2 n Z ′ 2,nR−1(θf )Z2,n. (3.27) To estimate ρ0, we notice that Z3(t) ≜ Z1(t) + Z2(t) is also a mean 0 univariate OU process with covariance function cov(Z3(t), Z3(s)) = (σ2 01 + σ2 02 + 2σ01σ02ρ0)e−θ0|s−t|. (3.28) As a result, we define ˆρ ≜ Z ′ 1,nR−1(θf )Z2,n (cid:113) Z ′ 1,nR−1(θf )Z1,n · Z ′ 2,nR−1(θf )Z2,n . (3.29) The following theorem establishes the strong consistency and asymptotic normality of our estimator for (σ2 01θ0, σ2 02θ0, ρ0)′. Theorem 3.3. Under the sampling design Eq. (3.13),       (cid:100)σ2 1θ (cid:100)σ2 2θ ˆρ       a.s.−−−→ n→∞       σ2 01θ0 σ2 02θ0 ρ0       , (3.30)    2θ2 0σ4 01 2θ2 0σ2 01σ2 02ρ2 0 2θ2 0σ2 01σ2 02ρ2 0 2θ2 0σ4 02 θ0ρ0σ2 01(1 − ρ2 0) θ0ρ0σ2 02(1 − ρ2 0) θ0ρ0σ2 01(1 − ρ2 0) θ0ρ0σ2 02(1 − ρ2 0) (1 − ρ2 0)2               (3.31) .               and √ n                                              (cid:100)σ2 1θ (cid:100)σ2 2θ ˆρ σ2 01θ0 − σ2 02θ0                ρ0                               d−→ N                0,               Proof. The strong consistency in Eq. (3.30) directly follows Theorem 3.2. In the following, we deal with the asymptotic normality part. Let σ2 03 be the true but unknown marginal variance of the OU process Z3(t) which is the sum of Z1(t) and Z1(t), that is, Meanwhile, define σ2 03 ≜ σ2 01 + σ2 02 + 2σ01σ02ρ0. Z3,n = (Z3(t1), · · · , Z3(tn))′. 53 To simplify notations, for i = 1, 2, 3, we write Zi,n as Zi for the rest of the proof. For i ∈ {1, 2, 3}, and k ∈ {1, . . . , n}, let Wi,k = zi,k − e−θ0∆kzi,k−1 0i(1 − e−2θ0∆k))1/2 , (σ2 where ∆k ≜ tk − tk−1, and zi,k ≜ Zi(tk). From the proof of Theorem 2 in Ying (1991), we see that for i = 1, 2, 3, as n → ∞, √ n (cid:18) Z ′ iR−1(θf )Zi n (cid:19) θf − σ2 0iθ0 = n−1/2σ2 0iθ0 n (cid:88) k=2 (W 2 i,k − 1) + Op(n−1/2); (3.32) (see also, Du, Zhang, and Mandrekar 2009, (B.36)). Since 2Z T 1 R−1(θf )Z2 = Z T 3 R−1(θf )Z3 − Z T 1 R−1(θf )Z1 − Z T 2 R−1(θf )Z2, there is √ n (cid:18) Z T 1 R−1(θf )Z2 n θf − σ01σ02ρ0θ0 (cid:19) √ n 2 = (cid:34)(cid:18) Z ′ 3R−1(θf )Z3 n θf − σ2 03θ0 (cid:19) − (cid:18) Z ′ 1R−1(θf )Z1 n (cid:19) θf − σ2 01θ0 − (cid:18) Z ′ 2R−1(θf )Z2 n θf − σ2 02θ0 (cid:19)(cid:35) = θ0 √ 2 n =n−1/2 n (cid:88) k=2 n (cid:88) k=2 (cid:2)σ2 03 (cid:0)W 2 3,k − 1(cid:1) − σ2 01 (cid:0)W 2 1,k − 1(cid:1) − σ2 02 ˜W3,k − ˜W1,k − ˜W2,k 2 + Op(n−1/2), (3.33) (cid:0)W 2 2,k − 1(cid:1)(cid:3) + Op(n−1/2) where ˜Wi,k = σ2 0iθ0 (cid:0)W 2 i,k − 1(cid:1) for i = 1, 2, 3. Denote (cid:32) ˜Wnk = ˜W1,k, ˜W2,k, ˜W3,k − ˜W1,k − ˜W2,k 2 (cid:33)′ , then for any fixed n, it is not hard to verify that (cid:110) ˜Wnk (cid:111)n k=2 are independent random vectors with zero mean and finite variance. Since for any n and i, Wi,k ∼ N (0, 1), we have var (cid:0)W 2 i,k − 1(cid:1) = 2, cov (cid:0)W 2 i,k − 1, W 2 j,k − 1(cid:1) = 2[cov (Wi,k, Wj,k)]2 ≜ 2c2 i,j. 54 From the definition, it is straightforward that c1,2 = ρ0. Moreover, c1,3 = cov (W1,k, W3,k) (cid:18) (cid:20) W1,k, = cov σ01 σ03 σ01 σ03 var(W1,k) + (cid:19)(cid:21) σ02 σ01 W2,k W1,k + σ02 σ03 c1,2 σ01 + σ02ρ0 (cid:112)σ2 01 + σ2 02 + 2σ01σ02ρ0 . = = Similarly, Thus, where c2,3 = σ01ρ0 + σ02 (cid:112)σ2 01 + σ2 02 + 2σ01σ02ρ0 . cov( ˜Wnk) = θ2 0       2σ4 01 2σ2 01σ2 02ρ2 0 α 2σ2 01σ2 02ρ2 0 α 2σ4 02 β β γ       , (3.34) α = (2σ2 01σ2 03c2 1,3 − 2σ4 01 − 2σ2 01σ2 02c2 1,2)/2 = 2σ3 01σ02ρ0, β = (2σ2 02σ2 03c2 2,3 − 2σ4 02 − 2σ2 01σ2 02c2 1,2)/2 = 2σ01σ3 02ρ0, γ = 1 4 (2σ4 03 + 2σ4 01 + 2σ4 02 − 4σ2 01σ2 03c2 1,3 − 4σ2 02σ2 03c2 2,3 + 4σ2 01σ2 02c2 1,2) = σ2 01σ2 02(1 + ρ2 0). By the multivariate Lindeberg-Feller CLT, as n → ∞, n−1/2 n (cid:88) k=2 ˜Wnk d−→ N (0, V ), where V = cov( ˜Wnk) = θ2 0       2σ4 01 2σ2 01σ2 02ρ2 0 2σ2 01σ2 02ρ2 0 2σ4 02 2σ3 01σ02ρ0 2σ01σ3 02ρ0 2σ3 01σ02ρ0 2σ01σ3 02ρ0 σ2 01σ2 02(1 + ρ2 0) By Eqs. (3.32), (3.33) and (3.35), we have that as n → ∞, (3.35) (3.36)       ,       θf       √ n Z ′ 1R−1(θf )Z1/n Z ′ 2R−1(θf )Z2/n σ2 01θ0 σ2 02θ0             d−→ N (0, V ). (3.37) Z ′ 1R−1(θf )Z2/n σ01σ02ρ0θ0             − 55 Define function g(x1, x2, x3) = (x1, x2, x3/ x1x2)′ : R3 → R3. Then √ (cid:16) (cid:17)′ 1θ, (cid:100)σ2 (cid:100)σ2 2θ, ˆρ = g(cid:0)θf Z ′ 1R−1(θf )Z1/n, θf Z ′ 2R−1(θf )Z2/n, θf Z ′ 1R−1(θf )Z2/n(cid:1). Using the multivariate Delta method, √ n                   (cid:100)σ2 1θ (cid:100)σ2 2θ ˆρ −       σ2 01θ0 σ2 02θ0 ρ0             d−→ N (0, (∇g)V (∇g)′), (3.38) where and (∇g)V (∇g)′ =       (∇g)i,j = ∂gi ∂xj , 2θ2 0σ4 01 2θ2 0σ2 01σ2 02ρ2 0 2θ2 0σ2 01σ2 02ρ2 0 2θ2 0σ4 02 θ0ρ0σ2 01(1 − ρ2 0) θ0ρ0σ2 02(1 − ρ2 0) .       θ0ρ0σ2 01(1 − ρ2 0) θ0ρ0σ2 02(1 − ρ2 0) (1 − ρ2 0)2 Remark 3.1. Note that the construction of our estimator is based on each component of the bivariate process Z(s), and it has the same asymptotic matrix as that in Theorem 3.1. This echos the fact that cokriging is identical to kriging under the bivariate OU model in Zhang and Cai (2015). 3.3 Simulations In this section, we investigate the finite sample performance of the estimator (cid:16) (cid:17)′ (cid:100)σ2 1θ, (cid:100)σ2 2θ, ˆρ with θf = 1 in Eqs. (3.26), (3.27) and (3.29). The simulations are conducted over an irregular grid within [0, 1] with n = 200, 250, . . . , 1000. Specifically, we have considered a regular grid with increment 1/n over [0, 1]. Then the grid points have been perturbed, adding a uniform random value on [−1/(5n), 1/(5n)]. We simulate, using the R package RandomFields developed by Schlather et al. (2017), and then we estimate with 1000 realizations from a zero mean bivariate OU process where σ2 1 = σ2 2 = 1, ρ = 0.2 and θ = 0.2, 0.5. 56 Figure 3.1 Simulated estimation for the bivariate OU process with θ = 0.2. Plots in the first row present the averaged absolute value of bias for each sample size and each parameter; plots in the second row present the distribution of normalized bias when n = 1000, where the red curve is the density of N (0, 1). Figure 3.2 The setting is the same as Fig. 3.1 except θ = 0.5. 57 CHAPTER 4 THE POWERED EXPONENTIAL FIELD 4.1 Introduction In this chapter, we consider a stationary, isotropic Gaussian random field {X(t) : t ∈ Rd} with mean 0 and covariance function C(u − v) ≜ cov {X(u), X(v)} = σ2e−θ∥u−v∥α , ∀ u, v ∈ Rd, (4.1) where α ∈ (0, 2), θ > 0 and σ > 0. The corresponding spectral density is defined as the Fourier transform of Eq. (4.1): fσ,θ(ω) = 1 (2π)d (cid:90) Rd e−iω′tC(t)dt. (4.2) Note that fσ,θ : Rd → R is actually the probability density function of a sub-Gaussian random vector (see, e.g., Samorodnitsky and Taqqu 2017, Proposition 2.5.5). By the Bochner’s theorem (see, e.g., Adler and Taylor 2007, Theorem 5.4.1), Eq. (4.1) is a covariance in Rd for any d. Another popular covariance function in spatial statistics is the Matérn covariance function CM (t) = σ2 (β∥t∥)ν 2ν−1Γ(ν) Kν(β∥t∥), ∀ t ∈ Rd, (4.3) where ν, β and σ2 are strictly positive, and Kν is the modified Bessel function of the second kind. Two covariance functions C(t) and CM (t) exhibit the same behavior at the origin if α = 2ν. Specifically, σ2 − C(t) ≍ σ2 − CM (t) ≍ ∥t∥α, (4.4) as ∥t∥ → 0 given α = 2ν. As a result, the realizations of X(t) have fractal dimension (d + 1 − α/2) with probability 1 (Gneiting and Schlather 2004; Adler 2010). However, Contrary to the fact that for d ≤ 3, σ2 and β in Eq. (4.3) cannot be estimated consistently under fixed-domain asymptotics proven by Zhang (2004), Anderes (2010) showed that both σ2 and θ in Eq. (4.1) can be consistently estimated given α ∈ (0, d/2) \ {1} under fixed- domain asymptotics. This shows that the powered exponential model and the Matérn model 58 have different statistical properties when d ≤ 3. As noted in Zhang (2004), the spectral densities corresponding to the powered exponential model do not have a closed form except in some special cases, which brings challenges when deriving results on the equivalence of Gaussian measures. The remainder of this chapter is organized as follows. In Section 4.2, we characterize the equivalence of Gaussian measure under the powered exponential model when α ∈ (d/2, 2). In Section 4.3, we establish the strong consistency and asymptotic normality of the MLE for σ2θ. In Section 4.4, we provide simulations of the MLE under finite sample cases. 4.2 The equivalence of Gaussian measures To study sufficient conditions for the equivalence of two Gaussian measures for the pow- ered exponential class, we use Theorem A.1 in M. L. Stein (2004): Let Pi, i = 1, 2, be two probability measures such that under Pi, the process X(s), s ∈ Rd, is stationary Gaussian with mean 0 and a second-order spectral density fi(v), v ∈ Rd. If, for some α > d, f1(v)|v|α is bounded away from 0 and ∞ as |v| → ∞, and for some finite c, (cid:90) |v|>c (cid:26) f1(v) − f2(v) f1(v) (cid:27)2 dv < ∞, (4.5) then for any bounded subset D ⊂ Rd, P1 ≡ P2 on the paths of X(s), s ∈ D. Since Eq. (4.1) and the characteristic function of a stable distribution have the same form, we will take advantage of the knowledge from stable distributions to verify conditions in Theorem A.1 in M. L. Stein (2004) for the powered exponential family, which is inspired by Formula (3.4) in Wang (2010). In Samorodnitsky and Taqqu (2017), a stable distribution is defined as follows: Definition 4.1. A random variable X is said to have a stable distribution Sα(σ, β, µ) if there are parameters 0 < α ≤ 2, σ ≥ 0, −1 ≤ β ≤ 1, and µ real such that its characteristic function has the follow form: E exp{iωX} =    exp(cid:8)−σα|ω|α(1 − iβ(sign ω) tan πα 2 ) + iµω(cid:9), if α ̸= 1, exp(cid:8)−σ|ω|(1 + iβ 2 π (sign ω) ln |ω|) + iµω(cid:9), if α = 1, (4.6) 59 where sign ω =    1, 0, if ω > 0, if ω = 0, −1, if ω < 0. Now consider X ∼ Sα(1, 0, 0), 0 < α < 1. By Eq. (4.6), it has ch.f. ψ(ω) = E exp{iωX} = exp{−|ω|α}. And its density function can be expanded into a convergent series as follows (see, e.g., Bergström 1952, Equation (4); Feller 1991, Lemma 1, p. 583) f (x) = 1 2π (cid:90) ∞ −∞ ψ(ω) exp{−iωx}dω = − 1 π ∞ (cid:88) k=1 (−1)k k! Γ(αk + 1) xαk+1 (cid:16) sin k (cid:17) , απ 2 x > 0. Obviously, for θ > 0, θ 1 α X has ch.f. Meanwhile ψθ(ω) = exp{−θ|ω|α}. fθ(x) = θ− 1 α x) α f (θ− 1 ∞ (cid:88) 1 π k=1 (−1)k k! = −θ− 1 α = θ Γ(α + 1) π sin (cid:16)απ 2 (cid:16) θ sin αk+1 α απ Γ(αk + 1) xαk+1 2 x−(α+1) − θ2 Γ(2α + 1) (cid:17) k 2π (cid:17) (4.7) sin(απ)x−(2α+1) + O(x−(3α+1)), as x → ∞. And for 1 < α < 2, we have the same asymptotic expansion Eq. (4.7) according to Eq. (5) in Bergström (1952). As for the multidimensional case, let Y be a centered d-dimensional isotropic stable random vector with characteristic function exp{−θ|u|α}. The amplitude of Y is defined by R = ∥Y∥ = (cid:113) Y2 1 + · · · + Y2 d. Then the density function of R has the following asymptotic series expansion: fR(r) = θ πΓ(d/2) θ2 2πΓ(d/2) − (cid:19) (cid:18) α + 2 Γ 2 (cid:19) (cid:18) α + d Γ 2 (cid:18) 2α + d 2 sin (cid:16)απ 2 (cid:19) (cid:17)(cid:18) 2 r (cid:18) 2 r Γ(α + 1)Γ sin(απ) (cid:19)α+1 (cid:19)2α+1 + O(r−3α−1), (4.8) 60 as r → ∞ for α ∈ (0, 1) ∪ (1, 2). Meanwhile for y ̸= 0, fY(y) = Γ(d/2) 2πd/2 |y|1−dfR(|y|); (4.9) (see Nolan 2005, Section 3.1 and 3.2 for details). And direct calculation shows that the asymptotic expression fY(y) is the same as Eq. (4.7) when d = 1. Lemma 4.1. Let Pi, i = 1, 2, be two probability measures such that under Pi, the process X(s), s ∈ Ω ⊂ Rd, is stationary Gaussian with mean 0 and a covariance function Ci(t) = i e−θi||t||α with parameters (σi, θi, α) where α ∈ (0, 2). If α ∈ (cid:0) d 2 , 2(cid:1), then for any bounded σ2 infinite set Ω ⊂ Rd, the Gaussian measures P1 and P2 are equivalent if and only if σ2 1θ1 = σ2 2θ2. Proof. Note that when α = 1 and d = 1, Theorem 2 in Zhang 2004 implies the theorem holds. Below, we only consider the case α ∈ (cid:0) d 2 , 2(cid:1) \ {1}. ” ⇐ ” Here we follow the proof of Theorem 2 in Zhang 2004. For powered exponential covariance function C(t) = σ2e−θ||t||α in Rd, the corresponding isotropic spectral density has the following asymptotic expansion: f (u) ≜ σ2fY((u, 0, · · · , 0)) = σ2 Γ(d/2) = σ2 Γ(d/2) (cid:32) 2πd/2 u1−dfR(u) θ πΓ(d/2) 2πd/2 u1−d θ2 2πΓ(d/2) 2α α + 2 πd/2+1 Γ( 2 22α−1 1θ2 πd/2+1 Γ(α + 1)Γ( 1 Γ(α + 1)Γ )Γ( − = σ2 1θ1 + σ2 (cid:18) α + 2 Γ 2 (cid:18) 2α + d 2 (cid:19) (cid:19) Γ (cid:18)α + d 2 (cid:18) 2 u sin(απ) (cid:19) sin (cid:16)απ 2 (cid:17)(cid:18) 2 u (cid:19)α+1 (cid:19)2α+1 (cid:33) + O(u−3α−1) (cid:17) u−(α+d) (cid:16) απ 2 ) sin α + d 2 2α + d 2 ) sin(απ)u−(2α+d) + O(cid:0)u−(3α+d)(cid:1), as u → ∞. Moreover, Eq. (4.5) can be expressed as (cid:90) ∞ c ud−1 (cid:26) f1(u) − f2(u) f1(u) (cid:27)2 du < ∞, 61 (4.10) (4.11) where f1 and f2 are isotropic spectral densities corresponding to P1 and P2. Assuming that 1θ1 = σ2 σ2 2θ2, then lim u→∞ uα+df1(u) = lim u→∞ uα+dσ2 1 Γ(d/2) 2πd/2 u1−d α + 2 2 )Γ( θ1 πΓ(d/2) α + d 2 ) sin Γ( α + 2 2 (cid:17) (cid:16) απ 2 . 2α πd/2+1 Γ( =σ2 1θ1 And )Γ( α + d 2 ) sin (cid:16) απ 2 (cid:17) ( 2 u )α+1 (cid:27)2 (cid:26) σ2 (cid:26) f1(u) − f2(u) f1(u) 1fR,1(u) − σ2 σ2 1fR,1(u) 2 − σ2 2θ2 2α−1 (σ2 σ2 1θ1 (cid:32) = ≲ 2fR,2(u) (cid:27)2 1θ2 1) Γ(α + 1)Γ( 2α+d Γ( α+2 2 )Γ( α+d 2 2 ) sin(cid:0) απ 2 ) sin(απ) (cid:1) u−α (4.12) (cid:33)2 , for all u big enough. So the integral in Eq. (4.11) is finite when α ∈ (cid:0) d 2 , 2(cid:1) \ {1}. Therefore, the two measures are equivalent. ” ⇒ ” If σ2 1θ1 ̸= σ2 tial covariance functions C(t; σ2 2θ2, let σ2 0 = σ2 2θ2/θ1. Then σ2 0, θ1, α) and C(t; σ2 2θ2 and the two powered exponen- 0θ1 = σ2 2, θ2, α) define two equivalent measures. 0, θ1, α) and We only need to show that for any countable infinite subset T ⊂ Ω, C(t; σ2 C(t; σ2 1, θ1, α) define two orthogonal Gaussian measures on σ{X(t), t ∈ T }. For Gaussian random variables {X(t), t ∈ T }, we re-index them as {X(n), n ∈ N}. Obviously, any finite many of them has a positive definite covariance matrix with respect to both Gaussian measures. Applying Gram-Schmidt orthogonalization to {X(n), n ∈ N} with the inner product in- duced by C(t; σ2 0, θ1, α), we get {ξi}∞ i=1 which is a sequence of independent standard Gaussian random variables with respect to the Gaussian measure induced by C(t; σ2 0, θ1, α). And {ξi}∞ i=1 is a sequence of independent Gaussian random variables with mean zero and variance σ2 1/σ2 0 with respect to the Gaussian measure induced by C(t; σ2 1, θ1, α). 62 Let ϵ ≜ |σ2 1/σ2 0−1| 2 , then by law of large numbers, we have lim n→∞ Pσ2 0,θ1,α (cid:32)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) k=1 (cid:12) (cid:12) ξ2 (cid:12) k/n − 1 (cid:12) (cid:12) (cid:33) > ϵ = 0. (4.13) Meanwhile, 1 = lim n→∞ Pσ2 1,θ1,α (cid:32)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) k=1 ξ2 k/n − (cid:12) σ2 (cid:12) 1 (cid:12) (cid:12) σ2 (cid:12) 0 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) < (cid:12) (cid:12)σ2 1/σ2 0 − 1(cid:12) (cid:12) − ϵ (cid:33) n (cid:88) k=1 ξ2 k/n − (cid:33) > ϵ (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) σ2 1 σ2 0 (cid:32) = lim n→∞ Pσ2 1,θ1,α (cid:12) (cid:12)σ2 1/σ2 0 − 1(cid:12) (cid:12) − ≤ lim n→∞ Pσ2 1,θ1,α (cid:32)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) k=1 (cid:12) (cid:12) ξ2 (cid:12) k/n − 1 (cid:12) (cid:12) (cid:33) . > ϵ Then lim n→∞ Pσ2 1,θ1,α (cid:32)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) k=1 (cid:12) (cid:12) ξ2 (cid:12) k/n − 1 (cid:12) (cid:12) (cid:33) > ϵ = 1. Finally, combining Eq. (4.13) and Eq. (4.14), we have Pσ2 thus on σ{X(t), t ∈ Ω}. 0,θ1,α ⊥ Pσ2 1,θ1,α (4.14) on σ{X(t), t ∈ T }, Remark 4.1. For d-dimensional isotropic Ornstein-Uhlenbeck process with covariance func- tion Vd(t, s; σ2, θ) = σ2e−θ[(cid:80)d i=1(ti−si)2] 1 2 , setting the smoothness parameter ν as 1 2 in the Matérn model, we have that for d ≤ 3, Vd(·; σ2 1, θ1) and Vd(·; σ2 2, θ2) induces two Gaussian measures that are equivalent if and only if σ2 1θ1 = σ2 2θ2 according to Theorem 2 in Zhang (2004). With Lemma 4.1, we see that σ and θ cannot be consistently estimated when α ∈ (d/2, 2) under fixed-domain asymptotics (see Zhang 2004, for a detailed explanation). However, as we will show in the next section, the MLE of σ2θ is consistent and asymptotic normal. 4.3 Asymptotic properties of the MLE Let the observation of Gaussian field X(t) be Xn = {X(t1), X(t2), · · · , X(tn)}′, (4.15) 63 where t1, t2, · · · , tn ∈ [0, T ]d are distinct locations, and [0, T ]d ⊂ Rd is the fixed cube with side T . Theorem 4.1. Let {X(s), s ∈ Rd}, d ∈ {1, 2, 3}, be stationary isotropic Gaussian field with mean 0 and covariance function C(t) = σ2 0e−θ0||t||α where σ2 0 and θ0 are unknown and α is known. Let Dn, n = 1, 2, · · · , be an increasing (nested) sequence of finite subsets of [0, T ]d ⊂ Rd, and Ln(σ2, θ) be the likelihood function when the process is observed at locations in Dn. For any fixed θ1 > 0, let ˆσ2 n maximize Ln(σ2, θ1). Then when α ∈ ( d 2 , 2), ˆσ2 nθ1 → σ2 0θ0, with P0 probability 1, where P0 is the Gaussian measure defined by the powered exponential covariance function corresponding to parameter values σ2 0, θ0 and α. Namely, X′ nΓ−1 (n,θ1)Xn n θ1 P0−a.s. −−−−→ n→∞ σ2 0θ0, (4.16) where Γ(n,θ0) is the true correlation matrix. Proof. The proof follows the same arguments of the proof of Theorem 3 in Zhang (2004). To derive the asymptotic normality, the key is to control the error between quadratic forms with the misspecified correlation matrix in Eq. (4.16) and the true one respectively. Before stating the main result on the asymptotic normality, we need to introduce several preliminaries that will be used in the main proof. Let σ2 1 = σ2 covariance function C(t; σ2 0θ0/θ1, and P1 is the Gaussian measure corresponding to mean 0 and the 1, θ1, α) in the following. First, we observe that by the eigende- composition (cid:0)σ2 0Γ(n,θ0) (cid:1)−1/2(cid:0)σ2 1Γ(n,θ1) (cid:1)(cid:0)σ2 0Γ(n,θ0) (cid:1)−1/2 = UnΛnU ′ n, (4.17) where Un is an orthogonal matrix consisting of the eigenvectors and Λn = diag(λ1,n, · · · , λn,n) is a diagonal matrix with the eigenvalues {λk,n}n k=1 . Define Yn = U ′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1/2Xn, (4.18) 64 then (cid:16) 1 √ n X′ n (cid:0)σ2 1Γ(n,θ1) (cid:1)−1Xn − X′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn (cid:17) = 1 √ n n (cid:88) k=1 (cid:0)λ−1 k,n − 1(cid:1)Y 2 k,n, where (Y1,n, · · · , Yn,n)′ = Yn ∼ N (0, In) under P0, and Yn ∼ N (0, Λn) under P1. Second, we plan to represent Yn in the frequency domain. The spectral representation theorem (see, e.g., Adler and Taylor 2007, Theorem 2.4.2) says X(t) has the following spectral representation X(t) fdd= (cid:90) Rd eiω′tM (dω), (4.19) where M is a Gaussian random measure on Rd whose control measure has the Radon- Nikodym derivative fσ,θ with respect to the Lebesgue measure. The fdd= means that both sides in Eq. (4.19) have the same finite-dimensional distributions. Subsequently, we introduce two isomorphic Hilbert spaces. Define and H 0 ≜ L0 ≜ (cid:40) n (cid:88) i=1 (cid:40) n (cid:88) i=1 riX(ti) : r1, · · · , rn ∈ R, t1, · · · , tn ∈ [0, T ]d (cid:41) , rieiω′ti : r1, · · · , rn ∈ R, t1, · · · , tn ∈ [0, T ]d (cid:41) . (4.20) (4.21) Eq. (4.19) implies that cov (cid:40) n (cid:88) j=1 ajX(tj), m (cid:88) k=1 (cid:41) bkX(tk) = (cid:32) n (cid:88) j=1 aj exp{iω′tj}, m (cid:88) k=1 (cid:33) bk exp{iω′tk} , (4.22) fσ,θ where the inner product (·, ·)fσ,θ is defined as (ϕ, ψ)fσ,θ Rd ϕ(ω)ψ(ω)dω. Therefore, two Hilbert spaces H(fσ,θ), as the closure of H 0 with respect to the inner product cov{·, ·}, and L(fσ,θ), as the closure of L0 with respect to the inner product (·, ·)fσ,θ e.g., M. L. Stein 1999a, chap. 2.6; Ibragimov and Rozanov 1978, chap. 3.1.3). Define , are isomorphic; (see, ≜ (cid:82) (ψ1, · · · , ψn)′ = U ′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1/2(exp{iω′t1}, · · · , exp{iω′tn})′. (4.23) Combining Eqs. (4.18) and (4.22), we see that for j, k = 1, · · · , n, (ψj, ψk)fσ0,θ0 = δj,k, (ψj, ψk)fσ1,θ1 = λj,nδj,k, 65 (4.24) where δj,k = 1 if j = k and is 0 otherwise. Then, to bound {λk,n}n k=1, we first introduce the entropy distance between two equivalent Gaussian measures (see M. L. Stein (1999a) and Ibragimov and Rozanov (1978) for a fuller exposition about the entropy distance). Lemma 4.1 implies P0 ≡ P1, and the entropy distance r between P1 and P0 is defined as (cid:20) r ≜ − E0 ln (cid:19) (cid:18) P1 P0 + E1 ln (cid:19)(cid:21) , (cid:18) P0 P1 (4.25) where P0 P1 and P1 P0 are the Radon-Nikodym derivative on the σ-algebra generated by the Gaussian process X. Moreover, the conditional expectation of P1 P0 given Xn is (cid:19) (cid:12) (cid:12) (cid:12) (cid:12) Xn pn ≜ E0 (cid:18) P1 P0 f1(Xn) f0(Xn) (cid:32)(cid:12) (cid:0)σ2 (cid:12) (cid:12) (cid:0)σ2 (cid:12) = = (cid:33)1/2 (cid:1)(cid:12) (cid:12) (cid:1)(cid:12) (cid:12) e− 1 2 X′ n(σ2 1Γ(n,θ1))−1 Xn+ 1 2 X′ n(σ2 0Γ(n,θ0))−1 (4.26) Xn, 0Γ(n,θ0) 1Γ(n,θ1) where f1(x) and f0(x) are the probability density functions of Xn with respect to P1 and P0 respectively. Define rn as rn ≜ − E0 ln (pn) + E1 ln (pn) (cid:16)(cid:0)σ2 E0 X′ n = 1 2 1Γ(n,θ1) (cid:16)(cid:0)σ2 0Γ(n,θ0) (cid:1)−1 − (cid:0)σ2 0Γ(n,θ0) (cid:1)−1 − (cid:0)σ2 (cid:1)−1(cid:17) Xn (cid:1)−1(cid:17) Xn + 1 2 1 2 = = 1 2 (cid:104) E1 X′ n (cid:16)(cid:0)σ2 (cid:18) 1 λk,n tr n (cid:88) k=1 1Γ(n,θ1) (cid:1)−1(cid:0)σ2 0Γ(n,θ0) (cid:19) , + λk,n − 2 1Γ(n,θ1) (cid:1)(cid:17) + tr (cid:16)(cid:0)σ2 0Γ(n,θ0) (cid:1)−1(cid:0)σ2 1Γ(n,θ1) (cid:1)(cid:17) (cid:105) − 2n (4.27) where {λk,n}n k=1 are defined in Eq. (4.17). Lemma 4.2. There exist two positive constants 0 < Ci ≤ Cs < ∞ such that Ci ≤ λk,n ≤ Cs for all n and 1 ≤ k ≤ n. 66 Proof. Theorem 4 in M. L. Stein (1999a, p.117) implies 0 ≤ r < +∞ in Eq. (4.25). And by Jensen’s inequality, rn ≤ r for all n, which implies that (cid:18) 1 λk,n 1 2 (cid:19) + λk,n − 2 ≤ r, for all n and 1 ≤ k ≤ n. Let Ci = (2 + 2r) − (cid:112)(2 + 2r)2 − 4 2 , Cs = (2 + 2r) + (cid:112)(2 + 2r)2 − 4 2 , (4.28) then we have This proves Lemma 4.2. Ci ≤ λk,n ≤ Cs. Finally, we describe properties of the spectral density function fσ,θ(ω) in the following lemma. Lemma 4.3. The spectral density fσ,θ(ω) defined in Eq. (4.2) is bounded, uniformly contin- uous and strictly positive on Rd. Furthermore, fσ,θ(ω) ≍ ∥ω∥−(α+d) as ∥ω∥ → ∞. Proof. Since (cid:82) Rd C(t)dt < ∞, fσ,θ(ω) is bounded and uniformly continuous by the properties of Fourier transform. Furthermore, based on Formula (7.5) in Khoshnevisan, Xiao, and Zhong (2003), for ω ̸= 0, f1,1(ω) = ∥ω∥−d (cid:90) ∞ ν 0 (cid:18) s ∥ω∥2 (cid:19) g(α/2)(s)ds, (4.29) where the function ν is defined as ν(s) = (4πs)−d/2e− 1 4s , g(α/2)(s) is the density function of the random variable τ (1), and where τ = τ (t), t ≥ 0 is a stable subordinator of index α 2 . Finally, f1,1(0) = 1 (2π)d (cid:90) Rd C(t)dt > 0. So we have f1,1(ω) > 0. The fact that fσ,θ(ω) ≍ ∥ω∥−(α+d) as ∥ω∥ → ∞ has been verified in the proof of Lemma 4.1. This completes the proof. 67 Theorem 4.2. With the same notation and assumptions as in Theorem 4.1, we have √ with respect to P0. n(cid:0)ˆσ2 nθ1 − σ2 0θ0 (cid:1) d−→ N (cid:0)0, 2(σ2 0θ0)2(cid:1), (4.30) Proof. Note that Ying (1991) has proven the theorem when α = 1 and d = 1. Therefore, we only consider the case α ∈ (cid:0) d 2 , 2(cid:1) \ {1}, and we follow the arguments from Wang (2010) in the following. Splitting (ˆσ2 nθ1 − σ2 0θ0) into two parts, we have √ n(cid:0)ˆσ2 nθ1 − σ2 0θ0 (cid:1) = (cid:16) σ2 0θ0√ n X′ n (cid:0)σ2 1Γ(n,θ1) (cid:1)−1Xn − X′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn (cid:17) + (cid:16) σ2 0θ0√ n X′ n (cid:0)σ2 0Γ(n,θ0) (cid:17) (cid:1)−1Xn − n . (4.31) (cid:1)−1/2Xn consists of i.i.d. standard normal variables with respect to P0, (cid:1)−1Xn is the sum of i.i.d. variables having χ2 1 distribution. Hence, the Lindeberg- 0Γ(n,θ0) Since (cid:0)σ2 (cid:0)σ2 X′ n 0Γ(n,θ0) Feller theorem implies (cid:16) σ2 0θ0√ n X′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn − n (cid:17) d−→ N (cid:0)0, 2(σ2 0θ0)2(cid:1). Thus to prove Theorem 4.2, it suffices to show that (cid:16) σ2 0θ0√ n X′ n (cid:0)σ2 1Γ(n,θ1) (cid:1)−1Xn − X′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn (cid:17) P0−→ 0. (4.32) (4.33) Let m = (cid:4) α+d 2 (cid:5) + 1 and κ = α+d 4m ∈ (0, 1 2), where ⌊·⌋ denotes the greatest integer function. Define the integrable function c0 and its Fourier transform ξ0 as c0(x) = ∥x∥κ−dI(∥x∥ ≤ 1), x ∈ Rd, ξ0(ω) = (cid:90) Rd e−ix′ωc0(x)dx, x ∈ Rd. (4.34) (4.35) It follows from Lemma 4.5 that ξ0 : Rd → R is a continuous, isotropic and strictly positive function and ξ0(ω) ≍ ∥ω∥−κ as ∥ω∥ → ∞. Let c1 = c0 ∗ · · · ∗ c0 denote the 2m-fold convolution of the function c0 with itself. Then supp(c1) ⊂ {x : ∥x∥ ≤ 2m} and ξ1(ω) = ξ0(ω)2m. (4.36) 68 It follows from Lemma 4.3 that there exist positive constants 0 < cξ1 ≤ Cξ1 < ∞ such that for all ω ∈ Rd, Define cξ1 ≤ fσ0,θ0(ω) ξ1(ω)2 ≤ Cξ1. η(ω) = fσ1,θ1(ω) − fσ0,θ0(ω) ξ1(ω)2 . (4.37) (4.38) It follows from the proof of Lemma 4.1 that there exists a constant Cη such that for all ω ∈ Rd, |η(ω)| ≤ Cη 1 + ∥ω∥α , (4.39) and η : Rd → R is square integrable when α ∈ (cid:0) d 2 , 2(cid:1) \ {1}. From the theory of the Fourier transform of L2(Rd) functions, there exists a square integrable function g : Rd → R such that where (cid:90) Rd (η(ω) − ˆgk(ω))2dω k→∞−−−→ 0, ˆgk(ω) = (cid:90) Rd e−iω′xg(x)I{∥x∥max ≤ k}dx, k ∈ N. (4.40) (4.41) In order to illustrate the following Lemma, we need to introduce the approximate identity en(x) defined by Formula (35) of Wang and Loh (2011). Let {ϵn}∞ n=1 be a positive sequence such that ϵn → 0 as n → ∞. Define and where Ce = (cid:82) en(x) = 1 Ceϵd n (cid:18) x ϵn ˜c1 (cid:19) , x ∈ Rd, ˜ξ1(ω) = (cid:90) Rd e−ix′ω˜c1(x)dx, x ∈ Rd, (4.42) (4.43) Rd ˜c1(x)dx and ˜c1 = ˜c0 ∗ · · · ∗ ˜c0, the 2ma-fold convolution of ˜c0, with ˜c0(x) = (cid:5) + 1. Here a is an arbitrary positive constant. For the −dI{∥x∥ ≤ 1} and ma = (cid:4) a+d a+d 4ma ∥x∥ 2 Fourier transform of en, we have ˆen(ω) = ˜ξ1(ϵnω) Ce . 69 (4.44) Lemma 4.5 implies that there exists a constant Cˆe (independent of ω and n) such that for all ω ∈ Rd, |ˆen(ω)| ≤ Cˆe (1 + ϵn∥ω∥)(a+d)/2 . (4.45) Lemma 4.4. Let en and g be as in Eqs. (4.41) and (4.42) respectively and β0 is any constant satisfying 0 < β0 < min{2, 2α − d}. Then there exists a constant Cβ0 such that (cid:90) Rd |en ∗ g(x) − g(x)|2dx ≤ Cβ0ϵβ0 n . (4.46) Proof. The Plancherel theorem implies that for y ∈ Rd, (cid:90) Rd |g(x − y) − g(x)|2dx = = ≤ (cid:90) Rd (cid:90) 1 (2π)d 1 (2π)d 22−β0∥y∥β0 (2π)d Rd (cid:12) (cid:12)e−iω′yη(ω) − η(ω) (cid:12) (cid:12) (cid:12) (cid:12)(e−iω′y − 1)η(ω) (cid:12) (cid:12) (cid:12) (cid:90) 2 2 (cid:12) (cid:12) (cid:12) dω dω ∥ω∥β0|η(ω)|2dω. Rd (4.47) Hence (cid:20)(cid:90) Rd |en ∗ g(x) − g(x)|2dx (cid:21)1/2 [g(x − y) − g(x)]en(y)dy (cid:35)1/2 2 (cid:12) (cid:12) (cid:12) (cid:12) dx |g(x − y) − g(x)|2dx (cid:21)1/2 en(y)dy (4.48) (cid:34)(cid:90) = Rd (cid:12) (cid:12) (cid:12) (cid:12) (cid:90) ≤ (cid:90) ∥y∥≤2maϵn (cid:20)(cid:90) ∥y∥≤2maϵn Rd ≤ 2(2−β0)/2(2maϵn)β0/2 (2π)d/2 ≤ Cη2(2−β0)/2(2maϵn)β0/2 (2π)d/2 This proves Lemma 4.4. (cid:20)(cid:90) ∥ω∥β0|η(ω)|2dω (cid:21)1/2 Rd (cid:34)(cid:90) Rd ∥ω∥β0 (1 + ∥ω∥α)2 dω (cid:35)1/2 . 70 = = (cid:90) Rd (cid:90) Rd (cid:90) R2d (cid:90) = (2π)d = (2π)d R2d (cid:90) + (2π)d R2d (cid:90) = (2π)d R2d (cid:90) + (2π)d Let x, y ∈ [0, T ]d, and observing that supp(c1) ⊂ [−2m, 2m]d, we obtain b(x, y) ≜ Efσ1,θ1 [X(x)X(y)] − Efσ0,θ0 [X(x)X(y)] ei(x−y)′ω[fσ1,θ1(ω) − fσ0,θ0(ω)]dω ei(x−y)′ωη(ω)ξ1(ω)2dω g(s − t)c1(x − s)c1(y − t)dsdt en ∗ g(s − t)c1(x − s)c1(y − t)dsdt (4.49) [g(s − t) − en ∗ g(s − t)]c1(x − s)c1(y − t)dsdt en ∗ g(s − t)c1(x − s)c1(y − t)dsdt h∗ n(s, t)c1(x − s)c1(y − t)dsdt, where R2d h∗ n(s, t) = [g(s − t) − en ∗ g(s − t)]I{∥s + t∥max ≤ 4m + 2T }, ∀ s, t ∈ Rd. Let η∗ n : Rd → C denote the L2(Rd) Fourier transform of g − en ∗ g. This implies that (cid:90) Rd (cid:0)η∗ n(ω) − ˆg∗ n,k(ω)(cid:1)2dω k→∞−−−→ 0, (4.50) where ˆg∗ n,k(ω) = (cid:90) Rd e−iω′x[g(x) − en ∗ g(x)]I{|x|max ≤ k}dx, k ∈ N. (4.51) Thus as in (24) in Wang and Loh (2011), we have (cid:90) (2π)d R2d (cid:90) =(2π)−d h∗ n(s, t)c1(x − s)c1(y − t)dsdt ei(ω′x−ν′y)η∗ n (cid:18) ω + ν 2 (cid:19) (cid:18) ω − ν κ 2 (cid:19) ξ1(ω)ξ1(ν)dωdν, R2d (4.52) where κ(x) = 2−d (cid:82) Rd e−ix′tI{∥t∥max ≤ 4m + 2T }dt, x ∈ Rd. We observe that since I{∥t∥max ≤ 4m + 2T } ∈ L1 ∩ L2, κ is continuous and κ ∈ L2(Rd). Now we define h∗∗ n (s, t) = (cid:90) ∥u∥max≤2m+2ma+T en(s − u)g(u − t)du, ∀ s, t ∈ Rd. 71 Then the function h∗∗ n : Rd → C is square-integrable and (cid:90) (2π)d R2d (cid:90) =(2π)d R2d (cid:90) =(2π)−d en ∗ g(s − t)c1(x − s)c1(y − t)dsdt h∗∗ n (s, t)c1(x − s)c1(y − t)dsdt ei(ω′x−ν′y)ξ1(ω)ξ1(ν) e−i(ω′u−ν′u)ˆen(ω)η(ν)du dωdν. (cid:33) (4.53) R2d (cid:32)(cid:90) × ∥u∥max≤2m+2ma+T It follows from Eqs. (4.52) and (4.53), that for x, y ∈ [0, T ]d, b(x, y) = (2π)−d (cid:90) R2d (cid:90) + (2π)−d ei(ω′x−ν′y)η∗ n (cid:18) ω + ν 2 (cid:19) (cid:18) ω − ν κ 2 (cid:19) ξ1(ω)ξ1(ν)dωdν R2d (cid:32)(cid:90) × ei(ω′x−ν′y)ξ1(ω)ξ1(ν) (4.54) e−i(ω′u−ν′u)ˆen(ω)η(ν)du dωdν. (cid:33) ∥u∥max≤2m+2ma+T Let {ψ1, · · · , ψn} be as in Eq. (4.23), and then for k = 1, · · · , n, (ψk, ψk)fσ1,θ1 − (ψk, ψk)fσ0,θ0 = λk,n − 1 = ν∗ k,n + ν∗∗ k,n, (4.55) where and ν∗ k,n = (2π)−d (cid:90) R2d ψk(ω)ψk(ν)η∗ n (cid:18) ω + ν 2 (cid:19) (cid:18) ω − ν κ 2 (cid:19) ξ1(ω)ξ1(ν)dωdν, (cid:90) R2d k,n = (2π)−d ν∗∗ (cid:32)(cid:90) × ψk(ω)ψk(ν)ξ1(ω)ξ1(ν) e−i(ω′u−ν′u)ˆen(ω)η(ν)du dωdν. (cid:33) ∥u∥max≤2m+2ma+T 72 Using Bessel’s inequality, we have n (cid:88) k=1 (cid:12) (cid:12)ν∗∗ k,n (cid:12) (cid:12) ≤ (2π)−d × (cid:12) (cid:12) (cid:12) (cid:12) (cid:90) Rd ≤ 1 2(2π)d (cid:90) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:40)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) (cid:90) k=1 ∥u∥max≤2m+2ma+T e−iω′uψk(ω)ξ1(ω)ˆen(ω)dω (cid:12) (cid:12) (cid:12) (cid:12) Rd eiν′uψk(ν)ξ1(ν)η(ν)dν du (cid:90) n (cid:88) ∥u∥max≤2m+2ma+T k=1 (cid:90) Rd e−iω′uψk(ω) ξ1(ω) f 1/2 σ0,θ0 (ω) ˆen(ω) × f 1/2 σ0,θ0 (ω)dω 2 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) + (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:90) Rd eiν′uψk(ν) ξ1(ν) (ν) η(ν)f 1/2 σ0,θ0 (ν)dν 2(cid:41) du (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (4.56) (cid:90) (cid:90) ∥u∥max≤2m+2ma+T (cid:90) Rd (cid:90) ≤ 1 2(2π)d + 1 2(2π)d (2m + 2ma + T )d 2πd ∥u∥max≤2m+2ma+T (cid:26) ξ2 1(s) fσ0,θ0(s) sup s∈Rd Rd f 1/2 σ0,θ0 ξ2 1(ω) fσ0,θ0(ω) ξ2 1(ν) fσ0,θ0(ν) (cid:27) (cid:90) Rd |ˆen(ω)|2dωdu |η(ν)|2dνdu |ˆen(x)|2 + |η(x)|2dx, ≤ and n (cid:88) k=1 (cid:12) (cid:12)ν∗ k,n (cid:12) (cid:12) 2 ≤ (2π)−d (cid:26) sup s∈Rd ξ2 1(s) fσ0,θ0(s) (cid:27)2 (cid:90) R2d (cid:12) (cid:12) (cid:12) (cid:12) η∗ n (cid:18)ω + ν 2 (cid:19) (cid:18) ω − ν κ 2 2 (cid:19)(cid:12) (cid:12) (cid:12) (cid:12) dωdν = π−d (cid:26) ξ2 1(s) fσ0,θ0(s) sup s∈Rd (cid:27)2 (cid:90) Rd |η∗ n(ω)|2dω (cid:90) Rd |κ(ν)|2dν. From Eqs. (4.37) and (4.45), there exists constant C2 independent of n such that n (cid:88) k=1 (cid:12) (cid:12)ν∗∗ k,n (cid:12) (cid:12) ≤ C2 ϵa+d n . (4.57) (4.58) Meanwhile, Combining Lemma 4.4 and Eq. (4.37), we observe that there exists constant C1 independent of n such that n (cid:88) (cid:12) (cid:12)ν∗ k,n (cid:12) (cid:12) 2 ≤ C1ϵβ0 n . Moreover, Jensen’s inequality implies that k=1 n (cid:88) k=1 (cid:12) (cid:12)ν∗ k,n (cid:12) (cid:12) ≤ (cid:32) n n (cid:88) k=1 (cid:33)1/2 (cid:12) (cid:12)ν∗ k,n (cid:12) (cid:12) 2 (cid:113) C1nϵβ0 n . ≤ So we conclude that n (cid:88) k=1 |λk,n − 1| ≤ n (cid:88) (cid:0)(cid:12) (cid:12)ν∗ k,n (cid:12) + (cid:12) (cid:12) (cid:12)ν∗∗ k,n (cid:12) (cid:12) (cid:1) k=1 (cid:113) C1nϵβ0 n + ≤ C2 ϵa+d n . 73 (4.59) (4.60) (4.61) Finally for any constant ϑ > 0, using Markov’s inequality, Lemma 4.2 and Eq. (4.61) we obtain (cid:18) 1 √ n P0 (cid:32) 1 √ n =P0 (cid:12) (cid:12)X′ (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) n (cid:88) k=1 (cid:0)σ2 1Γ(n,θ1) n n (cid:1)−1Xn − X′ (cid:12) (cid:33) (cid:12) (cid:12) (cid:12) (cid:12) > ϑ k,n (cid:0)λ−1 k,n − 1(cid:1)Y 2 (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn (cid:19) (cid:12) (cid:12) (cid:12) > ϑ |λk,n − 1| (4.62) (cid:12) (cid:12)λ−1 k,n − 1(cid:12) (cid:12) max 1≤j≤n (cid:18)(cid:113) λ−1 j,n (cid:27) n (cid:88) k=1 C1nϵβ0 n + ≤ ≤ ≤ = n (cid:88) k=1 (cid:26) 1 √ ϑ n 1 √ ϑ n 1 √ Ciϑ C 1/2 1 n ϵβ0/2 n Ciϑ + C2 Ciϑn1/2ϵa+d . (cid:19) C2 ϵa+d n n Choose ϵn such that ϵn → 0 and n1/2ϵa+d n → ∞ as n → ∞. It follows from Eq. (4.62) that (cid:16) σ2 0θ0√ n X′ n (cid:0)σ2 1Γ(n,θ1) (cid:1)−1Xn − X′ n (cid:0)σ2 0Γ(n,θ0) (cid:1)−1Xn (cid:17) P0−→ 0. (4.63) This proves Theorem 4.2. Lemma 4.5. Let ξ0 be as in Eq. (4.35). Then ξ0 : Rd → R is a continuous, isotropic and strictly positive function and ξ0(ω) ≍ ∥ω∥−κ as ∥ω∥ → ∞. Proof. We split the proof into two parts. When d = 1, we follow the proof of Lemma 5 in Du, Zhang, and Mandrekar (2009). Notice that ξ0 is continuous and real symmetric since c0 ∈ L1(R) and c0 is real symmetric. Meanwhile, ξ0(ω) = (cid:90) 1 −1 e−iωx|x|κ−1dx = 2 (cid:90) 1 0 cos(ωx)xκ−1dx. (4.64) When ω > 0, let u = ωx. We have ξ0(ω) = 2ω−κ (cid:90) ω cos(u)uκ−1du. (4.65) Notice that (cid:82) ∞ 0 cos(u)uκ−1du = Γ(κ) cos(cid:0) πκ 2 0 (cid:1). So ξ0(ω) ≍ |ω|−κ, as |ω| → ∞. 74 2kπ+5π/2 (cid:90) 5π/2 3π/2 (cid:90) π/2 cos(u)u−δdu As for the strict positiveness, first notice that ξ0(0) = (cid:82) 1 suffices to show −1 |x|κ−1dx > 0. For any ω > 0, it y(ω) = (cid:90) ω 0 cos(u)u−δdu > 0, (4.66) where δ = 1 − κ ∈ [1/2, 1). Note that y′(ω) = cos(ω)ω−δ and y′′(ω) = − sin(ω)ω−δ − δ cos(ω)ω−δ−1. Therefore, the minimum points over (0, ∞) are {2kπ + 3π/2, k = 0, 1, · · · }. We first claim y(3π/2) is the global minimum. Actually (cid:90) 2(k+1)π+3π/2 cos(u)u−δdu cos(u)u−δdu + (cid:90) 2(k+1)π+3π/2 cos(u)u−δdu (4.67) cos(s)(2kπ + s)−δds − cos(s)(2kπ + π + s)−δdu > 0. 2kπ+3π/2 (cid:90) 2kπ+5π/2 2kπ+3π/2 (cid:90) 5π/2 3π/2 = = And (cid:90) π/4 y(3π/2) = cos(u)u−δdu + 0 (cid:90) π + π/2 cos(u)u−δdu + π/4 (cid:90) 3π/2 π cos(u)u−δdu (cid:32) 1 − ≥ cos(π/4) (π/4)1−δ 1 − δ (cid:32) + √ 1 (π/2)δ (cid:33) √ √ 2 2 ≥ π + 1 − 2 π 2 2 (cid:114) 2 π − − (4.68) √ 2 2 1 √ π (cid:33) − 1 (π/2)δ − 1 πδ > 0. This completes the strict positiveness. When d ≥ 2, we follow the proof of Lemma 2 in the Supplement of Bevilacqua et al. (2019). Obviously, ξ0(0) > 0. Let Ud be the uniform probability measure on Sd−1 = {u ∈ 75 Rd : ∥u∥ = 1}. By isotropy, we have for all ω ∈ Rd \ {0}, (cid:90) (cid:90) ξ0(ω) = e−i∥ω∥x′u∥x∥κ−ddxUd(du) Sd−1 ∥x∥≤1 (cid:90) = ∥x∥≤1 Γ(d/2) = (2π)d/2∥ω∥ 2−d 2 (cid:18) 2 ∥ω∥∥x∥ (cid:90) 1 rκ− d (cid:19)(d−2)/2 J(d−2)/2(∥ω∥∥x∥)∥x∥κ−ddx 2 J(d−2)/2(∥ω∥r)dr (4.69) = (2π)d/2∥ω∥−κ 0 (cid:90) ∥ω∥ 0 rκ− d 2 J(d−2)/2(r)dr = 2κ−1πd/2Γ(d/2)−1 1F2(κ/2; κ/2 + 1, d/2; −(∥ω∥/2)2), where J(d−2)/2 is the Bessel function of the first kind, and the generalized hypergeometric function 1F2 is defined as 1F2(a; b, c; z) = ∞ (cid:88) k=0 (a)kzk (b)k(c)kk! , (4.70) with (q)k = Γ(q + k)/Γ(q) for k ∈ N ∪ {0}. See M. L. Stein (1999a, p.43) for the first and second equation in Eq. (4.69); the third equation is derived by the spherical coordinates transform; and see Prudnikov, Brychkov, and Marichev (1986, 1.8.1.1, p.37) for the last equation. From Theorem 2 in Fields and Ismail (1975), we get that ξ0 > 0 for d ≥ 2. Then ξ0 is a continuous, isotropic and strictly positive on Rd. Next, we prove the integral (cid:82) ∞ 0 rκ− d 2 J(d−2)/2(r)dr exists. Based on the series expansion (see, e.g., Abramowitz and Stegun 1968, p. 9.1.10) we have Jν(x) = (cid:16)x 2 (cid:17)ν ∞ (cid:88) k=0 (cid:0)− 1 4x2(cid:1)k k!Γ(ν + k + 1) , J(d−2)/2(r) = O(r(d−2)/2), as r → 0. So the integrand is integrable around the origin. Meanwhile, Based on the asymptotic expansion (see, e.g., Abramowitz and Stegun 1968, p. 9.2.1) Jν(x) = (cid:112)2/(πx) (cid:26) (cid:18) cos x − 1 2 νπ − (cid:19) 1 4 π (cid:27) + O(|x|−1) , as |x| → ∞, we have J(d−2)/2(r) = O(r−1/2), 76 as r → ∞. So the integrand is integrable around the infinity. Overall, (cid:82) ∞ 0 rκ− d 2 J(d−2)/2(r)dr exists. Therefore, ξ0(ω) ≍ ∥ω∥−κ as ∥ω∥ → ∞. This completes the proof. 4.4 Simulations In this section, we investigate the finite sample performance of the estimator ˆσ2 nθ1 of σ2 0θ0. We simulate, using the R package RandomFields developed by Schlather et al. (2017), and then we estimate with 1000 realizations from a zero mean Gaussian field with the covariance function Eq. (4.1) where • σ0 = 1 and θ1 = 1 are fixed. • θ0 = 2, 5 or 10. • For d = 2 and α = 1.2, we generated X(t) on a regular grid of [0, 1]2 with sample size n = 50m, that is, (cid:26) X (cid:19) (cid:18) i1 40 , i2 50 : 1 ≤ i1 ≤ m, 1 ≤ i2 ≤ 50 (cid:27) , where m = 20, 30 or 40. • For d = 3 and α = 1.8, we generated X(t) on a regular grid of [0, 1]3 with sample size n = 100m, that is, (cid:26) X (cid:18) i1 20 , i2 10 , i3 10 (cid:19) : 1 ≤ i1 ≤ m, 1 ≤ i2, i3 ≤ 10 (cid:27) , where m = 10, 15 or 20. Tables 4.1 and 4.2 summarize the percentiles of order 0.05, 0.25, 0.5, 0.75, 0.95, bias and sample standard deviation (SD) of ˆσ2 nθ1. Fig. 4.1 shows histograms of (cid:112)n/2(ˆσ2 nθ1/(σ2 0θ0) − 1) when d = 2 and n = 2000. Remark 4.2. As we see from Tables 4.1 and 4.2, the estimator ˆσ2 nθ1 has a larger bias when the fixed scale parameter θ1 is further away from the true parameter θ0. This phenomenon not only occurs under the powered exponential model but also other models. See, for example, 77 Table S-1 in the supplement of Kaufman and Shaby (2013) for the Matérn model and Table 2 in Bevilacqua et al. (2019) for the generalized Wendland model. θ0 2 5 10 θ0 2 5 10 n 1000 1500 2000 1000 1500 2000 1000 1500 2000 5% 1.853 1.887 1.903 4.677 4.761 4.797 9.533 9.635 9.718 25% 1.939 1.953 1.961 4.891 4.928 4.948 9.924 9.978 10.009 50% 2.006 2.007 2.007 5.056 5.056 5.056 10.243 10.245 10.236 75% 2.067 2.056 2.049 5.216 5.176 5.159 10.541 10.481 10.457 95% 2.157 2.126 2.112 5.458 5.368 5.320 10.983 10.850 10.795 Bias 0.004 0.006 0.006 0.058 0.055 0.054 0.240 0.239 0.240 SD 0.093 0.073 0.065 0.237 0.184 0.156 0.448 0.372 0.324 Table 4.1 Percentiles, Bias, and SD of ˆσ2 nθ1 when d = 2 n 1000 1500 2000 1000 1500 2000 1000 1500 2000 5% 1.933 1.956 1.965 5.453 5.506 5.506 13.152 13.233 13.298 25% 2.017 2.027 2.026 5.738 5.708 5.718 13.817 13.736 13.727 50% 2.085 2.075 2.076 5.925 5.875 5.857 14.297 14.118 14.039 75% 2.147 2.133 2.122 6.120 6.028 5.991 14.788 14.501 14.369 95% 2.250 2.217 2.191 6.392 6.246 6.166 15.498 15.083 14.875 Bias 0.085 0.079 0.076 0.927 0.873 0.852 4.315 4.138 4.055 SD 0.096 0.079 0.070 0.286 0.231 0.196 0.715 0.566 0.483 Table 4.2 Percentiles, Bias, and SD of ˆσ2 nθ1 when d = 3 78 Figure 4.1 Histograms of (cid:112)n/2(ˆσ2 0θ0) − 1) with α = 1.2 when d = 2 and n = 2000. The parameter θ0 is 2, 5 and 10 from left to right. The red curve is the density of N (0, 1). nθ1/(σ2 79 BIBLIOGRAPHY Abramowitz, Milton and Irene A Stegun (1968). Handbook of mathematical functions with formulas, graphs, and mathematical tables. Vol. 55. US Government printing office. Abry, Patrice and Gustavo Didier (2018). Wavelet estimation for operator fractional Brow- nian motion. In: Bernoulli 24 (2). Adler, Robert J (2010). The Geometry of Random Fields. Adler, Robert J and Jonathan E Taylor (2007). Random Fields and Geometry. Amblard, Pierre Olivier and Jean François Coeurjolly (2011). Identification of the mul- tivariate fractional brownian motion. In: IEEE Transactions on Signal Processing 59 (11). Anderes, Ethan (2010). On the consistent separation of scale and variance for Gaussian random fields. In: Annals of Statistics 38 (2). Bachoc, François, Emilio Porcu, Moreno Bevilacqua, Reinhard Furrer, and Tarik Faouzi (2022). Asymptotically equivalent prediction in multivariate geostatistics. In: Bernoulli 28 (4). Bayarr, M. J., James O. Berger, Eliza S. Calder, Keith Dalbey, Simon Lunagomez, Abani K. Patra, E. Bruce Pitman, Elaine T. Spiller, and Robert L. Wolpert (2009). Using statistical and computer models to quantify volcanic hazards. In: Technometrics 51 (4). Bergström, Harald (1952). On some expansions of stable distribution functions. In: Arkiv för matematik 2 (4). Bevilacqua, Moreno, Tarik Faouzi, Reinhard Furrer, and Emilio Porcu (Apr. 2019). Estima- tion and prediction using generalized Wendland covariance functions under fixed domain asymptotics. In: The Annals of Statistics 47 (2). Chan, Grace and Andrew T.A. Wood (2000). Increment-based estimators of fractal dimen- sion for two-dimensional surface data. In: Statistica Sinica 10 (2). Chilès, Jean Paul and Pierre Delfiner (2012). Geostatistics: Modeling Spatial Uncertainty: Second Edition. Cressie, Noel A. C. (Sept. 1993). Statistics for Spatial Data. Wiley. isbn: 9780471002550. Didier, Gustavo and Vladas Pipiras (Feb. 2011). Integral representations and properties of operator fractional Brownian motions. In: Bernoulli 17 (1). 80 Du, Juan, Hao Zhang, and V. S. Mandrekar (Dec. 2009). Fixed-domain asymptotic proper- ties of tapered maximum likelihood estimators. In: The Annals of Statistics 37 (6A). Feldman, Jacob (1958). Equivalence and perpendicularity of Gaussian processes. In: Pacific Journal of Mathematics 8 (4). Feller, William (1991). An introduction to probability theory and its applications, Volume 2. Vol. 81. John Wiley & Sons. Fields, Jerry L. and Mourad El-Houssieny Ismail (May 1975). On the Positivity of some 1F2’s. In: SIAM Journal on Mathematical Analysis 6 (3), pp. 551–559. Gardiner, C.W. (1985). Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences. Proceedings in Life Sciences. Springer-Verlag. isbn: 9783540113577. Gneiting, Tilmann, William Kleiber, and Martin Schlather (2010). Matérn cross-covariance functions for multivariate random fields. In: Journal of the American Statistical Associ- ation 105 (491). Gneiting, Tilmann and Martin Schlather (2004). Stochastic models that separate fractal dimension and the hurst effect. In: SIAM Review 46 (2). Gradshteyn, I. S., I. M. Ryzhik, Alan Jeffrey, Y. V. Geronimus, M. Y. Tseytlin, and Y. C. In: Journal of Biomechanical Fung (1981). Table of Integrals, Series, and Products. Engineering 103 (1). Hanson, D. L. and F. T. Wright (June 1971). A Bound on Tail Probabilities for Quadratic Forms in Independent Random Variables. In: The Annals of Mathematical Statistics 42 (3), pp. 1079–1083. Ibragimov, I. A. and Y. A. Rozanov (1978). Gaussian Random Processes. Springer New York. isbn: 978-1-4612-6277-0. Kaufman, C. G. and B. A. Shaby (June 2013). The role of the range parameter for estimation and prediction in geostatistics. In: Biometrika 100 (2), pp. 473–484. Kent, John T. and Andrew T.A. Wood (1997). Estimating the fractal dimension of a locally self-similar Gaussian process by using increments. In: Journal of the Royal Statistical Society. Series B (Methodological), pp. 679–699. Khoshnevisan, Davar, Yimin Xiao, and Yuquan Zhong (2003). Measuring the range of an additive Lévy process. In: Annals of Probability 31 (2). Khuri, Andre I. (Oct. 2009). Linear Model Methodology. Chapman and Hall/CRC. isbn: 9781420010442. 81 Lahiri, S N (2003). Central Limit Theorems for Weighted Sums of a Spatial Process under a Class of Stochastic and Fixed Designs. Lévy, M. Paul (1940). Le Mouvement Brownien Plan. In: American Journal of Mathematics 62 (1/4). Loh, Wei-Liem (Oct. 2005). Fixed-domain asymptotics for a subclass of Matérn-type Gaus- sian random fields. In: The Annals of Statistics 33 (5). Loh, Wei-Liem (2015). Estimating the smoothness of a Gaussian random field from irreg- In: Annals of Statistics 43 ularly spaced data via higher-order quadratic variations. (6). Loh, Wei-Liem, Saifei Sun, and Jun Wen (Dec. 2021). On fixed-domain asymptotics, param- eter estimation and isotropic Gaussian random fields with Matérn covariance functions. In: The Annals of Statistics 49 (6). Mallat, Stephane (2008). A Wavelet Tour of Signal Processing: The Sparse Way. Mardia, K. V. and R. J. Marshall (Apr. 1984). Maximum Likelihood Estimation of Models for Residual Covariance in Spatial Regression. In: Biometrika 71 (1), p. 135. Mason, J. D. and Yimin Xiao (Jan. 2002). Sample Path Properties of Operator-Self-Similar Gaussian Random Fields. In: Theory of Probability & Its Applications 46 (1), pp. 58–78. Nolan, John (2005). Multivariate stable densities and distribution functions: general and elliptical case Deutsche Bundesbank’s 2005 Annual Fall Conference Multivariate stable densities and distribution functions: general and elliptical case. Paulo, Rui (2005). Default priors for Gaussian processes. Peng, Chien Yu and C. F. Jeff Wu (2014). On the choice of nugget in kriging modeling In: Journal of Computational and Graphical for deterministic computer experiments. Statistics 23 (1). Prudnikov, A.P., Yury Brychkov, and O.I. Marichev (1986). Integrals and series: special functions. Vol. 2. CRC press. Sacks, Jerome, Susannah B. Schiller, and William J. Welch (Feb. 1989). Designs for Com- puter Experiments. In: Technometrics 31 (1), pp. 41–47. Sacks, Jerome, William J. Welch, Toby J. Mitchell, and Henry P. Wynn (1989). Design and analysis of computer experiments. In: Statistical Science 4 (4). Samorodnitsky, Gennady and Murad S. Taqqu (Nov. 2017). Stable Non-Gaussian Random 82 Processes. Routledge. isbn: 9780203738818. Schlather, Martin, Alexander Malinowski, Marco Oesting, Daphne Boecker, Kirstin Strokorb, Sebastian Engelke, Johannes Martini, Felix Ballani, Olga Moreva, Jonas Auel, Peter J Menck, Sebastian Gross, Ulrike Ober, Christoph Berreth, Katharina Burmeister, Ju- liane Manitz, Paulo Ribeiro, Richard Singleton, Ben Pfaff, and R Core Team (2017). RandomFields: Simulation and Analysis of Random Fields. R package version 3.1.50. Stein, Michael (June 1990). Uniform Asymptotic Optimality of Linear Predictions of a Random Field Using an Incorrect Second-Order Structure. In: The Annals of Statistics 18 (2). Stein, Michael L. (Mar. 1988). Asymptotically Efficient Prediction of a Random Field with a Misspecified Covariance Function. In: The Annals of Statistics 16 (1). Stein, Michael L. (Nov. 1989). [Design and Analysis of Computer Experiments]: Comment. In: Statistical Science 4 (4). Stein, Michael L. (Aug. 1993). A simple condition for asymptotic optimality of linear pre- dictions of random fields. In: Statistics & Probability Letters 17 (5), pp. 399–404. Stein, Michael L. (1999a). Interpolation of Spatial Data. Springer New York. isbn: 978-1- 4612-7166-6. Stein, Michael L. (1999b). Predicting random fields with increasing dense observations. In: Annals of Applied Probability 9 (1). Stein, Michael L. (June 2004). Equivalence of Gaussian measures for some nonstationary random fields. In: Journal of Statistical Planning and Inference 123 (1), pp. 1–11. van der Vaart, Aad (Oct. 1996). Maximum likelihood estimation under a spatial sampling scheme. In: The Annals of Statistics 24 (5). Velandia, Daira, François Bachoc, Moreno Bevilacqua, Xavier Gendre, and Jean Michel Loubes (2017). Maximum likelihood estimation for a bivariate Gaussian process under fixed domain asymptotics. In: Electronic Journal of Statistics 11 (2). Wang, Daqing (2010). Fixed domain asymptotics and consistent estimation for Gaussian In: Nat. Univ. random field models in spatial statistics and computer experiments. Singapore PhD thesis, Singapore. Wang, Daqing and Wei-Liem Loh (Jan. 2011). On fixed-domain asymptotics and covariance tapering in Gaussian random field models. In: Electronic Journal of Statistics 5 (none). Ying, Zhiliang (1991). Asymptotic properties of a maximum likelihood estimator with data 83 from a Gaussian process. In: Journal of Multivariate Analysis 36 (2). Ying, Zhiliang (Sept. 1993). Maximum Likelihood Estimation of Parameters under a Spatial Sampling Scheme. In: The Annals of Statistics 21 (3). Zhang, Hao (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. In: Journal of the American Statistical Association 99 (465). Zhang, Hao and Wenxiang Cai (May 2015). When Doesn’t Cokriging Outperform Kriging? In: Statistical Science 30 (2). Zhou, Yuzhen and Yimin Xiao (2018). Joint asymptotics for estimating the fractal indices of bivariate Gaussian processes. In: Journal of Multivariate Analysis 165. 84