Discussion Paper No. 46

Lévy-driven CARMA Random
Fields on RⁿPeter J. Brockwell

Yasumasa Matsuda

July 2015

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Lévy-driven CARMA Random Fields on Rn
Peter J. Brockwell∗

Yasumasa Matsuda†

Abstract
We define an isotropic Lévy-driven CARMA(p, q) random field on Rn as the integralof an isotropic CARMA kernel with respect to a Lévy sheet. Such fields constitutea parametric family characterized by an autoregressive polynomial a and a moving
average polynomial b having zeros in both the left and right complex half-planes.
They extend the well-balanced Ornstein-Uhlenbeck process of Schnurr and Woerner(2011) to a well-balanced CARMA process in one dimension (with a much richer
class of autocovariance functions) and to an isotropic CARMA random field on
Rn for n > 1. We derive second-order properties of these random fields and findthat CAR(1) constitutes a subclass of the well known Matérn class. If the drivingLévy sheet is compound Poisson it is a trivial matter to simulate the corresponding
random field on any n-dimensional hypercube. Joint estimation of CARMA kernelparameters and knots locations is proposed in cases driven by compound Poisson
sheets and is illustrated by applications to land price data in Tokyo as well as
simulated data.

Keywords: compound Poisson, convolution, CARMA random field, Gibbs sampling, knot
selection, Lévy sheet, Lévy noise, Matérn class.

1

Introduction

This paper introduces a class of scalar random fields Sn (t), t ∈ Rn , which, when n = 1,
reduces to a class of non-causal Lévy-driven CARMA processes whose properties havebeen studied by Brockwell and Lindner (2009).
Traditionally the modelling and analysis of spatial data involves the fitting of a regression model with spatially correlated errors, specified in terms of a parametric family∗

Statistics Department, Colorado State University, Fort Collins, Colorado 80523-1877, U.S.A.
pjbrock@stat.colostate.edu, Tel: 970 412 3861, Fax: 970 491 7895.
†
Graduate School of Economics and Management, Tohoku University, Sendai 980-8576, Japan.
matsuda@econ.tohoku.ac.jp

1

of covariance functions. A large number of families of such covariance functions have been
proposed in the literature, in many of which the covariance between the field at any two locations depends only on the Euclidean distance between them. The most frequentlly usedfamily is the Matérn class (Matérn, 1960), for which both the spectral density function
and the covariance function can be written explicitly (see e.g. Stein, 1999) . The Matérnclass was extended to a class of nonstationary and anisotropic covariances by Paciorek
and Schervish (2006).
When fitting covariance models to spatial data, the large sizes of the data sets give rise
to significant computational diﬃculties which mean in particular that maximum likelihood
estimation is usually not feasible.
To overcome some of the diﬃculties Higdon (2002) suggested specifying the field Sn
directly as a convolution of a kernel function with a continuous-parameter white noise
field. In this paper we utilize this approach, using the white-noise field associated witha second-order Lévy sheet L and a family of kernels gn which belong to a conveniently
parameterized class of kernels reducing, when n = 1, to a class of non-causal continuoustime ARMA kernels. The field can be expressed as∫
Sn (t) =gn (t − u)dL(u), for t ∈ Rn
Rn

where gn is specified in Definition 3.1. The resulting covariances between Sn (t) and Sn (t′ )depend only on the displacement vector d = t − t′ . Moreover if gn is isotropic (as it is in
this paper) then so is Sn , and the covariances will depend only on the magnitude of d. An
especially convenient special case occurs when L is a compound Poisson random measureon the Borel subsets of Rn . In this case Sn can be expressed as
Sn (t) :=∞∑gn (t − xi )Yi ,i=1

where xi denotes the location of the ith unit point mass of a Poisson random measure on
Rn and the sequence {Yi } is a sequence of i.i.d. random variables independent of {xi }.
We shall refer to Sn in this case as a compound Poisson random field with kernel gn . The
random points xi will be referred to as knots. The knots contained in any bounded subset
D ⊂ Rn , given their number, are uniformly and independently distributed on D.
As pointed out by Higdon (2002), the convolution model permits dimension reductionin the following sense. If we consider only those contributions to the field from a finite
number, M , of knots contained in a bounded subset of Rn , then the necessary computations for estimation and kriging can all be carried out in terms of M × M matrices rather
than the typically much larger T × T matrices, where T is the number of observations.
Moreover the convolution approach can be extended readily to deal with nonstationary,2anisotropic, multivariate and spatio-temporal modeling (see e.g. Fuentes, 2002; Calderand Cressie, 2007; Majumdar et al, 2010; Sampson, 2010).
In the convolution approach the choice of kernels and knot locations plays a critical rolefor estimation and kriging performance. Higdon (2002) introduced a family of kernels thatgave standard classes of Gaussian, exponential and spherical covariances, and suggested
ad hoc choices for them in his ozone modeling examples. In most existing studies knot
locations are either regularly or randomly spaced points, the number being chosen so as
not to exceed computational limits. For given knot locations, the kernel parameters can be
estimated by either classical maximum likelihood or Bayesian methods, the specificationof knot locations and the estimation of kernel parameters being conducted separately.
In this paper we determine the first and second-order properties of the scalar random
fields defined above. The family of isotropic CARMA(p, q) kernels specified in Definition
3.1 generates a rich new class of spatial covariances which are not necessarily non-negativeand not necessarily monotonically decreasing. For model-fitting and kriging we use a
Bayesian Markov chain Monte-Carlo method to estimate the knot locations and CARMA
parameters simultaneously. We find in our examples that joint estimation of the kernelparameters and knot locations improves substantially on estimation and kriging in which
the knot locations are specified prior to kernel estimation. For the CAR(1) field it is foundthat the covariance function belongs to the Matérn class with smoothness parametern+ 1. For the general CARMA(p, q) field the spectral density is found explicitly and2
the covariance function is expressed as a modified Hankel transform (a one-dimensional
integral) of an explicitly defined function. It is evaluated explicitly for n = 1 and n = 3.
If the Lévy noise is derived from a compound Poisson sheet then it is a trivial matterto simulate a corresponding CARMA random field without any matrix operations. We
make use of this in our simulated examples. We observe also that if the CARMA kernel
is replaced by a Matérn kernel with smoothness parameter ν then the resulting field has
a Matérn covariance function with smoothness parameter ν + n2 .
The stationary and isotropic CARMA random fields introduced in this paper suggest further studies, in which we plan to extend the models to represent nonstationary
and anisotropic random fields. The paper is organized as follows. Section 2 introducessome background tools relating to Lévy sheets and the multivariate Fourier transform of
an isotropic function. In Section 3 we evaluate the first and second-order moments of a
CARMA random field driven by a second-order Lévy sheet. Section 4 discusses the joint
estimation of CARMA parameters and knot locations. Section 5 demonstrates the performance of CARMA random fields fitted to both simulated data and to a real-life example
involving Tokyo land price data.

3

2Preliminaries

Lévy noise and Lévy sheets have recenty received much attention in connection withstochastic partial diﬀerential equations (see e.g. Koshnevisan and Nualart (2008)). In thispaper we make use of them to define multiparameter versions of single-parameter CARMA(continuous-time autoregressive moving-average) processes.
•
If L is a scalar n-parameter Lévy sheet and L is the corresponding Lévy noise as definedby Koshnevisan and Nualart (2008)), Section 2, then we can follow their construction todefine an n-parameter scalar CARMA random field via the integral∫
gn (t − u)dL(u), t ∈ Rn ,(1)
Sn (t) =Rn

where the kernel gn : Rn → R is a suitably chosen analogue of the one-parameter CARMA∏
kernel. Kernels of the form gn (t) = ϕ(∥t∥) and gn (t) = ni=1 ψi (ti ) are of particular interest
but in this paper we shall restrict attention to the former and refer to the corresponding
random field as isotropic.
•
Important cases arise (i) when L is a Brownian sheet, in which case L is Gaussian
white noise and (ii) when L is a compound Poisson sheet, in which case L is the randommeasure on the Borel subsets of Rn given by,
L(A) =∞∑

Yi 1xi (A), A ∈ B(Rn ),

(2)

i=1

where xi denotes the location of the ith unit point mass of a Poisson random measure on
Rn and the sequence {Yi } is a sequence of i.i.d. random variables independent of {xi }.
A variety of fields are generated by (1) depending on the choice of L and of the kernel
gn . Our prime concern will be with second order fields and their first and second orderproperties which we shall apply to the modelling of multiparameter data sets. If L is asecond-order Lévy sheet then there exist constants µ ∈ R and σ > 0 such that for all
real-valued functions g and h in L1 (Rn ) ∩ L2 (Rn ),∫∫Eg(t) dL(t) = µg(t) dt
(3)
Rn

[∫

and
Cov

Rn

∫
g(t) dL(t),

Rn

]
∫
2
h(t) dL(t) = σg(t)h(t) dt.

(4)

Rn

Rn

If the kernel gn in definition (1) belongs to L1 (Rn )∩L2 (Rn ) then the first and second-orderproperties of the field Sn are given by∫ESn (t) = µgn (u)du
(5)
Rn

4

and
γn (t) := Cov(Sn (t), Sn (0)) = σ 2 gn ∗ gn (t), t ∈ Rn ,(6)i.e. σ 2 times the multi-dimensional convolution of the kernel gn with itself.
In order to compute the autocovariance function (6) we shall first compute the multivariate Fourier transform g̃n = Fn f of the kernel gn , defined, for gn ∈ L1 (Rn ), by∫−n/2
g̃n (ω) = Fn gn (ω) := (2π)
e−i⟨ω,t⟩ gn (t) dt, ω ∈ Rn .
(7)
Rn

Then by (6), provided (g̃n )2 ∈ L1 (Rn ), the autocovariance function γn is (2π)n/2 times theinverse Fourier transform of σ 2 (g̃n )2 , i.e.
∫2ei⟨ω,t⟩ g̃n (ω)2 dω, t ∈ Rn ,(8)
γn (t) = σRnso that the multivariate spectral density function fn of Sn isfn (ω) = σ 2 g̃n (ω)2 , ω ∈ Rn .

(9)

If the function gn is radial, i.e. if there is a univariate function ϕ such thatgn (t) = ϕn (∥t∥), t ∈ Rn ,(10)

then it is well-known (see e.g. Nowak and Stempak (2014)) that the Fourier transform of
gn is also radial and can be expressed asg̃n (ω) = ϕ̃n (∥ω∥) := H n2 −1 ϕn (∥ω∥), ω ∈ Rn ,where Hm denotes the modified Hankel transform,∫ ∞Jm (xy) 2m+1
Hm g(x) =g(y)ydy, x > 0, m ≥ −1/2,
(xy)m
0

(11)

(12)

and Jm is the Bessel function of the first kind of order m. Clearly the function g̃n2 is alsoradial so that (8) can be reexpressed asγn (t) = σ 2 (2π)n/2 H n2 −1 ϕ̃2n (∥t∥), t ∈ Rn .

(13)

If, for each n ∈ N, gn is a radial function in L1 (Rn ) with gn (t) = ϕ(∥t∥), t ∈ Rn , where
ϕ : R → R is a function independent of n, then we can write the Fourier transform (7) asg̃n (ω) = Fn gn (ω) = Fn ϕ(∥ω∥), n ∈ N,

(14)

where, by a very elegant result of Grafakos and Teschl (2013), the mapping Fn+2 is related
to the mapping Fn by
1 d
Fn ϕ(r), n ∈ N.
(15)
Fn+2 ϕ(r) = −r dr5Equation (15) diﬀers by a factor 2π from equation (2) of Grafakos and Teschl owing
to their slightly diﬀerent definition of the Fourier transform. Thanks to this equation it
becomes a simple matter to compute Fn ϕ, n ∈ N, from F1 ϕ and F2 ϕ.
We shall say that the field Sn defined by (1) is isotropic if the kernel gn satisfies (10) forsome function ϕ : R → R. (The function ϕ may depend on n but in all of our applications
it will denote a function independent of n.)3

Isotropic Lévy-driven CARMA random fields∏
Definition 3.1. Let a∗ (z) = z p + a1 z p−1 + · · · ap = pi=1 (z − λi ) be a polynomial of degree
p with real coeﬃcients and distinct zeroes λ1 , . . . , λp having strictly negative real parts and∏let b∗ (z) = b0 + b1 z + · · · bq z q = qi=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 )q∏and b(z) =(z 2 − ξi2 ),

i=1

i=1

the isotropic L-driven CARMA(p, q) field with autoregressive polynomial a and moving
average polynomial b is
∫ ∑pb(λr ) λr ∥t−u∥edL(u), t ∈ Rn ,(16)
Sn (t) =′Rn r=1 a (λr )where a′ denotes the derivative of the polynomial a and ∥t − u∥ denotes the Euclidean
norm of the vector t − u.
2∑
Remark 3.2. The kernel g(t) = pi=1 b(λi )eλi ∥t∥ /a′ (λi ) appearing in (16) is the specialcase, when a(z) has distinct zeroes, of the kernel,∑[]
Resz=λ ez∥t∥ b(z)/a(z) ,g(t) =
λ:R(λ)<0

where the sum is over the distinct zeroes of a(z) with negative real parts and Resz=λ f (z)denotes the residue at z = λ of the function f . When n = 1, it is the kernel of the noncausal CARMA(2p, 2q) process with autoregressive and moving average polynomials a(z)and b(z) respectively (see Brockwell and Lindner (2009)). For simplicity of exposition andwith no essential loss of generality we shall assume throughout this paper that the zeroesof a(z) are distinct. In the special case when a(z) = z 2 − λ2 , λ < 0, and b(z) = 1, theprocess S1 defined by (16) reduces to (2λ)−1 times the well-balanced Ornstein-Uhlenbeck
process of Schnurr and Woerner (2011). The definition (16) can therefore be regarded intheir terminology as a well-balanced CARMA process when n = 1 and as the generalizationof such a process to a field on Rn when n > 1. .
6For the compound Poisson L defined by (2),Sn (t) :=p∞ ∑∑b(λr )a′ (λr )i=1 r=1

Yi eλr ||t−xi || , t ∈ Rn ,(17)

an expression which permits very simple simulation of Sn on any hypercube or on certain
more complicated bounded subsets of Rn .
Remark 3.3. It is interesting to note that the kernelp∑b(λr )a′ (λr )r=1

eλr ∥t∥ = γ(∥t∥),where γ is the autocovariance function of the one-dimensional CARMA(p, q) process withautoregressive polynomial a∗ (z),moving average polynomial b∗ (z) and driving Lévy process L such that V ar(L1 ) = 1. This result can be derived from the representation (65) of
γ in Brockwell (2014).
Our first theorem establishes the first and second-order properties of the isotropic
CAR(1) field with autoregressive polynomial a(z) = z 2 − λ2 , λ < 0.
Theorem 3.4. If L is a second-order Lévy sheet satisfying (3) and (4) and if the defining
polynomials in (16) are a(z) = z 2 − λ2 (where λ < 0) and b(z) = 1, we shall refer to
Sn as an isotropic CAR(1) field with kernel gn (t) = (2λ)−1 eλ∥t∥ = ϕ(∥t∥), t ∈ Rn , where
ϕ(r) := eλr /(2λ), r ≥ 0. Thenn

E[Sn (t)] =µπ 2 Γ(n + 1),2λ|λ|n Γ( n2 + 1)

(18)

the spectral density of Sn isfn (ω) = σ 2 c2n (∥ω∥2 + λ2 )−n−1 , ω ∈ Rn ,
where
cn =


−2n/2−1 Γ( n+1 )/√π

if n is odd,

−2−n/2 Γ(n)/Γ( n )if n is even,

2

2

(19)

and the autocovariance function of Sn isγn (t) = σ 2 c2n( π ) n22n∥λt∥ 2 +1
K n +1 (∥λt∥), t ∈ Rn
|λ|n+2 Γ(n + 1) 2

where K n2 +1 denotes the modified Bessel function of the second kind of order n2 + 1.

7

(20)

∫
µ
Proof. From (5) we find at once that E[Sn (t)] = 2λeλ∥t∥ dt. Rewriting the integral inRnpolar coordinates then gives∫ ∞n
n−1 n/2µπ 2 Γ(n + 1)
µ
π
λr nr
E[Sn (t)] =
e
dr =
.
2λ 0
Γ( n2 + 1)2λ|λ|n Γ( n2 + 1)
In order to compute the spectral density and autocovariance function of Sn we shallfirst determine g̃n (ω) = Fn gn (ω) = Fn ϕ(∥ω∥) = ϕ̃n (∥ω∥), ω ∈ Rn , for n = 1 and 2, andthen use these, in conjunction with (15), to determine g̃n , n ∈ N. The spectral density of
Sn will then be found from (9) and the autocovariance function from (13).
The case n = 1:
Substituting ϕ(r) = eλr /(2λ) in (7) and evaluating the integral givesg̃1 (ω) = − √1, ω ∈ R.
2π(ω 2 + λ2 )(21)

The case n = 2: From (11) and (12) we have
∫ ∞
1
1
1
λ(·)
g̃2 (ω) =, ω ∈ R2 ,H0 e (∥ω∥) =eλy J0 (∥ω∥ y)y dy = −
2
2
3/2
2λ
2λ 0
2(∥ω∥ + λ )(22)
where the last equality can be found from tables of the zero-order Hankel transform (seee.g. Piessens (2010), table 9.1, equation(4)).
From (21) and (22) it follows that F1 ϕ(r) = − √12π (r2 + λ2 )−1 and F2 ϕ(r) = − 12 (r2 +λ2 )−3/2 . Successive application of (15) with n = 1, 2, 3, . . . , then gives,
Fn ϕ(r) = ϕ̃n (r) = cn (r2 + λ2 )− 2 , n ∈ N,
n+1

(23)

with cn defined as in (19). These equations, with (9) and (14) immediately give equation(19) for the spectral density of Sn .
Substituting from (23) into equation (13) and using Piessens (2010), Table 9.2, equation (12), we find immediately that the autocovariance function of Sn is given by (20) as
claimed.
Remark 3.5. The values of Kr (z), r ≥ 0, z ∈ C, are available in Matlab as besselk(r, z).
The functions Kn+ 1 , n ∈ N, take a particularly simple form. From Abramowitz and2
Stegun (2012), equations (10.2.17), we have√
π/(2z)e−z (1 + z −1 )(24)

π/(2z)e−z (1 + 3z −1 + 3z −2 ).

(25)

K 3 (x) =2andK 5 (x) =2√

8

Moreover, if we definehn (z) := (−1)n+1√
π/(2z)Kn+ 1 (z) , n ∈ N,2

we can readily calculate Kn+ 1 (z), for integers n > 2 from (24), (25) and the recursions
2(Abramowitz and Stegun (2012), equation (10.2.18)),
hn+1 (z) = hn−1 (z) − (2n + 1)z −1 hn (z).
Remark 3.6. Direct application of Theorem 3.4 gives the following results for the important cases n = 1, 2, and 3.
µE[S1 ] = − 2 ,λf1 (ω) =σ2,2π(ω 2 + λ2 )2

γ1 (t) =σ 2 −|λt|e(1 + |λt|).
4|λ|3

πµ
,
λ3

f2 (ω) =σ2,
4(∥ω∥2 + λ2 )3γ2 (t) =σ2π∥λt∥2 K2 (∥λt∥),
4
16λ

2σ 2
4πµ
,
f
(ω)
=
,
3
λ4
π(∥ω∥2 + λ2 )4

γ3 (t) =σ 2 π −∥λt∥e(3 + 3 ∥λt∥ + ∥λt∥2 ).
12|λ|5

E[S2 ] =E[S3 ] = −Remark 3.7. The general Matérn autocovariance function has the form (see e.g. Stein(1999))
γn (t) = σ 2 (∥at∥)ν Kν (∥at∥), a > 0, ν > 0.
If we define Sn as in (1) with a Matérn kernel, i.e. withgn (t) = (∥at∥)ν Kν (∥at∥), ν > 0,(26)

then we find, by arguments analogous to those used in the proof of Theorem 3.4, that thecorresponding mean, spectral density and autocovariance function of Sn are,()2( ) n2σ 2 (2a2 )2ν 2n−2 Γ n2 + ν
4πnν−1E[Sn ] = µ 2
Γ( + ν), fn (ω) =anda22(∥ω∥2 + a2 )n+2ν
σ2
γn (t) =
2

(

2π
a2

)n/2

(
)2
Γ n2 + νn(∥at∥) 2 +2ν K n2 +2ν (∥at∥).
Γ(n + 2ν)In other words the Matérn kernel with index ν generates a sheet with Matérn autocovariance function having index n2 + ν. The results of Theorem 3.4 can be derived from thesemore general results since11(2λ)−1 eλ∥t∥ = (2λ)−1 e−∥λt∥ = √ ∥λt∥ 2 K 1 (∥λt∥).
2
λ 2π

9

We next consider the isotropic CARMA(p, q) field defined by (16).
Theorem 3.8. If L is a second-order Lévy sheet satisfying (3) and (4) and if the isotropic∑r ) λr ∥t∥
CARMA field Sn is specified as in Definition 3.1, i.e. with kernel gn (t) = pr=1 ab(λ,′ (λ ) ernt ∈ R , thennp∑b(λi ) π 2 Γ(n + 1)
E[Sn (t)] = µ,(27)′ (λ ) |λ|n Γ( n + 1)
a
i
2
i=1
and the spectral density of Sn is


)2
( n−1

σ 2 2n−2 ψ1( 2 ) (r) ,πfn (ω) = ρn (∥ω∥) =
( n−2
)2

σ 2 2n−3 ψ2( 2 ) (r) ,if n is odd,
(28)
if n is even,

where r = − ∥ω∥2 , the superscripts on the functions ψ1 and ψ2 denote order of diﬀerentiation with respect to r, andψk (r) = −p∑

2λi b(λi )

2 k+1
′
2
i=1 a (λi )(−r + λi ), k = 1, 2.

(29)

The autocovariance function of Sn isγn (t) = (2π)n/2 H n2 −1 ρn (∥t∥),(30)

where Hm denotes the modified Hankel transform of order m as defined in (12) and ρn isdefined as in (28).
∑r ) λr ∥t∥
Proof. For the field Sn with kernel gn (t) = pr=1 ab(λ, t ∈ Rn , equation (5) gives′ (λ ) erESn (t) = µ∫p∑b(λi )a′ (λi )i=1

eλi ∥u∥ du.

Rn

The integral on the right was evaluated in the proof of Theorem 3.4. Referring to equation(18) we see at once that (27) holds.
n+1
The Fourier transform of eλi ∥t∥ , t ∈ Rn , is, from (23), 2λi cn (∥ω∥2 + λ2i ) 2 , ω ∈ Rn ,where the coeﬃcients cn are as in (20). The Fourier transform of gn is thereforeg̃n (ω) = cnp∑

2λi b(λi )

2
2 n+1
′
2
i=1 a (λi )(∥ω∥ + λi ).

With ψ1 and ψ2 defined as in (29) it is a straightforward calculation to show thatn−1(n−2)/2 2 √π ψ1( 2 ) (r), if n is odd,g̃n (ω) =



n−2

2(n−3)/2 ψ ( 2 ) (r), if n is even,210where r = − ∥ω∥2 . The spectral density of Sn is then fn (ω) = σ 2 g̃n (ω)2 , i.e. (28).
The autocovariance function of Sn is (2π)n/2 times the inverse Fourier transform of
the spectral density and, because the spectral density is radial, it can be expressed as the
one-dimensional integral (30).
A fast algorithm for the numerical computation of Hankel transforms was developedby Knockaert (2000). If the dimensionality n is odd it is possible to express the spectral
density of Sn in a simpler form and to use this to compute the autocovariance functions
explicitly in the particular cases when n = 1 and n = 3.
Corollary 3.9. The function ψ1 in (28) simplifies to
β(r)
,
α(r)

ψ1 (r) =(31)∏∏
where β(z) := qi=1 (z − ξi2 ) and α(z) := pi=1 (z − λ2i ). This implies, by (28) that whenn is odd the spectral density fn (ω) of Sn is a rational function of ∥ω∥2 , easily calculated
from (28).
Corollary 3.10. In the special cases n = 1 and n = 3 we have
[
]
p
2
∑
2
z|t| b(z)γ1 (t) = σ
Resz=λi e
a(z)2
i=1
and
γ3 (t) = −σ

2 2π

p
∑

∥t∥ i=1[
Resz=λi(a′ (z)b(z) − a(z)b′ (z))2ez∥t∥za(z)4].

Proof. From (28) and (31) with n = 1,(

σ2
f1 (ω) =
2π
Hence
σ2
γ1 (t) =2πβ(−ω 2 )α(−ω 2 )
(

∫ ∞
e

iωt

−∞

)2

β(−ω 2 )α(−ω 2 ).
)2dω.

Changing the variable of integration to z = iω, using the relations α(z 2 ) = a(z) andβ(z 2 ) = b(z) and using contour integration to evaluate the integral gives the required
expression for γ1 .
From (28) and (31) with n = 3,
f3 (ω) =2σ 2δ(− ∥ω∥2 ),π

where
δ(r) = α(r)−4 (α(r)β ′ (r) − α′ (r)β(r))2 .
11

Hence2σ 2γ3 (t) =π∫R3ei⟨ω,t⟩ δ(− ∥ω∥2 )dω.

Expressing this integral in three-dimensional polar coordinates and integrating out theangular components gives∫ ∞4σ 2
γ3 (t) =Imρδ(−ρ2 )2 eiρ∥t∥ dρ.
∥t∥−∞Rewriting δ in terms of the polynomials a and b and again using contour integration toevaluate the resulting integral gives the required expression for γ3 .
Example 3.11. A two-parameter CAR(2) random field.
If L is a two-dimensional second-order Lévy sheet satisfying (3) and (4) and S2 (t), t ∈ R2 ,is the two-parameter random field defined by (16) with a(z) = (z 2 − λ21 )(z 2 − λ22 ), λ1 =−1 + 2i, λ2 = −1 − 2i and b(z) = 1, then the kernel g2 is readily calculated as
(
)
e−∥t∥
1
g2 (t) = ϕ2 (∥t∥) =cos(2 ∥t∥) + sin(2 ∥t∥) .
20
2
From (28) and (29) we then find that the spectral density of S2 is
f2 (ω) = ρ2 (∥ω∥) =[]2σ2(∥ω∥2 + λ22 )−3/2 − (∥ω∥2 + λ22 )−3/22
2 22(λ1 − λ2 )and, from (30) and (12),∫ ∞γ2 (t) = 2πH0 ρ2 (∥t∥) = 2πyρ2 (y)J0 (ty) dy.
0

The latter integral is easily evaluated using the Matlab function besselj to compute J0 (ty).
Numerical evaluations of correlation functions of two-parameter CAR(1) and CAR(2)random fields will be demonstrated in Sec. 5.

4

Estimation and kriging

4.1

Parameter estimation

This section considers estimation of the parameters of the CARMA kernel when thedriving Lévy sheet is compound Poisson as in (17). The data consist of observationsof Sn (s) at the finite number of sampling points s1 , . . . , sT , scattered irregularly over abounded region D ⊂ Rn .

12

Our model for the data will be a modified form of the model (17) which allows inparticular for measurement errors. We shall write the CARMA kernel, with parameter
vector θ = (λ1 , . . . , λp )′ asgθ (s, x) =p∑b(λr )a′ (λr )r=1

eλr ||s−x|| .

Our model is then,Sn (sj ) = Zj β +M∑Yi gθ (sj , xi ) + εj , j = 1, . . . , T,

(32)

i=1

Yi ∼ iidN (0, τ 2 ), εj ∼ iidN (0, δ 2 ).
This model diﬀers from (17) in several respects. First, the mean is accounted for by theterm Zj β, corresponding to the regressor vector Z = (Z1′ , . . . , ZT′ )′ . Secondly, independent
normally distributed noise terms εj have been included to represent measurement errors,
known in the literature as nugget eﬀects. Finally, the ’knot points’, xi , are the points of
the underlying spatial Poisson process falling in the region D and their number, n(D),has the Poisson distribution. Conditionally on the value of n(D), the knot points are
independently and uniformly distributed on D. The number of knots and their locations
are crucial in the estimation procedure as discussed in section 4.2.
The model (32) resembles a standard regression model with the values of the kernel
function gθ (s, xi ) as independent variables. There are however diﬃculties which are notencountered in classical regression analysis. (i) The knots xi , and hence the independent
variables gθ (s, xi ), are stochastic. (ii) The independent variables depend on the parameterθ to be estimated. (iii) The regression coeﬃcients Yi are not fixed but random. Thesediﬃculties prohibit classical maximum likelihood estimation and so a diﬀerent approachis needed.
We start by evaluating the conditional likelihood function of the parameters (β, δ 2 , θ, τ 2 )given the number, M , and locations, {xi }, of the knots in the model (32). The conditionalvariance of Sn = (Sn (s1 ), . . . , Sn (sT ))′ is then()τ 2 Vθ Vθ′ + δ 2 IT = δ 2 IT + κ2 Vθ Vθ′ = δ 2 R,
say, where Vθ is the T by M matrix whose (p, q)th component is Vpq = gθ (sp , xq ) andκ2 = τ 2 /δ 2 . The conditional likelihood is()
Sn ∼ N Zβ, δ 2 R .
Calculation of the inverse and determinant of R requires an infeasible O(T 3 ) operations.
We apply Corollaries 18.1.2 and 18.2.10 of Harville (1997) to calculate them in O(M 3 )13

operations as{}(δ 2 R)−1 = δ −2 IT − Vθ (κ−2 IM + Vθ′ Vθ )−1 Vθ′ ,andδ 2 R = δ 2T κ2 Vθ′ Vθ + IM .
These expressions yield, in O(M 3 ) operations, the conditional log likelihood of Sn (afunction (of Ω = (β, δ 2 , θ, κ2 , {xi })) in O(M 3 ) operations as
log L(Sn |Ω) = −T
1
1log δ 2 − log |R| − 2 (Sn − Zβ)′ R−1 (Sn − Zβ),
2
2
2δ

(33)

It is infeasible to estimate (β, δ 2 , θ, κ2 ) by maximizing (33) unless x1 , . . . , xM , are fixedand known, which leads us to take a Bayesian approach.
For given M , the knot locations, {xi }, are independently and uniformly distributed on
D. Taking this as the prior distribution for {xi } opens the way to Bayesian estimation of
the model. Bayesian inference requires us to specify priors for the other model parameters.
The Gibbs sampler is used to draw samples of parameters from the following posteriorp(Ω|Sn ) ∝ L(Sn |Ω)p(β)p(δ 2 )p(θ)p(κ2 )p({xi }).
We assign the conjugate priors N (β0 , Σ0 ) and Ga(a1 , a2 ) for β and δ −2 , respectively. We
update β from N (µβ , Σβ ), where−1−2 ′ −1Σ−1
β = Σ0 + δ Z R Z,

and()µβ = Σβ δ −2 Z ′ R−1 Sn + Σ−10 β0 ,
and update δ −2 from the gamma distribution with parameters a1 + T /2 and a2 + 1/2(Sn −Zβ)′ R−1 (Sn − Zβ). For the other parameters without conjugate prior distributions, wedraw samples using the Metropolis-Hastings (MH) algorithm. We employ the vague priordistribution, uniform on the parameter set for (θ, κ2 ), and use random walk chains driven
by normal variables with mean 0 to update them in the MH algorithm. Priors for the knot
locations and the updating procedure will be discussed in detail in the following section.

4.2

Knot specification

Specification of the knot locations {xi } plays a critical role in estimation and kriging using
CARMA models. Gelfand et al. (2012) proposed a way in the predictive process approach14

to specify them separately from parameter estimation. Zhang et al. (2015) employeda procedure for joint estimation for knots and parameter in the context of composite
approach. We shall follow the latter approach with a modification to allow for our assumed
uniform distribution of knot locations.
Suppose we have a set U = {u1 , . . . , uN } ⊂ D of knot candidates, which are independent uniformly distributed points in D. We take the prior distribution of the knots to bethat of points selected by independent Bernoulli sampling from U with success probability p. This distribution approximates that of an independent uniformly distributed set ofpoints K0 in D with the expected number of points, EM0 = N p. It approximates the
distribution of points in D generated by an n-dimensional Poisson process, conditionalon M0 points falling in D. The value of p is chosen to make operations with matrices ofdimension EM0 computationally feasible.
Given an initial knot set K0 chosen from U as described, we use a Metropolis-Hastings
algorithm to obtain successive knot sets K1 , K2 , . . . . For any subset K ⊂ U we shall use
the notation K c for the set U \ K and n(K) for the cardinality of K. The recursive
construction of the knot sets can then be described as follows for t = 1, 2, . . . .
c(1) Draw A and B from Kt−1 and Kt−1
by Bernoulli trials with success probabilities q1c
and q2 (t) = q1 n(Kt−1 )/n(Kt−1), respectively, and define K ′ by (Kt−1 \ A) ∪ B. (Ifn(Kt−1 ) = 0 or q2 (t) ≥ 1 then K ′ is obtained by independent Bernoulli trials withsuccess probability p from U ).

(2) Compute the acceptance ratio()
L(Sn |K ′ )p(K ′ )J(K ′ → Kt−1 )α = min 1,,
L(Sn |Kt−1 )p(Kt−1 )J(Kt−1 → K ′ )where L(Sn |K) is the conditional likelihood in (33) given the knot set K, and p(K)is the prior distribution assigned for the knot set K, i.e.,cp(K) = pn(K) (1 − p)n(K ) ,and J(K → K ′ ) is the probability of a transition from K to K ′ , i.e.,′n(K ′ \K)J(K → K ′ ) = (1 − q1 )n(K∩K ) q1′c′c(1 − q2 (t))n(K ∩K ) q2 (t)n(K \K ) ,ccslightly modified if q2 (t) ≥ 1. The other quatities p(K ′ ), J(K ′ → K) are defined in
the same way.
(3) Set Kt equal to K ′ with probability α and to Kt−1 with probability 1 − α.

15We note that p is chosen so that the expected number of initial knot points is within
computational feasibility (i.e. so that N p is of the order of one hundred). The probability
q1 is chosen to be small (typically between .01 and .03) so as to give reasonable acceptanceratios and the probability q2 (t) is designed to make En(Kt ) equal to En(K0 ) = N p for
all t.

4.3

Kriging

Kriging, i.e. prediction of the field at points for which no observations are available, is conducted in the course of the Gibbs sampling jointly with the estimation of the parameters.
Here we demonstrate kriging at a point s0 based on the model in (32). For a given Ω at a
step in the Gibbs sampling, we draw Y = (Y1 , . . . , YM )′ from the posteriors N (µY , δ 2 ΣY ),
where
−2
′
Σ−1
Y = κ IM + Vθ Vθ ,andµY = ΣY Vθ′ (Sn − Zβ).
Using the sampled Y , we construct the kriged value at s0 byŜn (s0 ) = Z0 β + Vθ (s0 )Y,
for the regressor Z0 at s0 . Conducting the kriging in the course of the Gibbs sampling,
we obtain posterior samples of the kriged value at s0 . Notice that all the procedures forthe kriging as well as in the Gibbs sampling are conducted using operations with M × M
matrices.

5

Empirical studies

5.1

Data generation

Simulation of the restriction of a compound Poisson CARMA random field to any boundedmeasurable subset D ⊂ Rn can be carried out by truncation of the infinite series in (17)as described in the following paragraph. Provided the compound Poisson sheet has finite
second moments, the covariance function of the CARMA field can be computed from
Theorem 3.8.
If the intensity parameter of the compound Poisson sheet is c then the number n(D′ )of knots contained in any measurable set D′ with volume |D′ | is simulated as a Poisson
16

random variable with mean c|D′ | and the knot locations are then simulated as n(D′ )independent and uniformly distributed points in D′ . If D′ is taken to be a suﬃciently
large set containing D, e.g. a large hypercube containing D, then the influence of theknots outside D′ will have little influence on the CARMA field restricted to D and canbe neglected in the sum (17). (The required size of D′ relative to that of D depends
on the zero of a(z) with smallest absolute value). The CARMA field restricted to D isthen obtained using the sum (17) truncated to include only the simulated knots in D′ . It
remains only to simulate the n(D′ ) i.i.d. random variables Yi to complete the evaluationof the simulated field at any point in D. If we take Yi to have the normal distribution
with mean zero and variance τ 2 then, in the notation of (3) and (4), µ = 0 and σ 2 = cτ 2 .
Using Remark 2.6 it is also a simple matter to simulate spatial data with the Matérn
covariance function,
γn (t) = ψ 2 (||at||)n/2+2ν Kn/2+2ν (||at||), ν > 0,
using the procedure in the last paragraph with the CARMA kernel in (17) replaced by
the Matérn kernel in (26), and with c and τ 2 chosen so that( 2 )n/2aΓ(n + 2ν)cτ = 2ψ.
2πΓ(n/2 + ν)2
2

5.2

2

Simulation studies

This section examines the empirical performance of our proposed estimation and kriging
procedure using simulated compound Poisson CARMA fields on R2 in the model (32). In
particular, we shall compare its performance with the corresponding procedure when theknots are fixed, rather than sampled.
We generated 100 sets of 1100 irregularly spaced points on the disc D = {(x, y) ∈2
R |(x − 50)2 + (y − 50)2 < 402 }. At each set of 1100 points, values of a CARMA random
field were generated to represent observations of the field at those points. The 1100 data
points so obtained were divided randomly into two sets of sizes 1000 and 100. The first 1000
points were used for estimation, while the remaining 100 were used for evaluating kriging
mean squared errors. The values of the field at the 1100 locations were generated from(32) with zero mean function. The values of the field were simulated using the proceduredescribed in Section 5.1 with D′ = [0, 100]2 ⊃ D. The intensity of the compound Poissonprocess was taken to be .02 so that the number n(D′ ) of knots uniformly and independently
distributed on D′ has the Poisson distribution with mean 200.
The following three CARMA kernel functions were employed in the the simulations.
All the kernels were normalized to be 1 at the origin to guarantee the identifiability of τ 2 ,the variance of Yj .
17

0.6

1.0

kernels

−0.2

0.2

Case1
Case2
Case3

0

5

10

15

20

1.0

correlations

−0.2

0.2

0.6

Case1
Case2
Case3

0

5

10

15

20

Figure 1: The kernel and resulted autocorrelation functions of the three CAR models usedin the simulation studies.
Model 1. CAR(1) kernel with defining polynomial a∗ (z) = z + α1 :g(s, x) = e−α1 ||s−x|| , α1 > 0.
Model 2. CAR(2) kernel with defining polynomial a∗ (z) = (z + α1 )(z + α2 ):()
g(s, x) = (α1 − α2 )−1 α1 e−α2 ||s−x|| − α2 e−α1 ||s−x|| , α1 > α2 > 0.
Model 3. CAR(2) kernel with defining polynomial a∗ (z) = (z + α1 + iα2 )(z + α1 − iα2 ):{}α1−α1 ||s−x||g(s, x) = ecos(α2 ||s − x||) +sin(α2 ||s − x||) , α1 > 0, α2 > 0.
α2
The true values for the parameters were α1 = 0.3 in Model 1, α1 = 0.5, α2 = 0.3 inModel 2 and α1 = 0.2, α2 = 0.4 in Model 3. We used the values δ 2 = 1 and κ2 = τ 2 /δ 2 =16 to simulate samples from (32). The kernel functions and corresponding correlation
functions in Models 1-3 are shown in Figure 1. These show that the CAR(2) kernels aresmoother at the origin than the CAR(1) kernel at the origin, and that the CAR(2) kernel
with complex roots can yield negative values for both the kernel and correlation functions.
We fitted the model (32) with unknown constant mean function µ to the first 1000
simulated values and evaluated the mean squared error of kriging for the remaining 100
simulated values. The prior distribution for our knot locations in D corresponded to
Bernoulli trials with success probability p = 0.1 applied to the first 1000 data points, sothat the expected initial number of knots was 100. We chose q1 = 0.01 for the updating ofknots in the MH algorithm with the initial knot set drawn as already described. We ran
1000 iterations of Gibbs sampling after a burn-in period of 1500 iterations. 100 posterior
samples for the parameters and kriging were collected by thinning, using every 10th18

selected knots
α1
α2
κ2
δ2
Model 1
0.3
16
1
0.23
10.41 1.06
medianquart. dev. 0.01
1.34 0.05
Model 2
0.5 0.3
16
1
median0.44 0.28 17.87 1.06
quart. dev. 0.04 0.03 2.67 0.04
Model 30.2 0.4
16
1
median0.19 0.39 14.33 1.11
quart. dev. 0.01 0.01 1.79 0.04

benchmarks by fixed knots
α1
α2
κ2
δ2
0.3
16
1
0.08
15.47
1.91
0.02
7.75
0.15
0.5 0.3
16
1
0.18 0.14 88.58
1.90
0.02 0.01 34.80
0.14
0.2 0.4
16
1
0.13 0.37 4.86
3.24
0.01 0.01 1.63
0.38

Table 1: The median and quartile deviations of the posterior medians of our proposed
estimators conducted for 100 sets of 100 samples in Models 1-3. The corresponding values
with the knots fixed at a randomly chosen 100 locations are shown for comparison.
iteration. The median and quartile deviations of the posterior medians for 100 sets of1000 simulated samples are shown in Table 1, and the mean squared errors of kriging
over 100 sets of 100 samples are listed in Table 2. Estimation and kriging with the knots
fixed at a randomly chosen 100 locations were conducted as benchmarks, and the mediansand quartiles are listed in Tables 1-2 for comparison with our algorithm. In Table 2, asa further benchmark, we included those of simple weighted averages, whose specific form
at u is
( T
)−1 T
∑
∑
wj
wj S(sj ),
(34)
j=1

j=1

for the weight function
wj = e−||u−sj || /(bdw) ,2

2

where bdw denotes the bandwidth controlling the rate of decay of the weights.
We see at once that our knot specification procedure improves dramatically on the useof fixed knots for both estimation and kriging. The biases of the estimators are greatly
reduced although in Model 1 there is a noticeable negative bias in the estimator of α1 . Thismay be due to the lack of smoothness of the CAR(1) kernel at zero. It can be reduced
by increasing the value of p to increase the number of knots. In general it seems that
kernels which are less smooth at the origin require more knots for accurate estimation.
The improvement in kriging performance achieved by the CAR models over the weighted
averages was more noticeable in Models 2 and 3 than in Model 1, which suggests that the19

CAR
selected fixed
Model 1
medianquart. dev.
Model 2
medianquart. dev.
Model 3
medianquart. dev.

weighted average
bdw=1
2
3

4

1.38
0.15

2.05
0.29

1.85
0.19

1.45
0.16

1.50
0.18

1.69
0.19

1.33
0.14

2.08
0.25

1.94
0.21

1.51
0.16

1.60
0.18

1.94
0.20

1.50
0.18

3.37
0.53

2.45
0.33

2.07
0.26

2.78
0.33

3.97
0.46

Table 2: The medians and quartile deviations of the mean squared errors of kriging. The
medians were the medians of posterior samples and are based on 100 sets of 100 samples
using each of the models 1, 2 and 3.. For comparison we also show the corresponding
results for weighted averages as defined in (34) with bandwdths 1, 2, 3 and 4, and for
CAR kriging with fixed knots..
weighted average procedure with suitably chosen bandwidth does a better job of capturing
the covariance structure of CAR(1) data than of the data generated by higher-order CAR
models.

5.3

Land price data analysis

This section demonstrates the application of the families of isotropic CAR models specified
in Models 1, 2 and 3 to the analysis of Tokyo land-price data. Public land prices in Tokyo,
sampled at points irregularly scattered over Tokyo’s 23 wards, are used. Since 1970 thegovernment of Japan has annually published land prices per square metre as of January
1st at hundreds of thousands of sampling points scattered irregularly all over Japan. They
are evaluated using a combination of transaction records, incomes and cost accounting
methods and do not therefore coincide with real market prices. They are published in
March every year on the government web page in order to help in the planning of public
works, inheritance tax and related economic issues.
This paper focuses on land prices in Tokyo’s 23 wards collected in 2015 fromhttp://nlftp.mlit.go.jp/ksj/old/old datalist.html (in Japanese).
The set of 1247 sampling points is shown in Figure 2. The price Pt (s) at location s in the20

35.80
35.75
35.70
35.55

35.60

35.65

latitude

●●
●
●
●
● ● ●
●
●● ●
●
● ●●
●
●●● ●
●●
●●
●● ●
● ●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ● ●●
●
● ●●
●●
●● ●
●
●
●
●
●
●
●
●●
●●
●
●
● ●
● ● ● ● ● ●● ●● ●
●●● ●
● ●
●
● ●●
●
●
●
● ● ● ●● ●●
●
● ● ●
●
●
●
●
●● ●
●● ●
●● ● ●●● ● ● ●●●●●● ● ●●
●
●
● ●
●
●
●
●
●
●● ● ●● ●● ●● ● ●● ●●● ●● ●●● ● ●
●
●
● ● ●
●
●
●
●
●
●
●
●
●
●● ●
● ●● ● ●●●●●●●● ● ●
● ●● ● ●●●●●●● ●●●●● ● ●● ●●● ● ●●●● ●● ●●
● ●●● ● ●●●● ●● ● ●●●●●● ● ●●● ●●● ● ●●●● ●● ● ● ●●● ●● ● ● ● ● ●● ●●● ● ●●●● ●●●● ● ● ● ● ●●● ●●●● ●●● ● ● ●● ●● ●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●● ●●●● ●
● ●
●
● ● ●
●●
●
● ●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
● ● ●●
● ●●
● ●
●●
●●
●
● ● ●
● ● ●● ●● ● ●●●●●●● ●● ● ● ● ●
● ●●
●● ●
●
● ● ●
●● ● ● ● ● ● ●●●●●●●● ● ● ●
●● ●
●
● ●●
●
● ● ● ●● ● ●
●
●●
●●●
● ●
●
●
●●●● ● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●● ●●●● ●●●●● ●● ● ●●●●● ●●● ●●●●
● ●● ● ● ●
●
●
●
●
●
● ●
●
●● ● ●● ● ●
● ● ●
●●●
●● ● ●● ●●●●●● ● ● ●●
●●● ●
●
●
●● ●
●●
●
●
● ●
●
●
●
●●
●
●● ●● ●● ● ●●
●● ●
●
●
●
● ●● ● ●●●● ●●
●
● ●
● ●●●
● ● ●●●● ●● ● ●●● ● ●●●●●● ●●● ●●●●● ●● ●●● ● ●●● ●●●●●●● ●●●●
● ●
●● ●
●●
●● ●
●
●●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●● ●
●
●
●● ●
●
●
●●●●●●●●●● ●● ●● ●● ●● ● ● ●●
●
● ● ●
●
● ●
●
●●●●●
●●
●
●●
●●● ●● ● ● ●●●●●●●●● ● ● ● ● ● ●
● ●●
●●●
●●
● ●
●
●●
●
●●
● ● ●●
●●
●●
● ●
●
●
●
●
●●
●
●
●
●● ● ● ● ●● ● ●●●●●● ● ●●
●
● ●
●●●●
●●
●●●
●
● ●
●
● ●
●
●●
●
●
● ●●
●
● ● ● ●● ● ●
●
●●● ●
●
●●●●● ● ● ● ● ●
●●
●
●
● ●●
●
●● ●●
● ● ●●●●●●●●● ●● ● ●● ●● ●● ●● ●● ●●● ● ● ●●●●
●
●
●
●●●
●
●
●● ●●
●
● ● ● ●●●● ●●●●● ● ●●
●● ●
●
● ●
●● ●
●
● ●
●●
● ●●● ●●
●●
●
●
●
●
●
●
●
●
● ●●
●
●
●● ● ●● ●● ● ● ●●●●● ● ●● ●●● ●●●
●●● ●
●●
●● ●
●●
● ●
●
● ●
●
●
●
●● ●● ● ●
●
●
● ●●
● ●
●
● ●●
●●
● ●
●
●● ● ●●●●●● ●●●● ● ● ● ● ●●●●●● ● ● ●●
● ●●
●
●●
●●
● ●● ●●●
●
● ●●
●
●
● ●●●
● ● ●
● ●● ● ● ●●●● ●● ● ●●
●●
●
● ●
●
●
● ●●●
●●●
● ●
● ●
●●
● ●● ● ●● ●●● ●● ● ●●● ●●●●● ● ●
●● ●
●
●
● ●
●● ●
●
●●●●
●
●
●● ● ●●●●●●● ●● ●
●
●
●
●
●
●
●
● ● ●
●
● ●
●●●
●
●
●● ● ●●
● ●●●
●● ●
●
●● ●
●
● ●
●

139.60139.65139.70139.75139.80139.85139.90longitude

Figure 2: 1247 sampling points for public land price data in Tokyo 23 wards collected in2015.
year t = 2015 was transformed to a log return rt (s), given byrt (s) = log Pt (s) − log Pt−1 (s), t = 2015.
The transformed great circle distance (Gneiting, 1999) was used to measure the distanced(s, s′ ) between the locations s and s′ . Thusd(s, s′ ) = 2r sin(∆ϕ/2),
where r is radius of the earth in kilometers and ∆ϕ ∈ [0, π] is the central angle between s
and s′ . We divided the 1247 locations randomly into two groups. The first group of 1147
was used for model fitting and the second group of 100 was used for evaluating kriging
mean squared errors.
We fitted the model (32) with constant mean function µ and CAR models 1, 2 and 3
to the returns at the first 1147 locations and evaluated the mean squared errors of kriging
at the remaining 100. We chose p = 100/1147 to give an expected number of knots equal
to 100 and q = 0.01 in the knots updating step. We ran 10000 iterations of the Gibbs
sampling after a burn-in period of 10000 iterations. 1000 posterior samples were collected
and these were thinned by using every 10th iteration.
21

α1

selected knots
α2
κ2
δ2

M

logL

Model 1
median1.03- 22.14 0.49 113 209.7
quart. dev. 0.04
- 3.64 0.01 4.5
7.7
Model 2
median15.10 1.02 18.41 0.54 101 168.3
quart. dev. 1.41 0.05 3.89 0.02 5.5 12.5
Model 3
median2.11 0.10 9.86 0.63 105 97.8
7.1
quart. dev. 0.06 0.07 1.24 0.01 3.5

benchmarks by fixed-knots
α1
α2
κ2
δ2
logL
0.22
0.02

-

2.25
0.54

1.07
0.01

-117.2
0.6

1.51
0.23

0.29
0.04

2.05
0.45

1.04
0.01

-112.2
0.7

0.55
0.04

0.08
0.05

1.83
0.38

1.04
0.01

-113.4
0.8

Table 3: CAR fitting by Models 1-3 for 1147 samples of log returns of land prices inTokyo. The median and the quartile deviation for the posterior samples were evaluatedin comparison with those of the cases of fixed knots with 100 randomly chosen locations.
In Table 3, we show the medians and quartile deviations of posterior samples with
their log likelihoods as given in (33). We conducted the fitting also with knots fixed at
100 randomly chosen locations and listed the results in Table 3 as benchmarks. In Table
4, we show the mean squared errors of kriging using the posterior medians of the kriging
samples as kriged values. The performance of these kriged values is compared with thatof the weighted average (34) and the corresponding kriged values obtained using fixedknots.
We find that in terms of log likelihood our knot selection procedure improves significantly on the fixed knot procedure and the CAR(1) model outperforms the CAR(2)models. The CAR(1) model with our knot selection procedure has the best kriging meansquared error of all considered. It is interesting that the improvement in kriging performance of the knot selection procedure over the fixed knot procedure is not nearly as
pronounced as in the simulation study.
The kernel and autocorrelation functions of the best fitting CAR(1) are shown in
Figure 3 and the corresponding smoothed values, obtained by kriging over the entire
sampling region of Tokyo’s 23 wards are shown in Figure 4. According to Figure 3 thelog returns of land prices have autocorrelations smaller than .05 for lags greater than
6 kilometres. Figure 4 clearly demonstrates the interesting shape of the returns surface,with the highest returns in the centre of Tokyo, gradually decreasing as we move towards
the remote suburbs.

22

MSE

selected knots
fixed knots
weighted average
Model1 Model2 Model3 Model1 Model2 Model3 bdw=1
2
3
4
0.93
0.96
1.10
1.02
1.03
1.05
1.05 1.03 1.08 1.17

0.8

1.0

Table 4: Kriging mean squared errors for randomly chosen 100 samples in log returns ofland price data in Tokyo, where kriging was conducted as the median of posterior samples.
The weighted averages with bandwidth of 1, 2, 3 and 4 and kriging with fixed knots wereconducted as benchmarks.

0.0

0.2

0.4

0.6

CAR(1) kernel
CAR(1) acf

0

2

4

6

8

km

Figure 3: The identified kernel and autocorrelation functions by Model 1 for log returnsof land prices in Tokyo 23 wards.

6

Discussion

This paper proposes the use of isotropic CARMA random fields driven by Lévy sheets
and specified by a kernel function defined in terms of autoregressive and moving average
polynomials of orders p and q respectively where 0 ≤ q < p. When the Lévy sheet has
finite second moments this class generates a rich source of isotropic covariance functions
on Rn with spectral densities which can be written explicitly. The covariance functions are
not necessarily non-negative or monotone. The covariance function of a CAR(1) field has
a simple explicit expression which belongs to the Matérn class with smoothness parameter
2. The covariances of more general CARMA fields can be expressed as modified Hankel
transforms which can be explicitly evaluated when n = 1 or n = 3.
Compound Poisson CARMA random fields are particularly useful as their restriction
to any bounded measurable subset of Rn can be approximated by a finite sum of terms
involving no matrix operations. Moreover estimation and kriging can be jointly conductedby Gibbs sampling requiring O(M 3 ) operations where M is the number of terms used in
the approximation to the random field. The Metropolis-Hastings algorithm can be used

23

8
35.80

35.75

6

35.70
4
35.65

2
35.60

35.55

0
139.60139.65139.70139.75139.80139.85139.90Figure 4: The smoothed figure of log returns of land price in Tokyo 23 wards by Model 1.
to specify the knot locations in the Gibbs sampling and the procedure works well in thesimulated examples and in terms of kriging mean squared error in the Tokyo land-priceexample.

References[1] Abramowitz, M. and Stegun, I. A. (2012). Handbook of Mathematical Functions: with
Formulas, Graphs, and Mathematical Tables. Dover Publications, New York.
[2] Brockwell, P.J. and Lindner, A. (2009). Existence and uniqueness of stationary Lévydriven CARMA processes, Stoch.Proc.Appl. 119, 2660–2681.
[3] Brockwell, P.J. (2014). Recent results in the theory and applications of CARMA
processes, Ann. Inst. Statistical Mathematics 66, 647–685.
[4] Calder, C.A. and Cressie, N. (2007). Some topics in convolution-based spatial modeling, In: Proceedings of the 56th Session of the International Statistical Institute,
Lisbon, International Statistical Institute. Voorburg, the Netherlands.

24[5] Fuentes, M. (2002). Spectral methods for nonstationary spatial processes. Biometrika,
89, 197-210.
[6] Gelfand, A. E., Banerjee, S. and Finley, A. O. (2012). Spatial Design for Knot Selection in Knot-Based Dimension Reduction Models, in Spatio-Temporal Design: Advances in Eﬃcient Data Acquisition (eds J. Mateu and W. G. Müller), John Wiley& Sons, Ltd, Chichester, UK.
[7] Gneiting, T. (1999). Correlation functions for atmospheric data analysis, Quarterly
J. Roy. Meteorological Soc. 125, 2449–2464.
[8] Grafakos, L. and Teschl, G. (2013). On Fourier transforms of radial functions anddistributions, J. Fourier Anal. Appl. 19, 167–179.
[9] Harville, D. A. (1997). Matrix Algebra from a Statistician’s Perspective. Springer,
New York.
[10] Higdon, D. (1998). A process-convolution approach to modeling temperatures in the
North Atlantic Ocean. Environmental and Ecological Statistics, 5, 173-190.
[11] Higdon, D. (2002). Space and space-time modeling using process convolutions. Quantitative Methods for Current Environmental Issues., 37-56.
[12] Khoshnevisan, D. and Nualart, E. (2008). Level sets of the stochastic wave equation
driven by a symmetric Lévy noise,Bernoulli 14, 899–925.
[13] Knockaert, L. (2000). Fast Hankel transform by fast sine and cosine transforms: theMellin connection,IEEE Transactions on Signal processing 48, 1695–1701.
[14] Majumdar, A., Paul, D. and Bautista, D. (2010). A generalized convolution model
for multivariate nonstationary spatial processes. Statistica Sinica, 20, 675-695.
[15] Matérn, B. (1960). Spatial Variation. Meddelanden frán Statens Skogsforskningsinstitut, 49, Almaenna Foerlaget, Stockholm. Springer-Verlag, Berlin.
[16] Nowak, A. and Stempak, K. (2014). A note on recent papers by Grafakos and Teschl,
and Estrada, J. Fourier Anal. Appl. 20, 1141–1144.
[17] Paciorek, C.J. and Schervish, M. J. (2006). Spatial modeling using a new class of
nonstationary covariance functions. Environmetrics, 17, 483-506.
[18] Piessens, R. (2010). The Hankel transform. Chapter 9 of The Transforms and Applications Handbook: Third Edition, ed. Alexander D. Poularikas, CRC Press, Boca
Raton.
25

[19] Sampson, P.D. (2010). Cnstructions for nonstationary spatial processes. Handbook of
Spatial Statistics. Edited by Gelfand, A.E. Diﬄe, P., Fuentes, M. and Guttorp, P.,
119-130. CRC Press, Boca Raton.
[20] Schnurr, A. and Woerner, J. H. C. (2011). Well-balanced Lévy-driven OrnsteinUhlenbeck processes, Statistics and Risk Modelling 4, 343–359.
[21] Stein, M. L. (1999). Interpolation of Spatial data: Some Theory for Kriging. Springer,
New York.
[22] Zhang, B., Sang, H. and Huang, J.Z. (2015). Full-scale approximations of spatiotemporal covariance models for large datasets. Statistica. Sinica 25, 99–114.

26