Discussion Paper No. 113
MULTIVARIATE CARMA RANDOM FIELDS

Yasumasa Matsuda
and Xin Yuan
April, 2020

Data Science and Service Research
Discussion Paper

Center for Data Science and Service Research
Graduate School of Economic and Management
Tohoku University27-1 Kawauchi, Aobaku
Sendai 980-8576, JAPANMULTIVARIATE CARMA RANDOM FIELDSYASUMASA MATSUDA AND XIN YUAN

Abstract. This paper conducts a multivariate extension of isotropic Lévydriven CARMA random fileds on Rd proposed by Brockwell and Matsuda(2017). Univariate CARMA models are defined as moving averages of a Lévy
sheet with CARMA kernels defined by AR and MA polynomials. We define
multivariate CARMA models by a multivariate extension of CARMA kernels
with matrix valued AR and MA polynomials. For the multivariate CARMA
models, we derive the spectral density functions as explicit parametric functions. Given multivariate irregularly spaced data on R2 , we propose Whittle
estimation of CARMA parameters to minimize Whittle likelihood given withperiodogram matrices and clarify conditions under which consistency and asymptotic normality hold under the so called mixed asymptotics. We finally introduce a method to conduct kriging for irregularly spaced data on R2 by multivariate CARMA random fields with the estimated parameters in a Bayesian
way and demonstrate the empirical properties by tri-variate spatial dataset ofsimulation and of US precipitation data.

1. Introduction
Continuous-time Autoregressive and Moving Average (CARMA) processes havebeen applied as a useful tool to analyze continuous time behaviors in physics andengineers for many years. Ornstein and Uhlenbeck process by Uhlenbeck and Ornstein [22] is such a typical example that is regarded as a CARMA (1, 0) process.
Doob [10] is one of several papers that examined basic properties and statisticalanalysis of CARMA processes. Recently CARMA models have been a resurgence ofinterest by growing needs to analyze high frequency observations in financial timeseries. CARMA modeling can be a tool to connect partial diﬀerential equations
with high frequency data in finance. Statistical properties of CARMA time seriesmodels, including stationary conditions, parametric forms of covariance and spectral density functions, estimation and prediction, have been investigated by many
authors. Brockwell [6] is a good review to see the recent progress in statisticalanalysis by CARMA processes.
Brockwell and Matsuda [7] extended CARMA models for continuous time series
to those for random fields on Rd , d ≥ 1, which we call CARMA random fields. The
main diﬃculty for the random field extension from time series is in no trivial order inspace like that in time series from past to future. Recalling that CARMA processes
are defined by moving averages of Lévy processes specified by CARMA kernels
constructed from AR and MA polynomials with negative real roots, Brockwell andMatsuda [7] employed bilateral moving averages by isotropic CARMA kernels on
Key words and phrases. Bayesian kriging, CARMA random fields, Irregularly spaced data,Lévy sheet. periodogram, Spectral density function, Whittle likelihood.
1

2

YASUMASA MATSUDA AND XIN YUAN

Rd , given byg(s) = g0 (||s||), s ∈ Rd ,where g0 is a time series CARMA kernel and || · || is Euclidean norm. to avoid
the diﬃculty. Several choices of AR and MA polynomials in the CARMA kernels
provide a flexible class of isotropic random fields.
This paper tries a multivariate extension of the univariate CARMA random
fields on Rd . The extension is conducted in a straightforward way of extending
AR and MA polynomials in the univariate models to AR and MA polynomials
with matrix coeﬃcients, which is similar to the extension from discrete ARMA
to vector ARMA models. The extension here is restricted to the cases when AR
polynomial matrices are diagonal in order to obtain explicit parametric expressions
for CARMA kernel matrices. The spectral density matrices are expressed in explicit
parametric forms as a result, which works for parameter estimation on the frequency
domain by the Whittle likelihood, while auto-covariance matrices are not explicitly
obtained except for d = 1, and 3 but are given by Hankel transform of the spectral
densities. It should be notified here that multivariate CARMA random fields target
multivariate irregularly spaced data where observation points for each component
are not supposed to be identical with those for the other components. We develop
a series of procedures to conduct statistical analysis for irregularly spaced databy multivariate CARMA random fields, following the standard way of time series
analysis in Box, Jenkins and Reinsel [4]. Namely, estimation of model parameters
and kriging, which is prediction in spatial data, by the estimated model as if theparameters are known, are developed. We employ frequency domain approach forthe parameter estimation by the Whittle likelihood, while we choose spatial domain
approach for the kriging.
The Whittle estimation is a classical technique in time series analysis which
has been employed by many time series researchers such as Dunsmuir [11], Robinson [21], Hosoya [15] and Dahlhaus [9]. Brockwell and Davis (Chapter 10, 1991)
provides an excellent introduction to the Whittle estimation. This paper employs
the periodogram extended for irregularly spaced data for the Whittle estimation,
which was originally proposed by Matsuda and Yajima [18] and Bandyopadhyay
and Lahiri [1], and has been applied to irregularly spaced data analysis by Bandyopadhyay and Subba Rao [3], Bandyopadhyay, Lahiri, and Nordman [2], Matsuda
and Yajima [19] and Subba Rao [20]. We define the Whittle likelihood function for
multivariate spatial data in a modified form from Matsuda and Yajima [18] to let
it be free from the extra nuisance estimation for distributions of sampling points.
We have established asymptotic normality of the Whittle estimator for multivariate
CARMA random fields which can be non-Gaussian with finite moments of all orders
under the so called mixed asymptotics, which is the asymptotic scheme where sample size and sampling region jointly diverge. The asymptotic results are regarded
as a non-trivial extension of the classical result for discrete stationary time series
by Dunsmuir and Hannan [12] and Dunsmuir [11] to those for continuous random
fields. The nontrivial diﬀerence between them is in the form of asymptotic variancematrix for random field cases which is not separated as those of CARMA kernel andnoise variance unlike time series cases (Remark 4, Dunsmuir [11]). The diﬀerence
comes from a feature of continuous processes for which Kolmogorov formula (seei.e. sec. 5.8 in Brockwell and Davis [5]) does not hold any more.

CARMA RANDOM FIELDS

3

Kriging, which is usually referred to as a minimum mean squared error method of
spatial prediction that depends on the second order properties of spatial processes(Cressie[8]), is one of main purposes in spatial data analysis. Kriging for multivariate spatial data, which is often called as cokriging, is a challenging topic andlots of methods have been proposed in the literatures. Gelfand and Banerjee [13]is a good review for kriging with multivariate spatial process models. Multivariate
CARMA random fields can be regarded as a multivariate extension of the kernel
convolution approach by Higdon [14] in spatial statistics literatures. Our approachre-expresses multivariate spatial observations following CARMA random fields asa form of spatial regression model. Assuming Gaussian for Lévy noise terms driving CARMA, we follow Bayesian approach to conduct kriging by regarding it as a
Bayesian hierarchical model. With the multivariate model involving a large numberof locations, the Bayesian regression requires too heavy computational burden towork for kriging in practice. We apply the technique of Zhang, Sang and Huang[24] to the spatial regression context, which can bee seen also as a kind of covariance tapering by Kaufman, Schervish and Nychka [16]. Specifically, we divide awhole region of spatial observations into several sub-regions to partition the spatial
regression into several sub-models, which lets the posterior computation feasible forlarge multivariate spatial dataset.
Let us start from the definition of multivariate CARMA random fields.
2. Multivariate extension of CARMA random fields
We define multivariate CARMA random fields in a formal way from univariate
ones without introductory arguments. For motivating implications on CARMA
models, see Brockwell [6] or Brockwell and Matsuda [7]. We start from the definitionof multivariate L’evy sheets necessary for the multivariate extension.
2.1. Multivariate Lévy sheet. Define an m-variate Lévy sheet L(x) for m ≥ 1,
which is necessary to multivariate extension of CARMA random fields, by
L(x) = L̇((0, x]), x = (x1 , . . . , xd )′ ∈ Rd+ ,for a random measure L̇, which satisfies(a) if A and B are disjoint Borel sets on Rd , L̇(A) and L̇(B) are independent,and(b) for every Borel set A on Rd with finite Lebesgue measure |A|,
E[exp{iθ′ L̇(A)}] = exp{|A|ψ(θ)}, θ ∈ Rm ,where ψ is the logarithm of the characteristic function of an infinity divisible
distribution.
d
We follow the tradition∫ in writing the integral of a deterministic function g on R
with respect to L̇ as Rd g(x)L(dx).
Let us introduce typical examples.
(a) If ψ(θ) = −θ′ Cθ/2 with a positive definite m×m matrix C, L is a m-variate
Brownian sheet.
(b) If, for a Borel set A on Rd ,L̇(A) =∞∑
i=1

Yi 1xi (A),4

YASUMASA MATSUDA AND XIN YUAN

where xi denotes the location of the ith unit point mass of a Poisson randommeasure on Rd with intensity λ and {Yi } is a sequence of IID random vectors
with distribution function F and independent of {xi }, L is a m-variatecompound Poisson sheet.
We shall restrict attention in this paper to second order Lévy sheets, i.e. thosefor which E[Li (t)2 ] < ∞, i = 1, . . . , m at t = (1, 1, . . . , 1), and then the first andsecond order moments for the sheets are determined by(1)

E{L̇(A)} = µ|A| and var{L̇(A)} = Σ|A|,for a m × 1 vector µ and m × m positive definite matrix Σ.
∑k
If h is a m × m matrix valued function on Rd of the form h(x) = i=1 Ci 1Ai ,where Ai , i = 1, . . . , k are disjoint Borel subsets on Rd with finite Lebesgue measurewith m × m matrix Ci , i = 1, . . . , k,∫h(x)dL(x) :=Rdk∑Ci L̇(Ai ).

i=1

This definition can be extended, by a standard construction, to include all matrixvalued functions h whose components are in L1 (Rd ) ∩ L2 (Rd ).
2.2. Multivariate CARMA random fields. Let us recall the definition of univariate CARMA random fields on Rd by Brockwell and Matsuda [7]. Let L(x) bean univariate Lévy sheet on Rd .
∏p
Definition 1. Let a∗ (z) = z p + a1 z p−1 + · · · ap = i=1 (z − λi ) be a polynomial of
degree p with real coeﬃcients and distinct zeros∏λ1 , . . . , λp having strictly negativeqreal parts and let b∗ (z) = b0 + b1 z + · · · bq z q = i=1 (z − ξi ) with real coeﬃcient bj
and 0 ≤ q < p. Suppose also that λi ̸= µj for all i and j. Then defininga(z) =

p∏(z 2 − λ2i ) and b(z) =

i=1

q
∏

(z 2 − ξi2 ),i=1

the univariate CARMA(p, q) random field driven by a Lévy sheet L is∫
Sd (x) =g(||x − u||)dL(u), x ∈ Rd ,Rdwhere ||x − u|| denotes the Euclidean norm of the vector x − u and g(s) is the
CARMA kernel defined byg(s) =p∑b(λj )a′ (λj )j=1

eλj s , s ∈ R,

where a′ denotes the derivative of the polynomial a.
We shall extend univariate models to m-variate CARMA random fields by extending the scaler functions a(z), b(z) in Definition 1 to the matrix ones with an
m-variate Lévy sheet on Rd described by (1). Here we restrict the cases when a(z)is extended trivially to matrix one by a(z)Im in order to obtain explicit expressions
for CARMA kernel matrix.

CARMA RANDOM FIELDS

5∏p
Definition 2. Let a(z) = j=1 (z 2 − λ2j ) with Re(λj ) < 0 be the polynomial in
Definition 1 and define the matrix polynomial with real m × m matrices B1 , . . . , Bq
byB(z) = z 2q Im + B1 z 2q−2 + · · · + Bq−1 z 2 + Bq ,where B(λi ) ̸= 0 for i = 1, . . . , p. The m-variate CARMA(p, q) random field driven
by a m-variate Lévy sheet L is∫G(||x − u||)dL(u), x ∈ Rd ,Sd (x) =Rd

where G(x) is the m × m CARMA kernel matrix defined byG(s) =p∑1

a′ (λj )j=1

B(λj )eλj s , s ∈ R.

Here we introduce two typical kernels for multivariate CARMA(1,0), CARMA(2,1).
Example. Defining a(z) and B(z) in Definition 2 as (z 2 −λ2 ) and Im , respectively,
we obtain CAR(1) kernel matrix given by1
Im eλs , Re(λ) < 0.
2λDefining a(z) and B(z) in Definition 2 as (z 2 − λ21 )(z 2 − λ22 ) and Im z 2 + B1 ,respectively, we have CARMA(2,1) kernel matrix given by
G(s) =(2)

G(s) =

λ21 Im + B1 λ1 s
λ2 Im + B1 λ2 se
+ 2 2e , Re(λ1 ) < Re(λ2 ) < 0.
22
2λ1 (λ1 − λ2 )2λ2 (λ2 − λ21 )2.3. Second order properties. We show the second order properties of multivariate CARMA(p, q) random fields on Rd , when the driving Lévy sheet has finite
variance. We show in Theorem 1 the spectral density matrix and autocovariancematrix without proof, since it is obtained in a straight forward way from the proof
in Theorem 2 in Brockwell and Matsuda [7]. Let us notice that the spectral density
matrix is explicitly expressed as a function of a, B and Σ, while the auto-covariancematrix is not obtained explicitly except for d = 1 and 3, but given by the Hankel
transform of the spectral density matrix. The explicit form of the spectral density matrix leads us to propose the Whittle likelihood estimation to estimate the
CARMA parameters in the next section.
Theorem 1. If L is a m-variate Lévy sheet with parameters µ and Σ in (1) andif the m-variate CARMA random field is defined as in definition 2, then the mean
vector isp∑1
π d/2 Γ(d + 1)
ESd (x) =B(λi )µ,a′ (λi ) |λ|d Γ(d/2 + 1)
i=1
and the spectral density matrix is(3)

fd (ω) = G̃d (ω)ΣG̃′d (ω), ω ∈ Rd ,whereG̃d (ω) = cd

p
∑

2λi

2
′
2
i=1 a (λi ) (||ω|| + λi )d+12

B(λi ),

6

YASUMASA MATSUDA AND XIN YUAN

with

(
) √
−2d/2−1 Γ d+1
2 ( /) π,−2−d/2 Γ(d)/Γ d2 ,{cd =if d is odd,
if d is even.

The autocovariance matrix isΓd (h) = (2π)d/2 Hd/2−1 Fd (||h||), h ∈ Rd ,where Fd (||ω||) := fd (ω) and Hr denotes the modified Hankel transform of order rdefined by, for Jr , the Bessel function of the first kind of order r,∫ ∞Jr (xy) 2r+1
1
Hr Fd (x) =Fd (y)
ydy, x > 0, r ≥ − .
r
(xy)
2
0
The autocovariance matrix, which is given by the Hankel transform of the spectral density matrix, is explicitly evaluated for d = 1 and d = 3 as{}p∑
exp(z||h||)tΓ1 (h) =B(z)ΣB (z) , h ∈ R
Resz=λi
a(z)2
i=1
and
2π ∑
Γ3 (h) = −Resz=λi||h|| i=1p{}exp(z||h||)t
D(z)ΣD (z) , h ∈ R3 ,za(z)4

respectively, whereD(z) = a′ (z)B(z) − a(z)B ′ (z).
3. Parameter estimation
3.1. Setup. This section focuses on parameter estimation for multivariate CARMA
random fields on R2 . The results here on R2 , however, can be extended to thoseon R3 easily. The point to prevent higher dimensional extension is the rate ofconvergence in Lemma 1, from which the approximation in (11) is not validated ford ≥ 4.
Let X(s) = (X1 (s), . . . , Xm (s))′ , s ∈ R2 be a m-variate random field on R2 ,given by∫
(4)
X(s) =
G(ψ; s − u)dL(u),R2where G is a m × m CARMA(p, q) kernel matrix defined in Definition 2 with theparameter ψ = (λ1 , . . . , λp , B1 , . . . , Bq ), and L is a m-variate Lévy sheet on R2 thatsatisfies (1) with µ = 0. The spectral density function is given in Theorem 1 by
(5)

f2 (θ; ω) = G̃2 (ψ; ω)ΣG̃′2 (ψ; ω), ω ∈ R2 ,for θ = (ψ, Σ) ∈ Rpdim with
1∑
2λi
B(λi ).
2 i=1 a′ (λi ) (||ω||2 + λ2 ) 32
p

G̃2 (ψ; ω) = −i

Let us propose the Whittle estimation for the parameter θ, which is extended
from the one for classical time series analysis, when observation points for eachcomponent of X(s) = (X1 (s), . . . , Xm (s))′ are not supposed to be identical. withthose for the other components.

CARMA RANDOM FIELDS

7

3.2. Whittle Estimation. Suppose we have observed irregularly spaced data that
follow multivariate CARMA random fields (4). Here we allow locations and sample
size of the observations for each component not necessarily to be identical. Namely,
we suppose that pth component has observations Xp (spj ), j = 1, . . . , np on locations
spj that can depend on p, which we assume distributes irregularly over a region in arectangular A = [0, A1 ] × [0, A2 ] on R2 . Let |A| be the area of A, ie., |A| = A1 × A2 .
.
First let us define the periodogram matrix whose (p, q)th element by
Ipq (ω) = dp (ω)dq (ω), ω ∈ R2 ,
where

√
np
′
|A| ∑
dp (ω) =Xp (spj )e−iω spj .
np j=1

For a grid point j = (j1 , j2 ) on Z2 , define a frequency ωj by()′
2πj1 2πj2
ωj =
,
.
A1
A2
For a symmetric compact set D on R2 such that −s ∈ D for s ∈ D, define{}
JD = j = (j1 , j2 ) ∈ Z2 |ωj ∈ D ,and |JD | by the number of elements in JD . The Whittle estimator θ̂ is defined by
the one that minimizes the Whittle likelihood:


{(
}
)
∑
−1
1
lw (θ) = log (6)trf (θ; ωj ) + η K̂I(ωj ) 
m|JD |
j∈JD
w
w
∑
1
w
w
log wf (θ; ωj ) + η K̂ w ,+m|JD |j∈JD

where η is the scaler nuisance parameter that is to be estimated jointly with θ, andK̂ is the matrix to compensate for the bias of the periodogram, which is definedby, for the set of observation points Sp = {spj , j = 1, . . . , np }, p = 1, . . . , m,∑|A|
Xp (sj )Xq (sj ), p, q = 1, . . . , m,K̂pq =
np nq
sj ∈Sp ∩Sq

and defined by 0 if Sp ∩ Sq is null.
Notify the following two points for the proposed likelihood function. First, it ismodified to be scale invariant in the sense that the spectral density matrix multiplied with any constant provides exactly same parameter estimate. In other words,
the variance matrix Σ requires a restriction such as Σ11 = 1 for the identifiability.
Second, the reason why we make the likelihood be scale free is because no-additional
estimation is necessary to let the Whittle estimation be consistent. In other words,
it would be necessary to estimate the quantity related with a density function ofsampling points, if we define the Whittle likelihood with a scale parameter τ by
}
w
w]
)−1
∑ [ {(ww
tr
τ f (θ; ωj ) + η̃ K̂
I(ωj ) + log wτ f (θ; ωj ) + η̃ K̂ w .
j∈JD

8

YASUMASA MATSUDA AND XIN YUAN

See Remark 1 in Matsuda and Yajima [18]. The scale-free likelihood in (6), which
is given by concentrating out the scale parameter, makes it possible to estimate
the parameter θ without any additional estimation at the sacrifice of the scaleparameter in Σ. For a practical discussion, see the beginning of Section 5.

3.3. Asymptotic Results. This section will show the asymptotic results of theWhittle estimator that minimizes (6). First we clarify the scheme under which theasymptotic results shall be derived, which is not trivial unlike time series cases. We
state it as assumption given as C1 below. Under the scheme in C1, we consider theasymptotic results for the Whittle estimator.
C1. The sample size np and the sampling region A = [0, A1 ] × [0, A2 ] diverge
jointly such that A1 → ∞, A2 → ∞, A1 /A2 = O(1) and |A|/np → 0, p =1, . . . , m for the area |A| = A1 × A2 . We shall employ a suﬃx k such as
np = npk , A = Ak to indicate explicitly that they diverge as k tends toinfinity.
C2. Let Sp , p = 1, . . . , m, be the set of sampling points of Xp in A = [0, A1 ] ×[0, A2 ]. We assume that elements in Sp are written as, for p = 1, . . . , m,
spj = (A1 u1,pj , A2 u2,pj ), j = 1, . . . , np ,where upj = (u1,pj , u2,pj ) is a sequence of independent and identically distributed random vectors with a probability density function g(x) supported
on [0, 1]2 which has continuous first derivatives.
C3. X(s), s ∈ R2 follows a random field in (4) driven by a zero-mean Lévy
sheet on R2 with finite moments of all orders. Every component of the
CARMA kernel is bounded and integrable and the spectral density matrix
has continuous second derivatives.
C4. Let Θ be a compact subset in Rpdim and D be a symmetric compact region
on R2 . The parametric spectral density matrix f (θ; ω) defined in (5) is
positive definite and has continuous second derivatives with respect to θ
on Θ × D. θ1 ̸= θ2 implies that f (θ1 ; ω) ̸= f (θ2 ; ω) on a subset of D withpositive Lebesgue measure.
Let us introduce the asymptotic results for the Whittle estimator as k tends to
be ∞ under the asymptotic scheme of C1.
Theorem 2. Under Assumptions C1-C4, the Whittle estimator θ̂k minimizing (6)converges to θ0 in probability as k tends to be infinity.
Theorem 3. If |Ak |3/2 /np → 0 for p = 1, . . . , m. and g(x), x ∈ [0, 1]2(has con)√tinuous second derivatives, in addition with Assumptions C1-C4, |Ak | θ̂k − θ0converges in distribution to()
N 0, bg (Ω1 − Ω2 )−1 (2Ω1 + Π)(Ω1 − Ω2 )−1 ,CARMA RANDOM FIELDS

9

as k tends to be infinity, where, for p, q = 1, . . . , pdim,
{∫
} {∫
}−2
|g(u)|4 dubg = (2π)2([0,1]2|g(u)|2 du,[0,1]2)∂f (θ0 ; ω) −1∂f (θ0 ; ω) −1
f (θ0 ; ω)f (θ0 ; ω) dω,
∂θp
∂θq
D
) (∫
)
(∫
∂ log ||f (θ0 ; ω)||1∂ log ||f (θ0 ; ω)||dωdω ,Ω2,pq =
m|D|
∂θp
∂θq
D
D
m
∑
Πpq =
κabcd ×∫

Ω1,pq =[∫

tra,b,c,d=1∂fG̃′2 (θ0 ; ω)
D

] [∫
]
(θ0 , ω)∂f −1 (θ0 , ω)′G̃2 (θ0 ; ω)dωG̃2 (θ0 ; ω)dω,G̃2 (θ0 ; ω)
∂θp
∂θq
D
ab
cd

−1

where κabcd is the fourth order cumulant of the Lévy sheet given bycum(La (du), Lb (du), Lc (du), Ld (du)) = κabcd du, a, b, c, d = 1, . . . , m.
Remark 1. There are several interesting diﬀerences in the asymptotic variancematrix from the classical one in discrete time series. First, Ω2 vanishes with respectto ψ in θ = (ψ, Σ) in discrete time series cases, because it holds for a spectral density
matrix for discrete time series that
w
w
∫ π
w 1 w
1
w
log ||f (ω)||dω = log w Σw
2π −π
2π w
by Kolmogorov formulra (Theorem 5.8.1, Brockwell and Davis, [5]), while it does
not hold any more in continuous random fields. Second, Π vanishes in discrete
time series except for the components with respect to Σ in θ = (ψ, Σ) (Remark 3,
Dunsmuir, [11]), while it does not in continuous random fields. In other words, theasymptotic variance matrix with respect to ψ and Σ are not separated for continuous
processes.
The twice diﬀerentiable assumption for g(x) in Theorem 3 as well as the diﬀerentiability in Theorem 2 are strict assumptions around edges of sampling points in
A = [0, A1 ] × [0, A2 ], which is rare to be satisfied. To avoid the diﬃculty, let us
propose the tapered periodogram. For a taper h(x), a continuous positive function
on [0, 1]2 , the tapered periodogram is defined by
Ih,pq (ω) = d˜p (ω)d˜q (ω),
for√d˜p (ω) =′|A| ∑Xp (spj )h(spj,1 /A1 , spj,2 /A2 )e−iω spj , p = 1, . . . , m.
np j=1np

The Whittle likelihood in (6) replaced with the tapered periodogram provides Theorems 1 and 2 under relaxed conditions on g(x), namely, the first and second
diﬀerentiability for g(x)h(x), not for g(x).
Let θ̃ be the estimator minimizing the modified Whittle likelihood with a taperh(x), x ∈ [0, 1]2 .

10

YASUMASA MATSUDA AND XIN YUAN

Theorem 4. Under C1-C4, where w(x) = g(x)h(x) has continuous first derivatives
instead of g(x) in C2, the Whittle estimator θ̃ constructed with a taper h(x) is consistent. If we assume more that w(x) = g(x)h(x) has continuous second derivatives
and that |A|3/2 /np → 0 for p = 1, . . . , m, the asymptotic normality in Theorem 3
still holds, where bg in the asymptotic covariance is replaced with
} {∫
}−2
{∫
4
2
|g(x)h(x)| dx|g(x)h(x)| dx.
bh =[0,1]2[0,1]2

4. CARMA Kriging
This section shall propose a way to conduct kriging by multivariate CARMA
random fields when the parameters are known, although they are estimated inpractice by the Whittle estimation we stated in the last section. The kriging is
conducted by following the way of time series forecasting by moving average models,
where the forecast is constructed with the error process recovered by the moving
average recursions. We regard the recovered error term in a Bayesian way as the one
obtained from the posterior distribution when iid assumption for the error terms
is seen as the prior. We extend the time series forecasting procedure to kriging
for multivariate spatial data, where it should be noticed that Gaussian assumptionthat is not necessary for the Whittle estimation is imposed in this section.
4.1. Setup as a spatial regression. We assume that the Lévy sheet driving am-variate CARMA random field is a compound Poisson sheet, which is as a result
given by
Y (s) =∞∑G(s − uj )Zj , s ∈ R2 ,j=1

where {uj }, j = 1, . . ., which are called as knots, are randomly distributed overR2 and Zj follows iid with mean 0 and variance matrix Σ. The restriction tocompound Poisson sheet does not lose generality in terms of covariances. In otherwords, the class of covariances by CARMA models driven by Lévy sheets and of
those by compound Poisson sheets are identical. We suppose the following special
but practical situations under which we shall introduce a method for kriging.
a. We truncate the range of the knots within a compact region C with thenumber of knots M , which follows a Poisson distribution. In addition,
inserting a constant term and an iid measurement error, we employ thefollowing empirically modified CARMA model for kriging by∑(7)
Y (si ) = µ +G(si − uj )Zj + εi , si ∈ R2 ,{uj ,j=1,...,M }⊂C

where we denote the set of knot points as KM = {uj , j = 1, . . . , M } ⊂ C.
b. The parameter Σ is a diagonal matrix without loss of generality, because,2
for the Cholesky decomposition of Σ given by Ldiag(σ12 , . . . , σm
)L′ with thelower triangular matrix L with Lii = 1, i = 1, . . . , m, and the re-defined
G by GL, we obtain the CARMA model driven by the compound Poisson
2
sheet with the diagonal variance matrix Σ = diag(σ12 , . . . , σm).
c. The CARMA kernel G is assumed to be known, although it is in practiceestimated by the Whittle estimation.

CARMA RANDOM FIELDS

112
2
d. The diagonal variance matrices Σ = diag(σ12 , . . . , σm
) and ∆ = diag(δ12 , . . . , δm)
for Zj and εi , respectively, are both unknown parameters to be estimatedin the kriging procedure.
e. We observe samples of Yp (s) at Sp = {spj , j = 1, . . . , np } for p = 1, . . . , m,and aim to krig unknown values of Yp (s) at S̃p = {vpj , j = 1, . . . , ñp } for
p = 1, . . . , m. Denote Vp = Sp ∪ S̃p with Np = np + ñp .
Under the setup, let us express (7) in a componentwise way as

Yp (si ) = µp +m∑∑

Gpq (si − uj )Zjq + εip , si ∈ Vp ,q=1 uj ∈KMwhich we stack together for p = 1, . . . , m, rewriting in a vector form as(8)

Y = Gu Z + E,

where
Y = (Y1′ , . . . , Ym′ )′ for Yp = Yp (si ), si ∈ Vp ,

1N1 0N1 · · · 0N1 Gu,11 Gu,12
 0N2 1N2 · · · 0N2 Gu,21 Gu,22

Gu =  .
..
.
..
..
..
 ..
. ..
.
.
.

···
···
..
.

Gu,1m
Gu,2m
..
.

···

···

Gu,mm

0Nm

0Nm

1Nm

Gu,m1

Gu,m2







for Gu,pq = Gpq (si − uj ), si ∈ Vp , uj ∈ KM ,′ ′Z = (µ′ , Z1′ , . . . , Zm
) for µ = (µ1 , . . . µm )′ and Zq = (Z1q , . . . , ZM q )′ ,E = (ε′1 , . . . , ε′m )′ for εp = (ε1p , . . . , εNp ,p )′ .
4.2. Gibbs sampling. We shall recover Z in (8) in order to construct the krigged
value. We evaluate the posterior of Z by
(9)p(Z, Σ, ∆|Y, u) ∝ p(Y |Z, u, ∆)p(Z|Σ)p(Σ)p(∆),

where we assume Gaussian distributions for the first and second terms in the−2 ′right hand side with the precision matrices given by diag(δ1−2 1′N1 , · · · , δm1Nm ),−2 ′′−2 ′and diag(0m , σ1 1M , · · · , σm 1M ) respectively. Assuming the independent noninformative inverse Gamma distributions as the priors of σi2 and δi2 for i = 1, . . . , m,
we obtain the posterior samples for Z by the Gibbs sampler. Iteratively samplingfrom the posteriors in (9) given knots u simulated independently at each iteration,
we aim to obtain the marginal posterior samples given Y .
0. Initialize Y at S̃p and σp2 , δp2 , p = 1, . . . , m.
1. Simulate knots u = (u1 , . . . , uM ) uniformly distributed over C for M that
follows a Poisson distribution.
2. Update Z for given Y, u, Σ and ∆ by the posterior N (Ωa, Ω) with−1
a = G′u DE
Y,( ′ −1)−1 −1Ω = Gu DE Gu + DZ.

3. Simulate Ỹ for given u, Z by the posterior N (Gu Z, DE ) and update Y withỸ only at S̃p , p = 1, . . . , m.
∑M
2
4. Update σp−2 and δp−2 for given u, Z, E = Y −Gu Z by the posteriors Ga(M/2, j=1 Zjp
/2)
∑Np 2
and Ga(Np /2, j=1 Ejp /2), respectively, for p = 1, . . . , m.

12

YASUMASA MATSUDA AND XIN YUAN

5. return to 1.
As the result of the iterations, we obtain the posterior samples in Step 3 for Y atS̃p , p = 1, . . . , m, by which we can construct interval as well as point estimation as
kriging.
Remark 2. We just simulate the prior samples for u at Step 1 in the Gibbs sampler
unlike Zhang, Sang and Huang [24], since knots are not parameters to generate
CARMA processes but are randomly distributed ones independent of observations.
It means that we construct the krigged values by the marginal posteriors at Step 3
when knots are distributed uniformly over C in (7).
Step 2 requires the inverse of m(M + 1) × m(M + 1) dense matrix, which isinfeasible when M is large. In the empirical example shown later, we will consider
US precipitation data in which CARMA models with M = 6, 000 and m = 3 are
fitted, where the step 2 requires 18, 000 × 18, 000 matrix inversion. In order toavoid the diﬃculty of huge dimensional matrix inversion, we propose a sub-chain
to approximate the inversion in the step 2 with a lower dimensional one.
2-a. We divide the region C randomly into several sub-regions as C1, . . . , Ck .
Following the division, we partition Y, u, Z, Gu , DZ and DE into the ones
corresponding with the devision, denoted as, for i = 1, . . . , k, Y (i) , u(i) ,(i)
Z (i) , Gu , and DZ(i) and DE(i) . We include the constant term µ(i) for all
the divisions Ci , i = 1, . . . , k. In other words, we allow µ to be dependent
on the sub-regions.
2-b. Initialize Z.
2-c. For i = 1, . . . , k, update Z (i) with N (Ω(i) a(i) , Ω(i) ) for
{
}
−1
(i)
a(i) = G′(i)− G(−i)Z (−i) ,u DE(i) Y
u
(
)−1
−1
−1
(i)
Ω(i) = G′(i),u DE(i) Gu + DZ(i)(−i)

where Gu and Z (−i) are the sub-components of Gu and Z excluding theones corresponding with Ci .
Iterating 2-c in the sub-chain for step 2, we obtain the posterior samples of Z by
(m + 1)Mi × (m + 1)Mi matrix inversion for Mi given roughly by M/k.
5. Empirical studies
5.1. tri-variate CARMA(2,1) kernel matrix. This section focuses on tri-variate
CARMA (2,1) random fields on R2 and demonstrates the empirical properties interms of estimation and kriging for simulated and real data. Although it is possible
to apply general class of CARMA(p, q) models, CARMA (2,1) shall be employedwith the two reasons. First, CARMA (2,1) provides wide enough class of covariancestructures to express flexible behaviors both at short and long lags by the CARMA
kernel in (10). The first and second terms in (10) control the behaviors at shortand long lags, respectively. Second, model selection issues from general CARMA
are in a challenging topic currently out of scope in this paper.
We shall employ the CARMA(2,1) kernel in (2) in a modified form. Normalizing
the kernel to satisfy G(0) = Im to guarantee the identifiability of Σ with a newparameter m × m matrix Φ, we obtain
G(s) = Φeλ1 s + (Im − Φ)eλ2 s , Re(λ1 ) < Re(λ2 ) < 0.

CARMA RANDOM FIELDS

true
mean
RMSE

λ1
−3.951−3.9380.325

λ2−0.619−0.648
0.083

ϕ11
0.822
0.825
0.027

ϕ22
0.864
0.858
0.145

ϕ33
0.825
0.805
0.069

ϕ21
1.595
1.547
0.204

ψ21
0.160
0.156
0.033

13
ϕ31
1.017
0.963
0.163

ψ31
0.032
0.030
0.033

ϕ32
0.608
0.694
0.293

ψ32
0.079
0.079
0.096

2
log σ20−0.0020.500

Table 1. The means and root mean squared errors of the Whittle
estimators for tri-variate CARMA (2,1) random field in (10), whichwere evaluated by 100 simulations.

We impose a restriction on Φ and Σ, the variance matrix of Lévy sheet, to increase
the identifiability of CARMA kernel by the Whittle likelihood in (6).
For the Cholesky decomposition for Σ, which is given byΣ = Ldiag(σ12 , σ22 , σ32 )L′ ,for the lower triangular matrix L with Lii = 1, i = 1, 2, 3, modify the CARMA(2,1)
kernel as
G(s)L = ΦLeλ1 s + (I3 − Φ)Leλ2 s .
We restrict the parameter matrix Φ to be within the class of lower triangular, which
improves significantly low identifiability of Φ by the Whittle likelihood. As a result,
the tri-variate CARMA(2,1) kernel we shall employ in this section is expressed as




ϕ11
0
0
1 − ϕ1100
 eλ2 s ,0  eλ1 s +  ψ21
1 − ϕ220(10) G(s) =  ϕ21 ϕ22
ϕ31 ϕ32 ϕ33
ψ31
ψ32
1 − ϕ33
with the diagonal variance matrix Σ for Lévy sheet, which is re-parametrized asΣ = σ12 diag(1, σ22 , σ32 ).
Recalling that the Whittle likelihood in (6) is scale invariant, we notice that theestimable parameters in the CARMA(2,1) model are λ1 , λ2 , ϕij , ψij and σ22 , σ32 .
5.2. Simulation studies. W examined the Whittle estimation for simulated trivariate data by the CARMA (2,1) kernel in (10) on irregularly spaced points overthe compact set A = [0, 50] × [0, 30]. More specifically, we employed the empiricallymodified expression in (7) driven by a Poisson Lévy sheet, where 4,000 knots uniformly distributed over [0, 60] × [0, 60] including A. We designed three independent
sets of 5,000 uniformly distributed points over A as observation points for eachcomponent of tri-variate data. In the notation in the last section, S1 , S2 and S3
were designed to have no intersections. We simulated 100 sets of the tri-variate
data under the setting, where the parameter values for Φ in Table 1 were takenfrom the empirical analysis for US precipitation data shown in the next section.
We estimated the CARMA (2,1) parameters Φ in (10) and σ22 , σ32 in Σ to minimize
the Whitttle likelihood for the 100 sets of simulated data, where the compact region
D in (6) was designed by {ω ∈ R2 , ||ω|| < 2π}. In Table 1, we showed the means
and root mean squared errors evaluated by 100 simulations.
We find from Table 1 that the Whittle estimation works overall for all the parameters in terms of bias and RMSE. The variance parameters has larger RMSE
than the other parameters. The Whittle likelihood evaluates the likelihood over thecompact region within D, namely it ignores the periodograms on higher frequency
regions than D. The ignorance on higher frequencies leads to poor estimation for2
log σ30−0.191
0.372

14

YASUMASA MATSUDA AND XIN YUAN

Figure 1. Weather stations in United States
the variance parameters in Lévy sheets in comparisons with those for the otherparameters in CARMA kernels.
5.3. Real example. This section demonstrates the applications of the tri-variate
CARMA (2,1) random field in (10) to real dataset of US precipitation. Monthly
total precipitation observed at weather stations all over US from 1895 through 1997
are available in the web page of Institute for Mathematics Applied to Geosciences(IMAG):
http://www.image.ucar.edu/Data/US.monthly.met/USmonthlyMet.shtml.
Around 6,000 weather stations, which are scattered over United Stats, are shownin Figure 1.
Regarding the monthly precipitation in November, December, 1996 and January, 1997 as tri-variate spatial observations, we fitted the tri-variate CARMA(2,1)
model to them to examine the estimation and kriging performances. Dividing thewhole dataset into the two sets of in-samples and out-of-samples, we estimated the
CARMA parameters to minimize the Whittle likelihood in (6) by the in-samples
and evaluated the mean squared errors of the kriging constructed with the estimated CARMA model by the out-of-samples. More specifically, the numbers of thedata points in Nov., Dec. and Jan. were 6,841, 6,838 and 6,463, respectively and
we divided each of them into randomly chosen 500 out-of-samples and the rest as
in-samples. The in-sample datasets in Nov. and Dec., in Nov. and Jan. and in
Dec. and in Jan. respectively have intersections of 5,772, 5,400 and 5,415 stations,
respectively. For the Whittle estimation, we took A = [0, 50] × [0, 30] to construct
the periodogram and chose D in (6) as the compact region of {ω ∈ R2 , ||ω|| < 2π}.
At Step 1 in the Bayesian kriging procedure, we designed as knots the points chosenrandomly from among the 6,838 stations in Dec. with M following the Poisson distribution with mean 6,000. To avoid the diﬃculty of huge matrix inversion in StepCARMA RANDOM FIELDS

est.
s.e.

λ1
−3.9510.318

λ2−0.619
0.068

ϕ11
0.822
0.023

ϕ22
0.864
0.020

ϕ33
0.825
0.029

ϕ21
1.595
0.083

ψ21
0.160
0.044

ϕ31
1.017
0.067

15
ψ31
0.032
0.032

ϕ32
0.608
0.027

ψ32
0.079
0.016

2
log σ20.879
0.088

Table 2. The Whittle estimation with the standard error for thetri-variate CARMA (2,1) random field in (10) fitted to US precipitation data in Nov., Dec., 1996 and Jan., 1997.

kerneluni-variate
tri-variatesmoother CARMA(2,1) CARMA(2,1)
Nov.
4.79
6.36
4.32
in-sample
Dec.
15.84
17.67
10.05
7.41
9.22
5.37
Jan.
Nov.
8.83
7.94
4.88
out-of-sample Dec.
26.43
27.94
13.84
Jan.
18.42
15.82
10.12
Table 3. Kriging MSE by the estimated tri-variate CARMA(2,1)
random field in comparisons with those of the benchmarks of kernelsmoother and univariate CARMA(2,1) random field.
MSE2, we employed the sub-chain of 2-a, b and c with randomly chosen 50 sub-regionsin the whole US continent and iterated Step 2-c four times. We iterated 200 times
Steps 1-3 in the kriging procedure and constructed krigigged values by the averagesof the last 100 posterior samples.
The Whittle estimators for the CARMA(2,1) parameters are in Table 2 and theidentified autocorrelation matrix by the formula in Theorem 1 is in Figure 2, while
the mean squared errors of the tri-variate CARMA kriging are shown in Table 3
with those of the Gaussian kernel smoother and univariate CARMA (2,1) model
as benchmarks, where the bandwidth for the kernel smoothing was optimized togive the best kriging performances. The standard errors in Table 2 were evaluatedby 2bg H −1 /|A|, where H is the Hessian matrix of the Whittle likelihood in (6)and bg was evaluated by the kernel density estimation, 4π 2 × 3.09. Notice that it
ignored the asymptotic variance terms of Ω2 and Π in Theorem 3, which would
cause negatively biased approximation.
We summarize the results in the three points. First, the comparison between
Tables 1 and 2 demonstrates the negatively biased standard errors in Table 2.
Namely, the standard errors obtained via the inverse Hessian are smaller than theones via the simulations, following the expectation from the asymptotic variancematrix in Theorem 3. Second, the univariate CARMA models and kernel smoothing
provide the close kriging MSEs. This is because kriging by uni-variate CARMA
models is regarded as a kind of smoothing with the kernel and bandwidth specified
by CARMA kernels. Finally, it is found from Figure 2 that correlation of Dec. andJan. is more durable than those of the other two especially in short lag distances,which accounts for the significantly better performances of the tri-variate CARMA
kriging over the univariate benchmarks. A local shock at a point that the univariate
benchmarks cannot account can be caught by the tri-variate CARMA models via
the correlations between Dec. and Jan. or Dec. and Nov., which resulted in thesignificant kriging improvements.

2
log σ3−0.9030.114

16

YASUMASA MATSUDA AND XIN YUAN

Figure 2. Autocorrelations as a function of lag distance identified
by the tri-variate CARMA(2,1) random field for US precipitationin Nov., Dec., 1996 and Jan., 1997.

6. Conclusion
This paper conducts a multivariate extension of univariate CARMA random
fields proposed by Brockwell and Matsuda [7], which is obtained as a continuous
analogue of discrete time moving average models. The features of the extension
are summarized in the following points. First, the extension is simply designed toprovide with explicit parametric expressions of multivariate CARMA kernel matrix and as a result with the spectral density matrix. Second, we propose theWhittle likelihood to estimate the CARMA parameters eﬃciently. Unlike usual
normal likelihood in spatial domains that requires huge dimensional matrix inversion of covariances, the Whittle likelihood is eﬃciently evaluated by the matrix
inversion of spectral density matrices. Third, we successfully derived the asymptotic normality of the Whittle estimation without imposing Gaussian assumption
but with the existence of all finite order moments. The asymptotic variance matrix is similar but diﬀerent from the one in traditional discrete time series analysis.
The diﬀerence comes from the feature of continuous random fields on which thecelebrated Kolmogorov formula does not hold any more. Fourth, we propose a
Bayesian way of kriging by multivariate CARMA random fields under Gaussianconditions. Although it requires huge matrix inversion, we give a way to avoid
the diﬃculty to provide kriging eﬃciently. Finally, our proposed methodology ofmultivariate CARMA random fields works well in practice. Multivariate CARMA
model captures well durable spatial correlations among components of multivariate
observations to provide better kriging than those by univariate modeling.
The summarized features are all related with second order properties of CARMA
random fields, which are resulted from the assumption that driving Lévy sheet has

CARMA RANDOM FIELDS

17

finite variance matrix. Our future study is to examine CARMA models driven byinfinite variance Lévy sheets, which would open fruitful applications in theories andpractices for random fields.
7. Proofs
Proof of Theorem 2. Let θ1 ∈ Θ be a parameter that is not equal to θ0 . By Lemmas
1 and 2,
[
]
∫
∫
{
}
1
1
tr f −1 (θ1 ; ω)f (θ0 , ω) dω +log ∥f (θ1 ; ω)∥ dω + log τ0 ,lw (θ1 ) → log
m|D| Dm|D| D:= l∞ (θ1 ),say, in probability as k tends to be infinity, where τ0 is defined in Lemma 2. Hence,by Jensen’s inequality,l∞ (θ1 ) − l∞ (θ0 ) =
[
]
∫
∫
w
w
{ −1
}
1
1
log
tr f (θ1 ; ω)f (θ0 , ω) dω −log wf −1 (θ1 ; ω)f (θ0 ; ω)w dω
m|D| Dm|D| D> 0.
It follows that, for any positive constant L(θ0 , θ1 ) that is smaller than l∞ (θ1 ) −l∞ (θ0 ),lim P (lw (θ0 ) − lw (θ1 ) < −L(θ0 , θ1 )) = 1.

k→∞

For any δ > 0, there exists an Hk,δ of the form,

−1 

 C ∑
  C2 ∑1δtr (I(ωj ))tr (I(ωj )) + C3
 |JD |
  |JD |

j∈JD

j∈JD

such that, for any θ1 and θ2 that satisfy |θ2 − θ1 | < δ,
|lw (θ2 ) − lw (θ1 )| < Hk,δ ,because of the non-negative definiteness of K̂ and of the mean value theorem. It isseen from the form of Hk,δ that there exists a δ > 0 such thatlim P (Hk,δ < K(θ0 , θ1 )) = 1.

k→∞

Applying Lemma 2 of Walker [23], we have the consistency.
Proof of Theorem 3. By Taylor series expansion,0=∂lw (θ0 ) ∂ 2 lw (θ∗ )∂lw (θ̂)=+(θ̂ − θ0 ),
∂θ
∂θ
∂θ∂θ′

where θ∗ is the mean value between θ0 and θ̂. Hence
}−1
{
}
) {
√ (
√
∂ 2 lw (θ∗ )∂lw (θ0 )|A| θ̂ − θ0 = m|D|
|A| −m|D|
.
∂θ∂θ′
∂θ

□

18

YASUMASA MATSUDA AND XIN YUAN

The (p, q)th element of the first factor, which is the Hessian matrix, is evaluatedas, for p, q = 1, . . . , pdim,][ −1 ∗|D| ∑∂f (θ ; ωj ) ∂f (θ∗ ; ωj )
−
tr
|JD |
∂θp
∂θq
j∈JD
[
]
[ −1 ∗]
1 |D| ∑
∂f −1 (θ∗ ; ωj ) I(ωj )|D| ∑∂f (θ ; ωj ) I(ωj )
tr
×
tr
−
m|D| |JD |∂θpτ̂ (θ∗ )|JD |∂θqτ̂ (θ∗ )
j∈JD

j∈JD

+ op (1),whereτ̂ (θ) =∑ []1
tr f −1 (θ; ωj )I(ωj ) .
m|JD |j∈JD∗

Noting that θ converges to θ0 , we find that the Hessian converges in probability
to Ω1 − Ω2 by Lemmas 1 and 2.
The pth element of the second factor, which is the score vector, is evaluated as,
for p = 1, . . . , pdim,
√
[{
} −1
]
|A||D| ∑I(ωj )∂f (θ0 ; ωj )tr
− f (θ0 ; ωj )+ op (1),|JD |τ̂ (θ0 )∂θpj∈JD

which is, by Lemma 1, equal to√[{}]|A||D| ∑I(ωj ) EI(ωj ) ∂f −1 (θ0 ; ωj )
tr
(11)
−
+ op (1).
|JD |
τ0
τ0
∂θp
j∈JD

By applying Lemma 3 to the first term, it is equivalent to show the asymptotic
distribution of
[{
}
]
√ ∫
I(ω) EI(ω) ∂f −1 (θ0 ; ω)Jp = |A|tr−dω, p = 1, . . . , pdim.
τ0
τ0
∂θp
D
For simplicity, we define∂f −1 (θ0 ; ω)
,
∂θp
∫
′
ϕ̃p (s) =ϕp (ω)e−iω s dω,ϕp (ω) =D

and re-express the random term for Jp as
(12)

Tp =

nb
na ∑
m ∑
m
∑
|A|3/2 ∑a=1 b=1

na nb c=1

τ0−1 Xa (sac )Xb (sbd )ϕ̃ba,p (sac − sbd ).

d=1

We shall show in Lemmas 4 and 5 that
Cov (Tp , Tq ) → bg (2Ω1,pq + Πpq ), p, q = 1, . . . , pdim,cum (Tp1 , . . . , Tpr ) → 0, p1 , . . . , pr = 1, . . . , pdim, for r ≥ 3,
as k tends to ∞, which proves the asymptotic normality in Theorem 3.

□

CARMA RANDOM FIELDS

19

8. Lemmas
Lemma 1.

(EIab (ω) = τ0 fab (ω) + O)
nab |A|−2−2, a, b = 1, . . . , m,
+ A1 + A2
na nb

where na , nb and nab are the number of elements in Sa , Sb and Sa ∩ Sb , respectively,
and∫|g(x)|2 dx.
τ0 = (2π)2[0,1]2

Proof. We find from Lemma 3 in Matsuda and Yajima [18] that the expectation isevaluated as
∑
(
)
|A|
−2
τ0 fab (ω) +EXa (sj )Xb (sj ) + O A−2,
1 + A2na nb
sj ∈Sa ∩Sb□which completes the proof.
Lemma 2. For a square integrable function ψ(ω), ω ∈ R2 ,



√
∑
|D|
Iab (ωj )ψ(ωj ) = O (1) , a, b = 1, . . . , m.
var
|A|


|JD |
j∈JD

Proof. Letψ̂(s) =′
|D| ∑
ψ(ωj )e−iωj s , s ∈ A = [0, A1 ] × [0, A2 ],
|JD |

j∈JD

which is extended periodically to [−A1 , A1 ] × [−A2 , A2 ]. Then the object for thevariance is evaluated asnbna ∑|A|3/2 ∑Xa (sac )Xb (sbd )ψ̂(sac − sbd ).
na nb c=1
d=1

The variance is given by(|A|3 ∑ ∑ ∑ ∑cum (Xa (sac1 ), Xb (sbd1 ), Xa (sac2 ), Xb (sbd2 ))E 2 2
na nb cd1 c2 d21)
+ γaa (sac1 − sac2 )γbb (sbd1 − sbd2 ) + γab (sac1 − sbd2 )γba (sbd1 − sac2 )× ψ̂(sac1 − sbd1 )ψ̂(sac2 − sbd2 )∫ ∫ ∫ ∫ (=cum(Xa (u1 ), Xb (v1 ), Xa (u2 ), Xb (v2 ))
A

A

A

A

)
+ γaa (u1 − u2 )γbb (v1 − v2 ) + γab (u1 − v2 )γba (v1 − u2 )× δ(u1 − v1 )δ(u2 − v2 )|A|−1 g(u1 /A)g(v1 /A)g(u2 /A)g(v2 /A)du1 dv1 du2 dv2 + o(1).
The first term is, by expressing the cumulant term with the cumulant spectrum:m∑(13) fabab (ω1 , ω2 , ω3 ) =κef gh G̃ae (ω1 )G̃bf (ω2 )G̃ag (ω3 )G̃bh (ω1 + ω2 + ω3 ),e,f,g,h=1

20

YASUMASA MATSUDA AND XIN YUAN

given by
∫ ∫ ∫

′

′

′

fabab (ω1 , ω2 , ω3 )eiω1 (u1 −v2 ) eiω2 (v1 −v2 ) eiω3 (u2 −v2 ) dω1 dω2 dω3∫ ∫ ∫ ∫
|A|−1
δ(u1 − v1 )δ(u2 − v2 )g(u1 /A)g(v1 /A)g(u2 )/Ag(v2 /A)du1 dv1 du2 dv2 ,
R2

R2

R2

A

A

A

A

which is, by Schwarz inequality, bounded by√∫ ∫ ∫
2∏(14)|fabab (ω1 , ω2 , ω3 )|Pj dω1 dω2 dω3 ,
R2

j=1

for
P1 = |A|

−1

R2

∫ ∫

′

2

′

ψ̂(u1 − v1 )g(u1 /A)g(v1 /A)eiω1 u1 eiω2 v1 du1 dv1 ,A

P2 = |A|−1

R2

∫ ∫

A

A

A

′

′

2

ψ̂(u2 − v2 )g(u2 /A)g(v2 /A)e−i(ω1 +ω2 ) v2 eiω3 (u2 −v2 ) du2 dv2 .

By applying Perseval’s equality to both terms, (14) is bounded by∫ A1 ∫ A2
2Cψ̂(s) ds,
(15)
−A1

−A2

which is evaluated as
4C|A||D|2 ∑|ψ(ωj )|2 < C ′|JD |2j∈JD∫|ψ(ω)|2 dω = O(1).
D

Also the second and third terms in the variance are bounded by a constant withthe same argument, which completes the proof.
□
Lemma 3. For a square integrable function ψ(ω), ω ∈ R2 ,√√ ∫|A||D| ∑Iab (ωj )ψ(ωj ) − |A|
Iab (ω)ψ(ω)dω = op (1).
|JD |
D
j∈JD

Proof. Let′|D| ∑ψ(ωj )e−iωj s , s ∈ A = [0, A1 ] × [0, A2 ],
|JD |
j∈JD
∫
′
ψ̃(s) =
ψ(ω)e−iω s dω, s ∈ R2 ,ψ̂(s) =

D

the first one of which is extended periodically to [−A1 , A1 ] × [−A2 , A2 ]. For δ(s) =ψ̂(s) − ψ̃(s), the diﬀerence is evaluated asa ∑b|A|3/2 ∑Xa (sac )Xb (sbd )δ(sac − sbd ).
na nb c=1

n

n

d=1

The variance, which is similarly evaluated till (15) in Lemma 2, is bounded by∫ A1 ∫ A2
2(16)C|δ(s)| ds.
−A1

−A2

CARMA RANDOM FIELDS

21

Notice that ψ̂(s), ψ̃(s) are square integrable, since∫ A1 ∫ A2∫∑|ψ̂(s)|2 ds = |A||D|2 |JD |−2|ψ(ωj )|2 < C|ψ(ω)|2 dω < ∞,
0

0

∫
R2

|ψ̃(s)|2 ds = (2π)2

D

j∈JD

∫

|ψ(ω)|2 dω < ∞.
D

It follows that, for any ε > 0, there exists a compact set BM = [−M1 , M1 ] ×[−M2 , M2 ] ⊂ [−A1 , A1 ] × [−A2 , A2 ] such that
∫ A1 ∫ A22
ψ̂(s) − ψ̂(s)IBM (s) ds < ε,−A1 −A2∫2ψ̃(s) − ψ̃(s)IBM (s) ds < ε.
R2Then (16) is bounded by{∫∫A1 ∫ A2
2
Cψ̂(s) − ψ̂(s)IBM (s) ds +
−A1

−A2

2

∫

ψ̂(s) − ψ̃(s) ds +

BM

R2

2

ψ̃(s)IBM (s) − ψ̃(s) ds< C ′ ε,□which completes the proof.
Lemma 4.
Cov(Tp , Tq ) → bg (2Ω1,pq + Πpq ), p, q = 1, . . . , pdim.
Proof. First, ϕ̃p (s) in Tp defined in (12) may be replaced withϕ̃M
p (s) = ϕ̃p (s)IBM (s), p = 1, . . . , pdim,

for a suﬃciently large compact set BM = [−M1 , M1 ] × [−M2 , M2 ], since the variance of the diﬀerence between them may be made arbitrary small by following theargument till (15) in Lemma 2. The covariance for the one replaced with ϕ̃Mp (s) isevaluated as∫ ∫ ∫ ∫
m
m
∑
∑
|A|−1
τ0−2 {cum(Xa1 (u1 ), Xb1 (v1 ), Xa2 (u2 ), Xb2 (v2 ))A

a1 ,a2 =1 b1 ,b2 =1

A

A

A

+ γa1 a2 (u1 − u2 )γb1 b2 (v1 − v2 ) + γa1 b2 (u1 − v2 )γb1 a2 (v1 − u2 )}M× ϕ̃Mp,b1 a1 (u1 − v1 )ϕ̃q,b2 a2 (u2 − v2 )g(u1 /A)g(v1 /A)g(u2 /A)g(v2/A)du1 dv1 du2 dv2+ o(1).
The first term is, by using the cumulant spectrum defined in (13), evaluated as∑ ∑∫ ∫ ∫fa1 b1 a2 b2 (ω1 , ω2 , ω3 )|A|−1 τ0−2 P1 P2 dω1 dω2 dω3
a1 ,a2 b1 ,b2

for

∫ ∫

R2

R2

′

P1 =
∫A ∫A
P2 =
A

R2

A

′

iω1 u1 iω2 v1
ϕ̃Meg(u1 /A)g(v1 /A)du1 dv1 ,p,b1 a1 (u1 − v1 )e′′}iω3 u2 −i(ω1 +ω2 +ω3 ) v2
ϕ̃Meg(u2 /A)g(v2 /A)du2 dv2 .
q,b2 a2 (u2 − v2 )e22

YASUMASA MATSUDA AND XIN YUAN

By change of variables by u1 − v1 = l1 , u2 − v2 = l2 and the compactness of the
supports of ϕ̃M , |A|−1 P1 P2 is evaluated asM
ϕMp,b1 a1 (ω1 )ϕq,b2 a2 (ω3 )F (ω1 + ω2 ) + o(1),

where
−1

∫

′

F (ω) = |A|−2
ϕM
p (ω) = (2π)

2g 2 (u/A)eiω u du ,
∫A

′

iω s
ϕ̃M
ds.
p (s)eBM

It follows by Lemma 1(c) in Matsuda and Yajima [18] that the first term converges
to
∫ ∫
m
m
∑
∑
M
bg
fa1 b1 a2 b2 (ω1 , −ω1 , ω3 )ϕM
p,b1 a1 (ω1 )ϕq,b2 a2 (ω3 )dω1 dω3 .
R2

a1 ,a2 =1 b1 ,b2 =1

R2

Replacing the cumulant spectrum with (13), we have the result arbitrary close tobg Πpq by taking BM large. By the same argument, the second and third terms
converge to bg Ω1,pq , which completes the proof.
□
Lemma 5. For r ≥ 3,()
cum (Tp1 , . . . , Tpr ) = O |A|−r/2+1 , p1 , . . . , pr = 1, . . . , pdim.

Proof. First, ϕ̃p (s) in Tp defined in (12) may be replaced withϕ̃M
p (s) = ϕ̃p (s)IBM (s), p = 1, . . . , pdim,
for a suﬃciently large compact set BM = [−M1 , M1 ] × [−M2 , M2 ], since the variance of the diﬀerence between them may be made arbitrary small by following theargument till (15) in Lemma 2.
The cumulant for the one replaced with ϕ̃Mp (s) is evaluated ascum(Tp1 , . . . , Tpr ) =
∫
×
×

···

a1 ,b1 =1

∫
···

A
r
∏

m
∑

m
∑

|A|−r/2

ar ,br =1

cum(Xa1 (u1 )Xb1 (v1 ), . . . , Xar (ur )Xbr (vr ))Aϕ̃M
pj ,bj aj (uj − vj )g(uj /A)g(vj /A)duj dvj + o(1).

j=1

Let Yi = (Yi1 , Yi2 ) be the two dimensional random vector defined by (Xai (ui ), Xbi (vi ))for i = 1, . . . , r. The cumulant term in the equation is expressed by the formula inLeonov and Shiryaev [17] asq∑ ∏cum(Y (Dp )),∪qp=1 Dp p=1

where the summation is taken over all the indecomposable partition ∪qp=1 Dp of thetwo way table of indices:CARMA RANDOM FIELDS

23(1, 1) (1, 2)
..
.. .
.
.
(r, 1) (r, 2)
Let us prove only for cum(Y11 , Y12 , . . . , Yr1 , Yr2 ), the highest order term, since theother cases are similarly evaluated. We express the highest cumulant term by the
2rth order cumulant spectrum ascum(Xa1 (u1 ), Xb1 (v1 ), . . . , Xar (ur ), Xbr (vr ))
∫
∫
r
r−1
∏
∏
′
′
=
···
fa1 b1 ···ar br (ω1 , . . . , ω2r−1 )eiωj (uj −vr )eiωj (vj −vr ) dω1 · · · dω2r−1 ,
R2

R2

j=1

j=1

where
m
∑

(17) fa1 b1 ···ar br (ω1 , . . . , ω2r−1 ) =κe1 ,...,e2r

e1 ,...,e2r =1

× G̃br ,e2r (ω1 + . . . + ω2r−1 )r∏G̃aj ,e2j−1 (ω2j−1 )

j=1

r−1
∏

G̃bj ,e2j (ω2j ).

j=1

By replacing the highest cumulant term with the spectrum expression, the summand of the corresponding term is given by|A|−r/2
×

∫

r−1
∏∫

A

R2

∫

A

j=1

∫ ∫

∫
···

R2

fa1 b1 ···ar br (ω1 , . . . , ω2r−1 )dω1 · · · ω2r−1′

A′iω2j−1 uj iω2j vj
ϕ̃Meg(uj /A)g(vj /A)duj dvj
pj ,bj aj (uj − vj )e′A′iω2r−1 ur −i(ω1 +···+ω2r−1 ) vrϕ̃Meg(ur /A)g(vr /A)dur dvr .
pr ,br ar (ur − vr )ewhich is, by changes of variables by uj − vj = lj , j = 1, . . . , r and the compactnessof the supports of ϕ̃M , evaluated as|A|−r/2

∫

∫
R2

···

R2

fa1 b1 ···ar br (ω1 , . . . , ω2r−1 )dω1 · · · ω2r−1

× ϕM
pr ,br ar (ω2r−1 )D(−ω1 − · · · − ω2r−2 )

r−1
∏

ϕM
pj ,bj aj (ω2j−1 )D(ω2j−1 + ω2j ) + o(1),

j=1

where
∫

′

g 2 (u/A)eiω u du.

D(ω) =A

24

YASUMASA MATSUDA AND XIN YUAN

Again by change of variables of ω2j−1 + ω2j = λj , j = 1, . . . , r − 1, and by replacing
the spectrum with (17), it is evaluated as∫∫m∑|A|−r/2
κe1 ,...,e2r···
G̃br ,e2r (λ1 + . . . + λr−1 + ω2r−1 )R2

e1 ,...,e2r =1

×r∏

G̃aj ,e2j−1 (ω2j−1 )

j=1

r−1
∏

R2

G̃bj ,e2j (λj − ω2j−1 )dω1 dλ1 · · · dω2r−3 dλr−1 dω2r−1

j=1

× ϕM
pr ,br ar (ω2r−1 )D(−λ1 − · · · − λr−1 )

r−1
∏

ϕM
pj ,bj aj (ω2j−1 )D(λj ),j=1

whose summand is bounded by
r ∫
r−1
∏
∏∫
C|A|−r/2+1
G̃aj ,e2j−1 (ω2j−1 )dω2j−1
× |A|−1

∫

j=1

R2

R2

j=2

R2

G̃bj ,e2j (λj − ω2j−1 )D(λj )dλjG̃br ,e2r (λ1 + . . . + λr−1 + ω2r−1 )G̃b1 ,e2 (λ1 − ω1 )D(λ1 )D(−λ1 − · · · − λr−1 )dλ1 ,which we find to be O(|A|−r/2+1 ) by Lemma 2 in Matsuda and Yajima [18].

□Acknowledgments We sincerely thank Professor Peter Brockwell of Colorado
State University for his kind advice to the derivation of the auto-covariance matrices
in Theorem 1.
References[1] Bandyopadhyay, S. and Lahiri, S. (2009). Asymptotic properties of discrete fourier transforms
for spatial data. Sankhya, Series A, 71, 221259.
[2] Bandyopadhyay, S., Lahiri, S. N. and Nordman, D. (2015). A frequency domain empirical
likelihood method for irregular spaced spatial data. Annals of Statistics, 43, 2. 519545.
[3] Bandyopadhyay, S., and Subba Rao, S. (2015). A test for stationarity of irregular spaced
spatial data. Journal of the Royal Statistical Society:Series B, 79, 95-123.
[4] Box, G. E. P., Jenkins, G. M. and Reinsel, G. C. (1994). Time Series Analysis: Forecasting
and Control. Englewood Cliﬀs: Prentice Hall.
[5] Brockwell, P. J. and Davis, R. A. (1991) Time Series: Theory and Methods, 2nd edn. SpringerVerlag.
[6] Brockwell, P. J. (2014). Recent results in the theory and applications of CARMA processes,
Ann. Inst. Stat. Math., 66, 647-685.
[7] Brockwell, P. J. and Matsuda, Y. (2017). Continuous auto-regressive moving average random
fields on Rn . Journal of the Royal Statistical Society: Series B, 79(3), 833-857.
[8] Cressie, N. (1993). Statistics for Spatial Data. Wiley, New York.
[9] Dahlhaus, R. (1997). Fitting time series models to nonstationary processes. Annals of Statistics, 16. 137.
[10] Doob, J. L. (1944). The elementary Gaussian processes. Annals of Mathematical Statistics,
25, 229-282.
[11] Dunsmuir, W. (1979). A central limit theorem for parameter estimation in stationary vector
time series and its application to models for a signal observed with noise. Annals of Statistics,
7, 490-506.
[12] Dunsmuir, W. and Hannan, E. J. (1976). Vector linear time series models. Adv. Appl. Probability, 8, 339-364.
[13] Gelfand, A. E. and Banerjee, S. (2010). Multivariate spatial process models. In: Gelfand, A.
E., Diggle, P. J., Fuentes, M. and Guttorp, P. (eds) Handbook of Spatial Statistics, CRC Press.
493-515.

CARMA RANDOM FIELDS

25[14] Higdon D. (2002). Space and Space-Time Modeling using Process Convolutions. In: Anderson
C.W., Barnett V., Chatwin P.C., El-Shaarawi A.H. (eds) Quantitative Methods for Current
Environmental Issues. Springer, London.
[15] Hosoya, Y. (1997). A limit theory for long-range dependence and statistical inference on
related models Annals of Statistics, .25, 1, 105-37.
[16] Kaufman, C., Schervish, M., Nychka, D. (2008). Covariance tapering for likelihood-based
estimation in large spatial data sets. J. Amer. Statist. Assoc. 103, 484, 1545-1555.
[17] Leonov, V. P. and Shiryaev, A. N. (1959). On a method of calculation of semi-invariants.
Theory of Probability and its Applications, 4(3), 319-329.
[18] Matsuda, Y. and Yajima, Y. (2009). Fourier analysis of irregularly spaced data on Rd. Journalof the Royal Statistical Society: Series B, 71(1), 191-217.
[19] Matsuda, Y. and Yajima, Y. (2018). Locally stationary spatio-temporal processes. Jpn. J.
Stat. Data Sci., 1, 41-57.
[20] Subba Rao, S. (2018). Statistical inference for spatial statistics defined in the Fourier domain.
Annals of Statistics, 46, 469-499.
[21] Robinson, P. M. (1995). Gaussian semiparametric estimation of long range dependence. Annals of Statistics, 23, 1630-1661.
[22] Uhlenbeck, G. E. and Ornstein, L. S. (1930). On the theory of Brownian Motion, PhysicalReview, 36, 823-841.
[23] Walker, A. M. (1964). Asymptotic properties of least-squares estimates of parameters of thespectrum of a stationary non-deterministic time-series. Journal of the Australian MathematicalSociety, 4, 363-384.
[24] Zhang, B., Sang, H. and Huang, J. Z. (2013). Full-scale approximations of spatio-temporal
covariance models for large datasets. Statist. Sin., 25, 99-114.
Graduate School of Economics and Management, Tohoku University, 27-1 Kawauchi,
Aoba ward, Sendai 980-8576, Japan
E-mail address: yasumasa.matsuda.a4@tohoku.ac.jp