Discussion Paper No. 30

WAVELET ANALYSIS OF IRREGULARLYSPACED DATA AND ITS SPATIO-TEMPORALEXTENSION
YASUMASA MATSUDA

November, 2014

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WAVELET ANALYSIS OF IRREGULARLY SPACED DATA
AND ITS SPATIO-TEMPORAL EXTENSIONYASUMASA MATSUDA

Abstract. This paper proposes an orthogonal basis expansion by Haar wavelets
for irregularly spaced spatial data and conducts its spatio-temporal extension
when spatial points are not necessarily the same at each temporal point. One
crucial contribution of this paper is discovery of algorithm to construct orthogonal basis for any irregularly spaced data via 2 dimensional Haar wavelets. We
construct a spatial model by imposing a prior distribution on coeﬃcients of
the orthogonal expansion. In conducting spatio-temporal extension, we fit asimple autoregressive model not to data itself but to the coeﬃcients of the orthogonal expansion, which makes it possible to handle spatial sampling points
that are not necessarily the same at each temporal point. The spatio-temporalmodel can express nonstationary spatial behaviors by allowing parameters dependent spatially. Advantage of the models is superb parameter estimationeﬃciency that is attained by independency of wavelet transform providing efficient Gaussian likelihood evaluation. Finally we demonstrate the wavelet
modeling for land price data in Tokyo from 2003 till 2014, yearly series of public land prices per m2 over about 8000 sampling points scattered irregularly
all over Tokyo areas.

1. Introduction
Spatial or spatio-temporal data have recently been attracting growing attentions
in areas of both natural and social sciences by rapid progress of data collection
technologies. This paper considers a modeling for spatial data observed at sampling
points irregularly scattered over so broad region that stationary covariances cannot
be assumed, and aims to conduct its spatio-temporal extension.
Classical geostatistical methods provides several kinds of spatial covariance models such as Spherical, Exponential, Gaussian, and Matérn classes to analyze spatial
data (Zimmerman and Stein, 2010). Estimation of parameters in the covariancemodels by Gaussian likelihood requires evaluation of determinant and inverse of
the covariance matrix whose dimension is as many as the sample size. It follows
that exact Gaussian likelihood is in practice diﬃcult to evaluate especially for large
data set.
In order to avoid the computational diﬃculties, several methods for approximations have been proposed in the literatures. Typical ones that represent them are
the covariance tapering by Kaufman et. al. (2008), the predictive process approachby Banerjee et. al. (2008) and the composite function approach by Bai et. al.
(2012). The three methods commonly employ likelihood functions approximated inclever ways to avoid the diﬃculties caused by large dimensional covariance matrices.
Although they work reasonably well for large spatial data set, trivial extensions to
Key words and phrases. forecasting. Haar wavelet. Kalman filter. kriging. state space model.
Whittle likelihood function.
1

2

YASUMASA MATSUDA

spatio-temporal data are not eﬃcient enough to approximate the likelihoods foras many as several ten thousand samples, which is normally supposed in spatiotemporal extensions.
Reduced-rank random eﬀects approach, a method of dimension reduction to let
Gaussian likelihood evaluation more eﬃcient, has been used for analysis of large
spatial data set in the literatures as an alternative (Wikle, 2010). Let Z be a ndimensional observation vector possibly for large n. For a n × p regressor matrix
H with p < n, the approach employs the following multiple regression models
Z = Hβ + ε,
where β is assumed to follow normal distribution with mean 0 and variance Σ, whichis possible to regard as a prior distribution for β. There are variety of choices for
H and Σ, which specify practical performances of modeling for large spatial dataset.
Let us introduce principal three approaches for the choice of H and Σ in spatial
data modeling. The first one is Karhunen-Loéve expansion approach. Assumingthat ε = 0, we haveV ar(Z) = HΣH ′ ,which means that H and Σ are obtained by a singular value decomposition for
V ar(Z) as orthogonal and diagonal matrices, respectively. Though the approachgives a natural interpretation for the selection of H and Σ in the same way as principal component analysis does, we have to tackle the following non-trivial problems
in practical applications. V ar(Z) is in practice unknown and should be estimated.
The estimation is not trivial except for the case such as independent realizationsof spatial data is available. The matrix H constructed by the eigenfunctions aredefined just over sampled points and requires some technique to apply it for kriging.
For papers handling those in Karhunen-Loéve expansion approach, see, e.g., Chen
and Huang (2013) and references therein.
The second one is non-orthogonal basis function approach that specifies the meanfunction H(s)β at a spatial point s withp∑h(s, rj )βj ,j=1

for some kernel function h. Since the Gaussian likelihood is calculated with ap × p matrix operation in the model ( see, e.g. Wikle, 2010), the computationalburden in the approach data will be greatly reduced, if p << n is chosen. As
long as p < n is designed to be moderate, the model can express flexible class of
nonstationary spatial covariance structures with eﬃcient estimation performances.
Adequate choice of kernel and Σ, which is diﬃcult in practice, would provide goodmodeling with eﬃcient performances. Papers such as Shi and Cressie (2007), Cressie
and Johannesson (2008) are typical ones categorized in this approach employing Wwavelets and local bisquare function as a basis matrix H, respectively. Stroud et.
al. (2007) attempted spatio-temporal modeling in this approach.
The final one is orthogonal basis expansion approach composed of orthogonal
vectors H = (h1 , . . . , hp ) that satisfies h′i hj = δij . The orthogonal basis expansion
for spatial data would provide Gaussian likelihood evaluation that requires one
dimensional operation, if we assume Σ is diagonal. As a result, p may be chosen as
many as sample sizes with suﬃcient approximation accuracy for the mean function

WAVELET ANALYSIS

3

term. The point in this approach is how we construct orthogonal basis for irregularly
spaced data. The traditional orthogonal basis such as Fourier, wavelets, Hermite
polynomials and so on can work only when sampling points are regularly spaced
to constitute mesh data. The orthogonal property of them is usually destructedunder irregular sampling. Paciorek (2007) employed Fourier basis as an orthogonalbasis by binning data locations, the modification of data locations to let them ingrid locations. One principal theme in this paper is how we find orthogonal basis
under irregular sampling without binning.
This paper adopts the final approach for choice of H. Specifically we aim toconstruct block-orthogonal basis H such that H ′ H is block diagonal by modifying
Haar wavelet basis under irregular sampling. Moreover we aim to conduct spatiotemporal extension of the Haar wavelet model from spatial data to spatio-temporal
data. In conducting the extension, we do not assume that sampling spatial points
are the same for each temporal point.
The contributions in this paper are summarized in the following three points.
The first one is construction of block-orthogonal basis matrix H for irregularly
spaced data, which is obtained by modifying Haar wavelets. The second one is
discovering eﬃcient way of estimation of parameters by wavelet transform withconsistency proof under fixed domain asymptotics. The use of block-orthogonalbasis expansion plays a key role in letting Gaussian likelihood evaluation be aneﬃcient form requiring at most three dimensional matrix operations. Final one is aspatio-temporal extension in the practical settings that allow spatial dependency ofparameters to express nonstationary spatial behaviors. Let us begin this paper fromthe first one, the construction of block-orthogonal basis under irregular sampling.
2. Empirical Haar Wavelets
Let sp , p = 1, . . . , n, . . . be possibly countably infinite sampling points for spatial
data, which we suppose are scattered irregularly over [0, 1]2 . For bounded and realvalued functions f and g defined over [0, 1]2 , let us define the inner product by(1)< f, g >:= lim n−1
n→∞

n
∑

f (sp )g(sp ).

p=1

Under the inner product, Fourier basis is no more an orthogonal one except for thecases when sampling points are regularly spaced to provide mesh data.
We construct an orthogonal basis under the inner product. Let us start from theintroduction of the one and two dimensional (1D and 2D) Haar wavelets that constitute an orthogonal basis over [0, 1] and [0, 1]2 , respectively, under the Lebesgue
measure. The mother wavelet over [0, 1] denoted as ψ(x), which is shown in theupper of Figure 1, produces the orthogonal basis for L2 ([0, 1]) given by
ϕk,j (x) := ψ(2k x − j),for k = 0, 1, . . . , j = 0, . . . , 2k − 1, under the Lebesgue measure (see Theorem 1.1
of Ogden (1997)). The index k is called as a resolution parameter, since it controls
the scale of the wavelet. The 1D Haar wavelets satisfy∫
ϕk,j (x)ϕk′ ,j ′ (x)dx = ck δkk′ δjj ′[0,1]

for the positive constant ck = 4−k .

4YASUMASA MATSUDA

1

1

1/2
-1

-1

-1
1

-1

1

1

-1

1

Figure 1. The mother wavelets for 1D and 2D Haar basis.
Similarly the mother wavelets on [0, 1]2 , which are given by the three functions
ψ1 (x), ψ2 (x) and ψ3 (x) defined over [0, 1]2 in lower part of Figure 1, produce theorthogonal basis on L2 ([0, 1]2 ) given by
ϕk,ij (x1 ,x2 ) :=()(2)
ψ1 (2k x1 − i, 2k x2 − j), ψ2 (2k x1 − i, 2k x2 − j), ψ3 (2k x1 − i, 2k x2 − j) ,for k = 0, 1, . . . , i, j = 0, . . . , 2k − 1, under the Lebesgue measure (see Ogden(1997,
page 167)). The index k also is called as a resolution parameter, since it controls
the scale of the 2D wavelet. The 2D Haar wavelets satisfy∫
ϕ′k,ij (x1 , x2 )ϕk′ ,i′ j ′ (x1 , x2 )dx1 dx2 = Ck δkk′ δii′ δjj ′[0,1]2for the diagonal matrix Ck = 4−k I3 .
Let us orthogonalize the 2D Haar wavelets ϕk,ij under the empirical measurein (1) for possibly countably infinite sampling points sp , p = 1, . . . , n, . . . scatteredover [0, 1]2 . Let r(D) be the ratio of the number of the sampling points included in asubregion D ⊂ [0, 1]2 to that for [0, 1]2 . We will orthogonalize ϕk,ij in the following
three disjoint cases specified by the sampling points over its support Dk,ij . Let us
denote the four disjoint equal-area sub-regions in Dk,ij generated from the divisionl, l = 1, 2, 3, 4.
by the lines x1 = (i + 0.5)/2k and x2 = (j + 0.5)/2k as Ek,ij
Definition 1. Let Dk,ij ⊂ [0, 1]2 be the region over which the Haar wavelet ϕk,ij
is defined. We categorize the region Dk,ij by the following three cases.
ll) > 0.
, l = 1, 2, 3, 4 have positive ratios r(Ek,ij
Case 1 All the four subregions Ek,ij
Case 2 Three of the four subregions have positive ratios.
Case 3 Two of the four subregions have positive ratios.

WAVELET ANALYSIS

5

Let us denote the two subregions of Dk,ij over which the mth component of theHaar wavelet ϕk,ij = (ϕ1k,ij , ϕ2k,ij , ϕ3k,ij ) takes values 1 and -1 as Fm and Gm , respectively, for m =values 1 and -1 of the original wavelet√√1, 2, 3. Modify the function
r(G
)/r(F
)
and
−
r(F
)/r(Gm ), respectively, and denote the
ϕm
to
take
m
m
m
k,ij
m
modified function as ϕ̃k,ij for m = 1, 2, 3. Let us define empirical Haar wavelets
ϕ̃k,ij .
Definition 2. The empirical Haar wavelet ϕ̃k,ij is defined by depending on the
category in Definition 1 as:(i) In Case 1, ϕ̃k,ij is defined by the three component function (ϕ̃1k,ij , ϕ̃2k,ij , ϕ̃3k,ij ).
(ii) In Case 2, ϕ̃k,ij is defined by the two component function (ϕ̃1k,ij , ϕ̃2k,ij ),where any two of (ϕ̃1k,ij , ϕ̃2k,ij , ϕ̃3k,ij ) can be chosen as the two components.
(iii) In Case 3, ϕ̃k,ij is defined by the one component function ϕ̃1k,ij , which is
the one that satisfy that r(Fm ) > 0 and r(Gm ) > 0 from among the threecomponents.
(iv) Otherwise ϕ̃k,ij is not defined.
Then a block orthogonality of ϕ̃k,ij in Definition 2 is guaranteed under the empirical measure in (1).
Proposition 1. For possibly countable infinite sampling points of s1 , . . . , sn , . . . ,scattered irregularly over [0, 1]2 , the empirical Haar wavelets ϕ̃k,ij in Definition 2
constitute a block orthogonal basis that satisfiesn∑
lim n−1ϕ̃′k,ij (sp )ϕ̃k′ ,i′ j ′ (sp ) = Ck,ij δkk′ δii′ δjj ′
n→∞

p=1

for the positive definite mk,ij by mk,ij matrix Ck,ij , where mk,ij is the number ofcomponents for ϕ̃k,ij specified in Definition 2.
The proof is given in Section 7.
3. Spatial Data Modeling3.1. Wavelet model for spatial data. We propose a model for spatial data byemploying the empirical Haar wavelets ϕ̃k,ij in Definition 2. Let us consider a modeling of a spatial process Z(s) when sampling is conducted possibly on countably
infinite points sp , p = 1, . . . , n, . . . over [0, 1]2 .
Let Pk be the set of (i, j)s that satisfy Cases 1, 2 or 3 in Definition 1 for a
resolution k. We propose the following model for the spatial data by(3)

Z(sp ) = µ(sp ) + f (sp ) + εp ,where µ is a deterministic mean function, εp is a sequence of independent andidentically distributed normal variables with mean 0 and variance σ 2 and f is therandom process over [0, 1]2 which is defined by
∞
∑
∑
(4)
f (s) =
ϕ̃k,ij (s)βk,ij ,k=0 (i,j)∈Pk

where βk,ij is a sequence of independent normal vectors with mean 0 and varianceC(k + 1)−d Imk,ij

6

YASUMASA MATSUDA

−6

−2

2 4 6

C=5,d=1.2

0.0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

−6

−2

2 4 6

C=5,d=3.6

0.0

0.2

0.4

0.6

Figure 2. Sample paths of the 1D wavelet model in (3) for thecases of (C, d) = (5, 1.2) and (5, 3.6).

for d > 1, where mk,ij is the number of components for ϕ̃k,ij .
The wavelet model has nonstationary covariance function. By simple calculation,
we haveCov (Z(sp ), Z(sq )) =∞∑∑ϕ̃k,ij (sp )ϕ̃′k,ij (sq )C(k + 1)−d + σ 2 δpq ,

k=0 (i,j)∈Pk

which shows nonstationarity by the nature of the Haar wavelets.
In Figure 3.1, we show sample paths of the 1D wavelet model in (3) for C =5, d = 1.2 and C = 5, d = 3.6 when µ = 0 and σ 2 = 0. It is confirmed that larger d
gives a smoother sample path because of more rapid decay of coeﬃcients caused by
the variance C(k + 1)−d . We call the parameter d as smoothness parameter since
it controls the smoothness of the sample path by the wavelet model in (3).
3.2. Whittle likelihood estimation. We consider estimation of the parametersof C, d and σ 2 when the mean function µ(s) is absent in the wavelet model. Suppose we have observed Z(s) that follows the model (3) on finite sampling points
s1 , . . . , sn .
Let ϕ̃k,ij be the n by mk,ij matrix given by (ϕ̃k,ij (s1 ), . . . , ϕ̃k,ij (sn ))′ and Z be
the n by 1 vector given by (Z(s1 ), . . . , Z(sn ))′ . The wavelet transform is definedWAVELET ANALYSIS

7

by the following mk,ij by 1 vector,
(
)−1
(5)
wk,ij = ϕ̃′k,ij ϕ̃k,ij
ϕ̃′k,ij Z,
which is equal to, by substituting the right hand side of (3) for Z and noting thatthe mean function is assumed to be absent,
(
)−1
(6)
βk,ij + ϕ̃′k,ij ϕ̃k,ij
ϕ̃′k,ij ε,where ε = (ε1 , . . . , εn )′ . The block orthogonality of the empirical Haar wavelets
makes the wavelet transform mutually orthogonal with mean 0 and varianceC(k + 1)−d Imk,ij + σ 2 (ϕ̃′k,ij ϕ̃k,ij )−1 ,(7)which makes it easy and eﬃcient to evaluate Gaussian likelihood for the wavelet
transform, because the likelihood requires at most 3 dimensional operation. However the likelihood function is not sensitive to distinguish the first term from thesecond in (7) with respect to the parameters. The likelihood is not expected towork well for estimation. To remedy the ineﬃciency we slightly modify the form of
the variance in (7). Choose two positive integers p < q and redefine the variances
for k = 0, 1, . . . , p as they are, and that for k = q as the one obtained by dropping
the first term in order to distinguish between the two terms. Namely, putting′
Wk,ij = wk,ij
wk,ij ,

Vk,ij (d) = (k + 1)−d Imk,ij ,()−1
Bk,ij = ϕ̃′k,ij ϕ̃k,ij,
we propose the Gaussian likelihood function
L(C, d, σ 2 ) =
p
∑
∑

{
}
log CVk,ij (d) + σ 2 Bk,ij + tr (CVk,ij (d) + σ 2 Bk,ij )−1 Wk,ij

k=0 (i,j)∈Pk

(8)

+

∑

log σ 2 Bq,ij + tr

{(
}
)−1
σ 2 Bq,ij
Wq,ij .

(i,j)∈Pq

The parameters C, d and σ 2 are estimated by the ones that minimize the likelihoodin (8), which we call Whittle estimators in honor of the name for likelihood of
Fourier transform in stationary time series.
3.3. Consistency under fixed domain asymptotics. We show the consistencyof Whittle estimators for the wavelet model in (3) in the asymptotic situation as
the sampling points s1 , . . . , sn tends to infinity over the fixed region [0, 1]2 . Thisis the kind of asymptotics called fixed domain asymptotics (Stein, 1999, page 62),
since sampling points are denser and denser over the fixed region [0, 1]2 as n tends
to be large.
To let the proof for the consistency easier to see, we redefine the parameters in(7) by}{C (k + 1)−d Imk,ij + τ (ϕ̃′k,ij ϕ̃k,ij )−1 .

8

YASUMASA MATSUDA

Concentrating C out of the Whittle likelihood yields the profile likelihoodl(d, τ ) =
log A−1
n

p
∑
∑


∑
{
}
{
}
tr (Vk,ij (d) + τ Bk,ij )−1 Wk,ij +tr (τ Bq,ij )−1 Wq,ij 

k=0 (i,j)∈Pk

(9)
+ A−1
n


p
∑
∑


(i,j)∈Pq∑

log |Vk,ij (d) + τ Bk,ij | +k=0 (i,j)∈Pk

log |τ Bq,ij |(i,j)∈Pq





,

where
An =

p
∑
∑

∑

mk,ij +

k=0 (i,j)∈Pk

mq,ij .

(i,j)∈Pq

The followings are assumed to prove the consistency.
Assumption 1. For the wavelet model in (3), we assume:(1) The wavelet model has no mean function.
(2) The sequences of random variables εp and βk,ij in the wavelet model are
Gaussian.
(3) The parameter spaces for d and τ are D = (1, δ0 ] and Θ = [δ1 , δ2 ], which
include d0 and τ0 as true values inside the regions, respectively, for δ0 > 1
and δ1 > τ0 /2.
(4) There exist sequences of p = pn and q = qn with pn < qn that diverge as
ntends to be large such that()−1
sup pδ0
ϕ̃′p,ij ϕ̃p,ij→ 0,(i,j)∈Pp




∑

−1 
mq,ij 2p∑(k + 1)3/2−δ0k=0(i,j)∈Pq√ ∑
mk,ij  → 0,(i,j)∈Pk

as n tends to be large.
Consistency for the Whittle estimatorsˆ τ̂ ) = arg
(d,
and




Ĉ = A−1
n

p
∑
∑

min

d∈D,τ ∈Θl(d, τ )


{
}
∑
{
}
ˆ + τ̂ Bk,ij )−1 Wk,ij +tr (Vk,ij (d)tr (τ̂ Bq,ij )−1 Wq,ij 

k=0 (i,j)∈Pk(i,j)∈Pq

are shown in the following Theorem 1.
ˆ τ̂ and Ĉ converge in
Theorem 1. Under Assumption 1, Whittle estimators d,
probability to d0 , τ0 and C0 , respectively, as n tends to be infinity, and hence σ̂ 2 =Ĉ τ̂ is also consistent.

WAVELET ANALYSIS

93.4. Kriging. Suppose we have identified the wavelet model in (3) for samples
over s1 , . . . , sn . Let us show the way to estimate the value of an unsampled pointu ∈ [0, 1]2 , which is called traditionally as kriging. Assume that the mean function
is 0 without loss of generality for simplicity. Our model provides an eﬃcient way
for kriging that requires at most 3 dimensional matrix operation.
Bayes theorem applies to kriging in the wavelet model (3), when we regard the
distribution for βk,ij as a prior. The wavelet transform wk,ij is from (6) normally
distributed independently with mean βk,ij and variance σ 2 (ϕ̃′k,ij ϕ̃k,ij )−1 , where theparameter βk,ij has the normal prior with mean 0 and variance C(k + 1)−d Imk,ij .
By Bayes Theorem, the posterior variance and mean are given by
(
)
−1
vk,ij
= C −1 (k + 1)d Imk,ij + ϕ̃′k,ij ϕ̃k,ij /σ 2 ,()
bk,ij = vk,ij ϕ̃′k,ij ϕ̃k,ij wk,ij /σ 2 ,respectively. It follows that the posterior distribution for Z(u) is normal with mean
and variance given by
(10)

∞
∑
∑

ϕ̃k,ij (u)bk,ij ,k=0 (i,j)∈Pk∞∑∑

ϕ̃k,ij (u)vk,ij ϕ̃′k,ij + σ 2 ,

k=0 (i,j)∈Pk

respectively. The mean square error as well as kriging are obtained eﬃciently withat most 3 dimensional matrix operations.

4. Spatio-Temporal Extension
4.1. Wavelet model for spatio-temporal data. We shall extend the wavelet
model in (3) for spatial data to that for spatio-temporal data. Suppose that, for
temporal points t = 1, . . . , T , Z(t, s) is onserved on st,m , m = 1, . . . , nt that are
scattered irregularly over [0, 1]2 . We allow st,p to depend on t to account for thefeature that sampling points may not be the same at each temporal point. Let St
be the set of the sampling points for t = 1, . . . , T and define S0 by ∩Tt=1 St . Let
Pk be the set of (i, j) that is included in any of Cases 1-3 in Definition 1 for theset S0 and k = 0, 1, . . .. For each t, we construct the empirical Haar wavelets for(i, j) ∈ Pk , k=0,1,2,. . . , by the ones that constitute orthogonal basis for s ∈ St .
The empirical Haar wavelets depends on t, since St depends on t. To show thedependency explicitly, we denote the empirical Haar wavelet for each t as ϕ̃k,ij (t, s).
We fit the spatial model in (3) to Z(t, st,m ) for each t and conduct temporal extension by regarding the regression coeﬃcients βk,ij as the state vector that follows
autoregressive models. We propose the state space model for spatio-temporal databyZ(t, st,m ) = µ(t, st,m ) + f (t, st,m ) + εm (t),10

YASUMASA MATSUDA

for t = 1, . . . , T and m = 1, . . . , nt , where µ is a mean function, and(11)f (t, s) =∞∑∑ϕ̃k,ij (t, s)βk,ij (t),k=0 (i,j)∈Pk

βk,ij (t + 1) = ρk,ij βk,ij (t) + uk,ij (t),for t = 1, . . . , T , where uk,ij (t) is a normal random vector with mean 0 and variancematrixAk,ij Imk,ij ,for |ρk,ij | < 1.
The model in (11) can express non-separable space time covariance functions.
By simple calculations, we have(12)

cov (Z(t, u), Z(t − h, v))∑=ϕ̃k,ij (t, u)ϕ̃′k,ij (t − h, v)Ak,ij ρhk,ij + σ 2 δuv δt,t−h ,(k,i,j)∈P0

from which the spatial and temporal covariances are not separated unless ρk,ij is aconstant.
4.2. Parametric modeling of space time covariance functions. Allowing C, din the spatail model in (3) to depend on k, i, j, we identify Ak,ij by(13)

Ak,ij = Ck,ij (k + 1)−dk,ij , d > 1.

Points are in how Ck,ij , dk,ij and ρk,ij are described by parametric functions.
The dependencies of C, d and ρ on k, i, j are interpreted as the local dependenciesof the scale, smoothness and temporal correlation in the space time covariancestructures in (12). It is of crucial importance how the dependencies are identifiedin parametric forms to express locally dependent spatial behaviors.
We shall identify Ck,ij , dk,ij and ρk,ij with functions f1 , f2 and f3 defined over[0, 1]2 through the integral averages
∫
−1
(14)
|Dk,ij |fm (x)dx,
Dk,ij

for m = 1, 2 and 3, respectively, where |Dk,ij | is the area of the domain Dk,ij . To
describe the spatial functions in a parametric form, we employ a Fourier expansion
∑{
}
(15)
fm (x) = a0 +
am,p cos(ωp′ x) + bm,p sin(ωp′ x) ,p∈Q

where ωp is a frequency which we suggest to design by (2πip , 2πjp )′ for the set Q
of mesh points (ip , jp ) except for the origin over the upper half plane of R2 . The
underlying functions f1 , f2 and f3 are interpreted as the locally dependent scale,smoothness and temporal correlation, respectively.
4.3. Whittle likelihood estimation. This section considers estimation for parameters that describe Ak,ij and ρk,ij by the parametric models in (15). To show
explicitly the dependency on parameters, we express the functions as Ak,ij (θ) andρk,ij (θ) for parameters θ ∈ Θ. We assume µ = 0 for simplicity.

WAVELET ANALYSIS

11

Exactly as in (5), we define the wavelet transform by()−1
wk,ij (t) = ϕ̃′k,ij (t)ϕ̃k,ij (t)ϕ̃′k,ij (t)Z(t),for t = 1, . . . , T , where ϕ̃′k,ij (t) = (ϕ̃′k,ij (t, st,1 ), . . . , ϕ̃′k,ij (t, st,nt ))′ and Z(t) =(Z(t, st,1 ), . . . , Z(t, st,nt ))′ .
By the block orthogonality of Φ̃k,ij (t), the state space model for Z(t, st,p ) in (11)reduces to, for (i, j) ∈ Pk , k = 0, 1, . . . , p,(16)

wk,ij (t) = βk,ij (t) + fk,ij (t),βk,ij (t + 1) = ρk,ij (θ)βk,ij (t) + uk,ij (t),where fk,ij (t) is the independent random vector with mean 0 and variance matrix
σ 2 (ϕ̃′k,ij (t)ϕ̃k,ij (t))−1 , which is independent for each k, i, j. It follows that thelikelihood function for the wavelet transform, which we call Whittle likelihood again,
is eﬃciently evaluated by following Kalman recursion.
Let us introduce an explicit form of Kalman recursion to evaluate the likelihood.
Putµk,ij (t) = E(βk,ij (t)|wk,ij (t), . . . , wk,ij (1)),
Qk,ij (t) = var(βk,ij (t)|wk,ij (t), . . . , wk,ij (1)),ηk,ij (t + 1) = E(βk,ij (t + 1)|wk,ij (t), . . . , wk,ij (1)),
Pk,ij (t + 1) = var(βk,ij (t + 1)|wk,ij (t), . . . , wk,ij (1)),
which are evaluated by the Kalman recursion when we initialize by ηk,ij (1) = Omk,ij−1
and Pk,ij(1) = Omk,ij , the diﬀuse prior, namely by
)
(
−1
−2
Q−1
ϕ̃′k,ij (t)ϕ̃k,ij (t) ,k,ij (t) = Pk,ij (t) + σ
)
}
{
(
−1
µk,ij (t) = Qk,ij (t) Pk,ij(t)ηk,ij (t) + σ −2 ϕ̃′k,ij (t)ϕ̃k,ij (t) wk,ij (t) ,(17)

ηk,ij (t + 1) = ρk,ij µk,ij (t),Pk,ij (t + 1) = ρ2k,ij Qk,ij (t) + Ak,ij Imk,ij ,()−1
Fk,ij (t + 1) = Pk,ij (t + 1) + σ 2 ϕ̃′k,ij (t + 1)ϕ̃k,ij (t + 1),

for t = 1, . . . , T and (i, j) ∈ Pk , k = 0, 1, 2, . . . , p.
Define the mk,ij × mk,ij periodogram matrix by′

Jk,ij (t) := (wk,ij (t) − ηk,ij (t, θ)) (wk,ij (t) − ηk,ij (t, θ)) .
Then the likelihood function for the wavelet transform is evaluated as(18)

log Lww (θ) :=

p
T ∑
)
}
∑
∑ { (
−1
tr Jk,ij (t)Fk,ij(t, θ) + log |Fk,ij (t, θ)| ,t=2 k=0 (i,j)∈Pk

which we call the wavelet version of Whittle likelihood for spatio-temporal data.
The parameter θ that describes Ak,ij and ρk,ij by (15) is estimated by minimizing
Whittle likelihood.

12

YASUMASA MATSUDA

4.4. Asymptotic distribution. Asymptotic distribution for Whittle estimators
minimizing (18) is derived when temporal size T tends to be large with finitely
fixed spatial sampling points. Unlike Whittle estimators for the spatial model in(3) under fixed domain asymptotics, traditional theories for maximum likelihood
estimation can apply to have the asymptotic distribution. Assume that spatiotemporal data follows (11) with supk=0,1,2,..., sup(i,j)∈Pk |ρk,ij | < 1. Let θ0 be thetrue value for the parameters θ. Then temporal stationarity of the state vector processes βk,ij (t) in (16) is guaranteed for k = 0, . . . , p, (i, j) ∈ Pk , and the distributionof θ̂ minimizing (18) is approximately, as T tends to be large,
N (θ0 , Ω(θ0 )) ,where(Ω(θ) =

1 ∂ 2 log Lww (θ)
2
∂θ∂θ′

)−1
.

4.5. Forecasting. Following the spatio-temporal data analysis literatures, we denote estimation for future values at t > T as forecasting. Forecasting is conducted
eﬃciently as a result of the Kalman recursion in (17). For h = 1, 2, . . ., definehηk,ij(T ) = ρhk,ij µk,ij (T ),hPk,ij(T ) = ρ2hk,ij Qk,ij (T ) +h−1∑

ρ2a Ak,ij Imk,ij

a=0

and the forecasting for a spatial point u at time T + h is
∞
∑
∑
h
(19)
ϕ̃k,ij (T, u)ηk,ij(T + 1)k=0 (i,j)∈Pk

for which the mean square error is evaluated as∞∑∑h
ϕ̃k,ij (T, u)Pk,ij(T )ϕ̃′k,ij (T, u) + σ 2 .
k=0 (i,j)∈Pk

5. Empirical examples
This section examines the empirical performances of the wavelet modeling byapplying to real data. Public land price data in Tokyo, which is composed of public
land prices for sampling points irregularly scattered over Tokyo, Kanagawa, Chiba
and Saitama prefectures, are employed for the demonstration.
Let us summarize briefly public land price data in Japan. Government of Japan
has been publicizing land prices in the form of price per m2 as of January 1st
over hundred thousands of sampling points scattered irregularly all over Japan.
They evaluate land prices by combinations of real record of transactions, income
approach and cost accounting method, which are therefore not necessarily equal toreal market prices. The purpose of the publications in March every year since 1970
in the government web page is to help evaluate public works, inheritance tax andso on.
This paper focuses on yearly series of land prices in Tokyo areas from 2003 till
2014 to demonstrate eﬃciency of our wavelet approach. The yearly series of land
prices has the following features. First one is large number of sampling points,
which is about 6000-8000 at each temporal point and amounts to 100,000 overall.

WAVELET ANALYSIS

13

year 03
04
05
06
07
08
09
10
11
12
13
no. 8554 8554 8378 8348 8025 7829 7509 7360 6881 6881 6878
Table 1. Number of sampling points for public land prices inTokyo area from 2003 to 2014.

35.8
35.6
35.2

35.4

latitude

36.0

36.2

Sampling points in 2003 are shown in Figure 5, fro which irregularity of sampling
points is confirmed. The second one is that sampling points at each temporal point
are not necessarily the same year by year. Table 1 shows the actual number ofsampling points in each year. The final one is clear exhibition of spatial nonstationarity between behaviors in central Tokyo and in suburbs of Tokyo. The three
features make it impossible to apply existing approaches such as vector ARMA
models in a trivial way and provide us with a good chance to demonstrate our
proposed approach by the empirical Haar wavelet expansions.

139.0

139.2

139.4

139.6

139.8

140.0

140.2

longitude

Figure 3. 8554 sampling points for public land price data inTokyo in 2003.
Let Pt (s) be the land price at spatial point s at temporal point t. Spatio-temporal
data of log return that is the first diﬀerences of log prices, namely(20)

rt (si ) = log Pt (si ) − log Pt−1 (si ), t = 2004, . . . , 2014,is our object for analysis, where si is a ith sampling point that has observations at
two consecutive years t − 1 and t.

14
5977

14

YASUMASA MATSUDA

year

dˆ

σ̂ 2

S/N

wavelet wavelet thresholding√
MSE thrshld=0
log nσ2004 3.01 3.06 8.17
5.47
7.56
6.40
2005 2.43 2.59 6.07
3.34
6.26
4.13
2006 1.36 2.17 3.94
2.15
2.61
3.55
2007 1.23 3.60 5.50
5.61
9.39
9.22
2008 1.15 4.37 5.62
2.69
5.08
5.50
2009 1.65 1.45 2.45
1.31
1.93
3.34
2010 1.46 1.71 2.76
1.31
1.81
3.19
2011 1.43 0.79 1.65
0.67
1.47
1.48
2012 1.78 0.61 2.60
0.59
1.05
0.77
2013 1.46 0.47 2.60
0.47
0.95
0.78
2014 1.72 0.66 2.75
0.66
1.61
2.74
Table 2. Whittle estimatots for d and σ 2 in (3) with the SN
ratios and MSEs of the wavelet kriging in comparison with thoseof wavelet thresholding, where MSEs were evaluated for randomly
chosen 100 points among the samples in each year.

First we fit the wavelet model for spatial data in (3) to the log return in
t = 2004, . . . , 2014. We conduct Whittle estimation for the parameters C, d, σ 2 ,minimizing Whittle likelihood with p = 6, q = 7 in (8), where the samples used forevaluating kriging mean square errors (MSEs) mentioned just below are removed.
With the estimated parameters, we conduct kriging in each year by the formula in(10) and evaluate MSEs of the kriging for randomly chosen 100 points from among
samples. Estimated results as well as MSEs of the krigings are shown in Table 2,
where MSEs of wavelet thresholding (Donoho and Johnstone, 1994) are also shown
as benchmarks. Wavelet thresholding is a kriging just by substituting wk,ij for βk,ijin (3) when |wk,ij | is larger than a threshold√ value and substituting 0 otherwise.
Here the threshold is specified with 0 or log nσ which is motivated by certain
properties of the normal distribution discussed in detail by Donoho and Johnstone(1994). Our method outperforms the benchmarks in terms of MSEs every year.
The SN ratio, which is var(f (s))/var(ε) in (3), declined suddenly in 2009 and has
been held low after Lehman shock. Interpreting var(f (s)) as general spatial tendency reflecting market activity, we find that land market activity in Tokyo becameshrink just after the shock and has not yet recovered yet even in 2014.
Next we fit the wavelet model for spatio-temporal datal in (11) to the yearly
series of log return in (20) separately before and after Lehman shock occurredin September of 2008, namely, separately to the series from 2004 to 2008 and tothat from 2009 to 2013, where samples in 2014 were excluded for use of validating
forecasting performances later. We will focus on the cases when Ck,ij , dk,ij in (13)and ρk,ij that describes the temporal correlations are specified by (14), where weset Q = {(i, j)||i| ≤ 2, 0 ≤ j ≤ 2} ∩ {(i, j)|i > 0 or j > 0} in (15). The identified
model has the 146 parameters for the data size of 75,166 as a result.
We conduct Whittle estimation for the parameters by minimizing the likelihoodin (18) with p = 7. We show the figures of the underlying functions in (14) for ρk,ij
only as the ones that identify interesting features of land prices before and after
Lehman shock. It is observed after Lehman shock that the temporal correlations

WAVELET ANALYSIS

15

declined to be nearly independent in central Tokyo, while they hold high in suburbsof Tokyo.

36.2
0.8
36.0

0.6
35.8

35.6
0.4

35.4
0.2
35.2
139.2

139.4

139.6

139.8

140.0

140.2

36.2
0.8
36.0

0.6
35.8

35.6
0.4

35.4
0.2
35.2
139.2

139.4

139.6

139.8

140.0

140.2

Figure 4. The identified temporal correlation structures before
and afyer Lehman shock, the underlying functions in (14) for ρk,ij .
To see the goodness of fit for the identified spatio-temporal model, we construct
forecasts by the formula in (19) for randomly chosen 500 points from samples in2014 and for points chosen from among them over which temporal correlations havebeen identified as higher than .6 and lower than .3. For comparisons, we consider
a benchmark which forecasts always 0. In Tabke 3, MSEs for both forecasts areshown. It is found that our wavelet method works better especially in highly correlated areas and hence it detects well temporal correlation structures and provides
reasonable forecasts based on the identified correlations.
6. discussion
This paper proposes a model for irregularly spaced data by orthogonal expansion with Haar wavelets and conducts its spatio-temporal extension. The spatiotemporal model can express nonseparable and spatially nonstationary behaviors in16

YASUMASA MATSUDA

no. samples wavelet null forecast
overall
500
2.79
3.65
correl.> .6
165
1.31
2.43
correl.< .3
133
6.33
6.93
Table 3. MSEs of the forecasts by the wavelet model and benchmark that predicts always 0 for randomly chosen 500 points fromsamples in 2014 and for the points from among them over which
temporal correlations are identified with higher than .6 and withlower than .3. The forecastsl are constructed with the model identified by the samples till 2013.

space time covariances by allowing the prior variances of the wavelet coeﬃcientsdependent spatially. The crucial contribution is the discovery of a way to construct orthogonal basis under irregular sampling via Haar wavelets, which makes
our approach unique in the literatures for spatial or spatio-temporal data analysis.
The orthogonal expansion by the empirical Haar wavelets provides us with eﬃcient
ways for kriging or forecasting as well as parameter estimation that requires at mostthree dimensional matrix operations, which opens a way to analyze huge spatiallynonstationary data set of several hundred thousand space time points.
We state some points that are to be studied in the futures. First one is atemporal nonstationary extension of the stationary temporal covariances that thispaper assumed in the nonseparable space time covariances in (12). In order toallow nonstationary in the temporal dimension, there are two possible ways. One
is to apply random walk model instead of the autoregressive model to the statevectors in (11). This approach is expected to be eﬀective for nonstationary datain both space and time when sampling is conducted regularly in time but irregularly in space, which are popular in social science areas. The other one is theuse of three dimensional Haar wavelets to construct an orthogonal basis under irregular sampling over three dimensional space. Then the model by the empirical
three dimensional Haar wavelets describes nonstationary covariances jointly in both
space and time. This approach is expected to be eﬀective for spatio-temporal datacollected irregularly in both space and time, which frequently appears in naturalscience areas.
Second one is possibility of jointly conducting estimation and kriging by fullBayes approach. The models proposed in the paper are Bayesian regression models
whose regression coeﬃcients have prior distributions. Our approach is considered as
an empirical Bayes method that estimates hyper parameters by marginal likelihood
functions, which is Whittle likelihood in our case, and conducts kriging or forecasting separately from the estimation. The jointly conducting approach by full-Bayes
method with MCMC sampling can improve the performances of our empirical Bayes
approach, though it is much more time consuming.
Final one is a multivariate extension of the univariate wavelet model, which we
expect is possible by a combination with the method of coregionalization (Banerjee
et. al. (2004)). It would make it possible to analyze a case when all componentsof multivariate observations are not necessarily observed.

WAVELET ANALYSIS

17

7. proof
7.1. proof of Proposition 1. It is trivial by simple algebra that, for (k, i, j) ̸=(k ′ , i′ , j ′ ),1∑
ϕ̃k,ij (sp )ϕ̃k′ ,i′ j ′ (sp ) = 0.
n→∞ n
p=1
n

lim

l
We will show that Ck,ij is positive definite in Case 1. Put R(Ek,ij) = rl for
l = 1, 2, 3, 4. Let us prove by contradiction. Assume that Ck,ij is singular. Sinceϕ̃k,ij is linearly dependent, there exist a non-zero constants k1 , k2 and k3 that satisfy

1
1
1
+ k2
+ k3
=0
r1 + r2
r1 + r3
r1 + r4
1
1
1
k1
− k2
− k3
=0
r1 + r2
r2 + r4
r2 + r3
1
1
1
−k1
+ k2
− k3
=0
r3 + r4
r1 + r3
r2 + r3
1
1
1
− k2
+ k3
= 0.
−k1
r3 + r4
r2 + r4
r1 + r4
k1

It follows that the matrix given by1
B=

r1 +r2
1
r1 +r21− r3 +r
4

1
r1 +r31− r2 +r
4
1
r1 +r31
r1 +r41− r2 +r3
1− r2 +r3

has a rank that must be less than 3. By elementary row operations, the matrix B
reduces to the upper triangular matrix given by
 1

1
1



r1 +r2

r1 +r3

r1 +r40+r2 +r3 +r4− (rr21 +r
4 )(r1 +r3 )+r2 +r3 +r4− (rr11 +r
4 )(r2 +r3 )00

1 +r2 +r3 +r4 )− 2r(r4 (r1 +r4 )(r2 +r3 ).

r4 must be 0 to let the matrix B be singular, which contradicts the assumptionthat r1 , r2 , r3 , r4 must be all positive.
In Cases 2 and 3, the positive definiteness is proved more easily by the same
argument.
7.2. proof of Theorem 1. We will employ the method by Hurvich et. al (2005)
who proved the consistency of the long memory parameter estimation in long memory time series by by local Whittle estimation.
ˆ Replacing Wk,ij and Wq,ij in (9) with Vk,ij (d0 ) + τ0 Bk,ij
7.2.1. Consistency of d.
and τ0 Bq,ij , respectively, we define l0 (d, τ ) by


p
∑
∑
∑

l0 (d, τ ) = log A−1(k + 1)d−d0 tr {αk,ij (d, τ )} + τ −1 τ0mq,ij nk=0 (i,j)∈Pk

(21)

+ A−1
n


p
∑
∑


k=0 (i,j)∈Pk(i,j)∈Pqlog (k + 1)−d Hk,ij (d, τ ) +∑(i,j)∈Pq

log |τ Bq,ij |





,

18

YASUMASA MATSUDA

where
Hk,ij (d, τ ) = Imk,ij + τ (k + 1)d Bk,ij ,−1αk,ij (d, τ ) = Hk,ij(d, τ )Hk,ij (d0 , τ0 ).

Since (d0 , τ0 ) minimizes l0 ,ˆ τ̂ ) − l0 (d0 , τ0 )0 ≤ l0 (d,ˆ τ̂ ) − l(d,ˆ τ̂ ) + l(d,ˆ τ̂ ) − l(d0 , τ0 ) + l(d0 , τ0 ) − l0 (d0 , τ0 )= l0 (d,≤2|l(d, τ ) − l0 (d, τ )|sup(d,τ )∈D×Θ=2

log 1 +sup(d,τ )∈D×Θp∑∑γk,ij (d, τ ) +k=0 (i,j)∈Pk

∑γq,ij (d, τ ) ,(i,j)∈Pq

where
{
(
)}
(k + 1)d−d0 tr αk,ij (d, τ ) Ek,ij − Imk,ij,
Mn (d, τ )()τ −1 τ0 tr Eq,ij − Imq,ij
γq,ij (d, τ ) =,Mn (d, τ )p∑∑∑(k + 1)d−d0 tr (αk,ij (d, τ )) + τ −1 τ0mq,ij ,Mn (d, τ ) =γk,ij (d, τ ) =k=0 (i,j)∈Pk(i,j)∈Pq−1
Ek,ij = C0−1 k d0 Hk,ij(d0 , τ0 )Wk,ij ,−1
Eq,ij = C0−1 τ0−1 Bq,ij
Wq,ij .

Define

p
∑

∑

∑




mk,ij + τ −1 τ0mq,ij(k + 1)x
k=0(i,j)∈Pk(i,j)∈Pq


p

 ∑
∑
∑
−1
m
.
− A−1
x
log(k
+
1)
m
+
log(τ
τ
)
q,ij
k,ij
0
n



Kn (x, τ −1 τ0 ) = log A−1nk=0(i,j)∈Pk(i,j)∈Pq

The function Kn (x, y) is a non-negative function and twice diﬀerentiable with respect to x with a positive second derivative that is bounded away from 0 for all(d, τ ) ∈ D × Θ, and there exists x0 such that ∂K(x0 , y)/∂x = 0 for a fixed y. Thus
there exist positive constants c1 and c2 such that for (d, τ ) ∈ D × Θ,
∂ 2 Kn (x̃, τ −1 τ0 )(d − d0 − x0 )2∂x2≥ Kn (x0 , τ −1 τ0 ) + c1 (d − d0 − x0 )2 ≥ c2 (d − d0 )2 ,Kn (d − d0 , τ −1 τ0 ) = Kn (x0 , τ −1 τ0 ) +where in the last inequality we used x0 = 0 when τ −1 τ0 = 1 and that Kn (x, y) is
non-negative, continuous and attains 0 if and only if x = 0, y = 1. Hence, defining

WAVELET ANALYSIS

19

Rn (d, τ ) = l0 (d, τ ) − l0 (d0 , τ0 ) − Kn (d − d0 , τ −1 τ0 ), we obtainˆ ≤ c−1 Kn (dˆ − d0 , τ −1 τ0 )0 ≤ (dˆ − d0 )2 1D (d)
2
−1
−1
ˆ
ˆ
≤ c−1
2 l0 (d, τ ) − c2 l0 (d0 , τ0 ) − c2 Rn (d, τ )


p
∑
∑
∑
≤ 2c−1suplog 1 +γk,ij (d, τ ) +γq,ij (d, τ )2(d,τ )∈D×Θ

+ c−1
2

k=0 (i,j)∈Pk(i,j)∈Pq|Rn (d, τ )| ,sup(d,τ )∈D×Θand{Rn (d, τ ) = log

(

∑

∑p

)

d−d0
k=0 (k + 1)(i,j)∈Pk tr αk,ij (d, τ ) − Imk,ij∑∑
1 + ∑pd−d−1 τ00
k=0 (k + 1)(i,j)∈Pk mk,ij + τ(i,j)∈Pq mq,ij

− A−1
n

p
∑
∑

}

log |αk,ij (d, τ )| .

k=0 (i,j)∈Pk

Under Assumption 1. (4), as n tends to be large,suptr(αk,ij (d, τ ) − Imk,ij ) → 0.

sup

k=0,...,p,(i,k)∈Pk (d,τ )∈D×Θ

Hence we obtainsup|Rn (d, τ )| ≤ K(d,τ )∈D×Θ

suptr(αk,ij (d, τ ) − Imk,ij ) → 0.

sup

k=0,...,p,(i,k)∈Pk (d,τ )∈D×Θ

Applying the argument of Proposition A.1 in Hurvich et al. (2005), summation
by parts and Assumtion 1. (4), we obtain that


p
∑
∑
∑
sup
γq,ij (d, τ ) → 0,γk,ij (d, τ ) +log 1 +(d,τ )∈D×Θ

k=0 (i,j)∈Pk(i,j)∈Pqˆ
in probability, by which we have the consistency of d.
7.2.2. Consistency of τ̂ . DefineJn (x) = log A−1
n

p
∑
∑

mk,ij + x

k=0 (i,j)∈Pk

∑

− A−1
n

∑


mq,ij (i,j)∈Pq

mq,ij log x.

(i,j)∈Pq

The function Jn (x) is twice diﬀerentialble, Jn (1) = 0, Jn′ (1) = 0 and Jn′′ (x) isbounded away from 0 for x ∈ [0, 2]. Thus there exists a constant c such that forτ ∈ Θ,()2
Jn (τ −1 τ0 ) ≥ c τ −1 τ0 − 1 .
ˆ τ ) − l0 (d0 , τ0 ) − Jn (τ −1 τ0 ), we obtain
Hence, defining Sn (τ ) = l0 (d,ˆ τ ) − c−1 l0 (d0 , τ0 ) − c−1 Sn (τ̂ )0 ≤ (τ̂ −1 τ0 − 1)2 1Θ (τ̂ ) ≤ c−1 Jn (τ̂ −1 τ0 ) = c−1 l0 (d,
≤ c−1

sup
(d,τ )∈D×Θ

|l0 (d, τ ) − l0 (d0 , τ0 )| + c−1 sup |Sn (τ )| ,
τ ∈Θ

20

YASUMASA MATSUDA

and



Sn (τ ) = log



(
)

ˆ 0
d−d
ˆ
(k
+
1)
tr
α
(
d,
τ
)
− mk,ij k,ij
k=0(i,j)∈Pk
∑
∑p ∑
−1 τ

0
k=0
(i,j)∈Pk mk,ij + τ(i,j)∈Pq mq,ij

∑p
1+

− A−1
n

∑

p
∑
∑

ˆ
ˆ τ) .
log (k + 1)d−d0 αk,ij (d,k=0 (i,j)∈Pk

By the consistency of dˆ and Assumption 1. (4), we obtainsup |Sn (τ )| → 0,τ ∈Θin probability, which proves the consistency of τ̂ .
7.2.3. Consistency of Ĉ. Define−1
ek,ij = C0−1 k d0 Hk,ij(d0 , τ0 )Wk,ij − Imk,ij ,−1eq,ij = C0−1 τ0−1 Bq,ij
Wq,ij − Imq,ij .

ThenĈ − C0 =


p
{
}
∑
∑
∑
ˆ
ˆ τ̂ )ek,ij +k d−d0 tr αk,ij (d,τ̂ −1 τ0 tr (eq,ij )
C0 A−1nk=0 (i,j)∈Pk
+C0 A−1np∑∑(i,j)∈Pq

k

ˆ 0
d−d


{
}
∑
( −1
)
ˆ τ̂ ) − Imtr αk,ij (d,+mq,ij τ̂ τ0 − 1 k,ij

k=0 (i,j)∈Pk(i,j)∈Pq→0in probability.
References[1] Bai, Y., Song, P. and Raghunathan, T. E. (2012) Joint composite estimating functions inspatio temporal models. J. R. Statist. Soc. B, 74, 799-824.
[2] Banerjee, S., Carlin, B. and Gelfand A. (2004) Hierarchical Modeling and Analysis for Spatial
Data. Chapman & Hall/CRC.
[3] Banerjee, S., Gelfand A., Finley, A. O. and Sang, H. (2008) Gaussian predictive process models
for large spatial data sets. J. R. Statist. Soc. B, 70, 825–848.
[4] Brockwell, P. J. and Davis, R. A. (1991) Time Series: Theory and Methods. Springer-Verlag,
2nd ed.
[5] Cressie, N. and Johannesson, G. (2008) Fixed rank kriging for very large spatial data sets. J.
R. Statist. Soc. B , 70, 209-226.
[6] Donoho, D. L. and Johnstone, I. M. (1994) Ideal spatial adaptation by wavelet shrinkage.
Biometrika, 81, 425-455.
[7] Hurvich, C., Moulines, E. and Soulier, P. (2005) Estimating long memory in volatility.
Econometrika, 73, 1283-1328.
[8] Kaufman, C. G., Schervish, M. J. and Nychka, W. (2008) Covariance tapering for likelihoodbased estimation in large spatial data sets. J. Am. Statist. Assoc., 103, 1545–55.
[9] Ogden, R. T. (1997) Essential Wavelets for Statistical Applications and Data Analysis.
BirkHäuser, Boston.
[10] Paciorek, C. J. (2007) Bayesian smoothing with Gaussian processes using Fourier basis functions in the spectral GP package. J. Statist. Software, 19, 1-28.
[11] Robinson, P. M. (1995) Gaussian semiparametric estimation of long range dependence. Annalsof Statistics, 24, 1630-1661.

WAVELET ANALYSIS

21[12] Shi, T. and Cressie, N. (2007) Global statistical analysis of MISR aerosol data: a massive
data product from NASA’s Terra satellite. Environmetrics, 18, 665-680.
[13] Stein, M. L. (1999) Interpolation of Spatial Data. New York: Springer-Verlag.
[14] Stroud, J., Mueller and Sanso, B. (2001) Dynamic models for spatial-temporal data. J. Royal
Statist. Soc. B., 63, 673-689.
[15] Wikle, C. K. (2010) Low-rank representations for spatial processes. Handbook of Spatial
Statistics. (Edited by A. E. Gelfand, P. J. Diggle, M. Fuentes and P. Guttorp), 107-118.
Chapman & Hall/CRC.
[16] Zimmerman, D. L. and Stein, M. (2010) Classical geostatistical methods. Handbook of Spatial
Statistics. (Edited by A. E. Gelfand, P. J. Diggle, M. Fuentes and P. Guttorp), 29-44. Chapman& Hall/CRC.
Graduate School of Economics and Management, Tohoku University, 27-1 Kawauchi,
Aoba ward, Sendai 980-8576, Japan
E-mail address: matsuda@econ.tohoku.ac.jp