Discussion Paper No.311
WAVELET ANALYSISOFSPATIO-TEMPORAL DATA
Yasumasa Matsuda
January 28, 2014

TOHOKU ECONOMICS RESEARCH GROUPGRADUATE SCHOOL OF ECONOMICS AND

MANAGEMENT TOHOKU UNIVERSITY

27-1 KAWAUCHI, AOBA-KU, SENDAI,
980-8576 JAPANWAVELET ANALYSIS OF SPATIO-TEMPORAL DATAYASUMASA MATSUDA

Abstract. This paper aims to provide a wavelet analysis for spatio-temporal
data which are observed on irregularly spaced stations at discrete time points,
where the spatial covariances show serious non-stationarity caused by local
dependency. A specific example that is used for the demonstration is US precipitation data observed on about ten thousand stations in every month. By
a reinterpretation of Whittle likelihood function for stationary time series, we
propose a kind of Bayesian regression model for spatial data whose regressors
are given by modified Haar wavelets and try a spatio-temporal extension bya state space approach. We also propose an empirical Bayes estimation forthe parameters, which is regarded as a spatio-temporal extension of Whittle
likelihood estimation originally defined for stationary time series. We conduct the extended Whittle estimate and compare mean square errors of theforecasts with those of some benchmarks to evaluate its goodness for the US
precipitation data in August from 1987-1997.

1. Introduction
This paper focuses on analysis of spatio-temporal data, which is observations ofregular time series at huge number of stations scattered irregularly over space. We
aim at proposing a method of modelling, estimation and kriging for the kind of
spatio-temporal data whose spatial covariances are not necessarily stationary.
Several kinds of models for space time covariances have been proposed to analyze
spatio-temporal data. The simplest one is a separable covariance that is givenby the product of spatial and temporal covariances, which makes it possible togive a covariance model separately in space and time. It provides an easy way
to identify models for spatio-temporal data by fitting the product of the temporaland spatial covariances given by time series models such as autoregressive andmoving average (ARMA) models (Brockwell and Davis, 1991) and spatial models
such as Matérn class (Banerjee et al., 2004), respectively. Separable covariances,
however, restrict covariance structures in the very narrow range in which temporalcorrelations on each spatial point must coincide. Nonseparable covariances canprovide a practical class of space time covariances, but require a careful treatment
to guarantee the positive definiteness of the covariance functions. Gneiting (2002)
provides a suﬃcient condition for the positive definiteness and propose a usefulclass of space time covariances that satisfy it.
After a choice of space time covariance models, estimation of the parametersin the space time covariance model is a next step. Since the Gaussian maximum
likelihood requires calculations of determinant and inverse of the covariance matrix
whose dimension is the number of space time points of the data, MLE is in practicediﬃcult to apply especially for huge data set. In order to avoid this computational
Key words and phrases. empirical Bayes estimate. Haar wavelet. Kalman filter. kriging.
periodogram. spectral density function. state space model. Whittle likelihood function.
1

2

YASUMASA MATSUDA

diﬃculty, several methods have been proposed in the literatures. Typical ones that
represents them are the covariance tapering by Kaufman et. al. (2008), predictive
process approach by Banerjee et. al. (2008) and the composite function approach by
Bai et. al. (2012). The composite approach, which is regarded as the combined oneof the other two methods, gave the most satisfactory performances in their empiricalstudies, while the number of space time points that the predictive approach cancope with is the largest of the three methods. In fact, the numbers of space time
points that the three papers analyzed in their empirical studies were about 15,000,
7,000 and 3,000, respectively, and we see that as many as several thousand spacetime points are supposed in the composite approach, which is regarded as the one
to represent the existing methods.
In this paper, we have an interest in analysis of huge data set observed at stations scattered irregularly over so broad space such as all over US continent thatstationary spatial covariances are not necessarily satisfactory models. Specifically,
our main object is to consider a method of modelling spatio-temporal data over as
many as several hundred thousand space time points that shows crucial nonstationarity in spatial covariances, for which growing attentions have been paid recently
in areas of both natural and social sciences by rapid progress of data collection
technologies. The existing methods, which basically assume stationary covariances
over as many as several thousand space time points, do not suppose the kind of
spatio-temporal data as their objects of modelling.
This paper proposes a modelling by wavelets, which is an approach completely
diﬀerent from the existing methodology of covariance model fitting by maximizing
a likelihood approximated in the clever ways to avoid large dimensional matrix
operations. We will employ a kind of Bayesian regression model with regressors
given by wavelets whose regression coeﬃcients have prior distributions, which isconsidered as an extension of Fourier analysis of time series (Brockwell and Davis,
1991). The use of wavelets was originally proposed for analysis of nonstationary
time series by Nason et. al. (2000). This paper extends the wavelet model fromtime series to spatio-temporal data by using the Haar wavelets that are modified
to let them empirically orthogonalized under irregular sampling.
The striking features of the wavelet model for spatio-temporal data are the following three points. The first one is that it can express nonseparable space timecovariances given with nonstationary spatial covariances and spatially dependent
temporal covariances. The second one is that the use of empirically orthogonalized Haar wavelet makes it possible to calculate the likelihood function, kriging
and forecasting eﬃciently by just three dimensional matrix operations in Kalman
recursions, which means it opens a way of analysis for huge data set larger thanseveral hundred thousand. The final one is that it allows any numbers of NAs byregarding the data as samples of continuous stochastic process expanded by themodified Haar wavelets.
US precipitation data observed at 11918 stations scattered over US continent
in August from 1987 to 1997 is a typical example that has the features we focus
in this paper. The space time points are about 70,000 with clear appearance of
crucial nonstationarity caused by local dependency of rain falls and the existence of
frequent NAs. We fit the wavelet models to the data and examine the goodness of fit
by performances of kriging and forecasting in order to demonstrate our methodology
by the wavelet modelling.

WAVELET ANALYSIS

3

2. Fourier analysis of stationary time series
Let {Xt }, t ∈ Z be stationary time series with the covariance function
γ(θ, h) = Cov(Xt , Xt−h ),for a parameter θ ∈ Θ. We assume that the spectral density function f (θ, λ) thatsatisfies∫ π(1)γ(θ, h) =exp(ihλ)f (θ, λ)dλ,−πexists.
To estimate the parameter θ ∈ Θ by the observations X = (X1 , . . . , Xn )′ , weusually use the maximum likelihood estimator (MLE) θ̂ that maximizes the log
likelihood function
1
1
log L(θ) = − |Γn (θ)| − X ′ Γn (θ)−1 X,
22where Γn (θ) is the n by n covariance matrix whose (i, j)th element is given byγ(θ, i − j). Since the calculation of the log likelihood function is sometimes time
consuming especially for large sample sizes because of large dimensionality of thecovariance matrix, we often use an approximated likelihood called Whittle likelihood function given by}[n/2] {
1 ∑I(ωk )
log Lw (θ) = −(2)
log f (θ, ωk ) +,2
f (θ, ωk )k=1

where ωk = 2πk/n, which is called Fourier frequency, and I(λ) is the periodogram
given by the squared discrete Fourier transform (DFT) on λ, namely by
n

(3)

2

1 ∑
I(λ) =
Xt exp(−iλt) .
2πn t=1

We will give an interpretation for the Whittle likelihood function to show intuitively how it approximates the exact likelihood function. Though many timeseries researchers have examined the asymptotic equivalence between the estimators that maximize log L(θ) and log Lw (θ) (see e.g. Proposition 10.8.3 of Brockwell
and Davis, 1991), the theories do not directly provide how the exact likelihood isrelated with the Whittle likelihood.
Let us consider what model the Whittle likelihood function is the exact likelihood
function for. For independent random variables αk , k = −n/2, . . . , n/2 with mean
0 and variance f (θ, ωk ), let us consider the following model
√
[n/2]
2π ∑
(4)
αk exp(iωk t),X̃t =n
k=−[n/2]

which is regarded as a Bayesian regression model whose regression coeﬃcients havethe prior distribution with no error terms. Since the regressors exp(iωk t) for k =−[n/2], . . . , [n/2] constitute an orthonormal basis that satisfies
1∑
exp(iωp t)exp(iωq t) = δpq ,
n t=1
n

4

YASUMASA MATSUDA

for Kronecker delta δpq , the DFT of X̃t given by
√
n
1 ∑
(5)
X̃t exp(−iωk t)
2πn t=1
reduces to αk . It follows naturally that the likelihood function for X̃t is equal to
the Whittle likelihood function (2) by considering the likelihood for DFT instead ofX̃t . In other words, the Whittle likelihood in (2) for Xt is regarded as the likelihoodfor X̃t in (4), and the use of the Whittle likelihood function for MLE is equivalent
to approximating Xt by X̃t .
The approximation of Xt by X̃t is justified in the sense that they have asymptotically equivalent covariance functions. In fact, by simple calculation, we have
(6)

2π
Cov(X̃t , X̃t−h ) =
n

∑

[n/2]

exp(ihωk )f (θ, ωk ),k=−[n/2]

which is a Riemannian approximation for γ(θ, h) in (1).
The frequency domain expression in (4) for time series, which we derive from aninterpretation for the Whittle likelihood, plays a significant role for extension fromtime series models to those for spatio-temporal data.
3. Extension to spatio-temporal data
This section proposes a model for spatio-temporal data by an extension of thefrequency domain expression for time series in (4). The crucial point for the extension is to find an orthogonal basis for irregularly spaced data points, as theorthogonality of the sinusoidal functions for time series plays a significant role inderiving the Whittle likelihood.
3.1. Modified Haar wavelets. For the extension of the model in (4) from timeseries to spatio-temporal data, it is necessary to find a new orthogonal basis under
the inner product defined by observation stations for spatio-temporal data. Let
sj , j = 1, . . . , n be the locations of observation stations of spatio-temporal data,
which we suppose are scattered irregularly over [0, 1]2 . For two functions f and g
defined on sj , let us define the inner product by(7)< f, g >:=n∑f (sj )g(sj ).

j=1

The sinusoidal function exp(iωs) no more constitutes an orthogonal basis under
the inner product except for the cases when the stations are regularly spaced toprovide a mesh data.
We construct an orthogonal basis under the inner product by modifying the twodimensional Haar wavelets. Let us start from the introduction of the one and two
dimensional 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] denotedas ψ(x), which is shown in the upper of Figure 1, produces the orthogonal basis on
L2 ([0, 1]) given by
ϕk,j (x) := ψ(2k x − j),WAVELET ANALYSIS

5

1

1

1/2
-1

-1

-1
1

-1

1

1

-1

1

Figure 1. The Haar wavelets on [0, 1] and on [0, 1]2 .
for k = 0, 1, . . . , j = 0, . . . , 2k − 1, under the Lebesgue measure (see Theorem 1.1 of
Ogden (1997)). In other words, the one dimensional Haar wavelets satisfy∫
ϕk,j (x)ϕk′ ,j ′ (x)dx = ck,j δkk′ δjj ′[0,1]

for the positive constant ck,j = 4−k .
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,i,j (x1 ,x2 ) :=()(8)
ψ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)). In other words, the two dimensional Haar wavelets satisfy∫
ϕ′k,ij (x1 , x2 )ϕk′ ,i′ j ′ (x1 , x2 )dx1 dx2 = Ck,ij δkk′ δii′ δjj ′[0,1]2for the positive definite diagonal matrix Ck,ij = 4−k I3 .
Let us modify the two dimensional Haar wavelets to let them constitute theorthogonal basis under the empirical measure in (7) for irregularly spaced stations
sj , j = 1, . . . , n scattered over [0, 1]2 . Let N (D) be the number of the stations
included in D ⊂ [0, 1]2 , and Dk,ij ⊂ [0, 1]2 be the support for the Haar wavelet
ϕk,ij , which is give by
Dk,ij = {(x1 , x2 )|i/2k ≤ x1 ≤ (i + 1)/2k , j/2k ≤ x2 ≤ (j + 1)/2k } ⊂ [0, 1]2 .
lDefinition 1. In the subregions Fk,ij, Glk,ij in Dk,ij over which the lth componentof the Haar wavelet ϕk,ij takes values 1 and -1, we modify the values 1 and -1 to6those given by

YASUMASA MATSUDA

v (
v (
)
)
u
u
u N Glu N Fl
k,ij
k,ij
u
u
) and − t (),
t (l
N Fk,ij
N Glk,ijlrespectively, which we denote as ψ̃k,ij(x1 , x2 ) for l = 1, 2, 3. Then the modified Haar
wavelet is defined by
(
)
1
2
3
(9)
ϕ̃k,ij (x1 , x2 ) = ψ̃k,ij(x1 , x2 ), ψ̃k,ij(x1 , x2 ), ψ̃k,ij(x1 , x2 ) .

By simple calculation, it is confirmed that the modified Haar wavelets constitutea block orthogonal basis under the empirical measure in (7), namely, satisfyn∑(10)
ϕ̃′k,ij (sp )ϕ̃k′ ,i′ j ′ (sp ) = Lk,ij δkk′ δii′ δjj ′ ,p=1

where Lk,ij is not a diagonal matrix except for mesh data cases, which means thatthe three components are not necessarily mutually orthogonal, and is not always
positive definite. Since the positive definiteness of Lk,ij is crucial for modelling of
spatio-temporal data, we introduce a suﬃcient condition for it.
Here we show a suﬃcient condition to guarantee the positive definiteness of
3 × 3 matrix Lk,ij in (10), or equivalently the linear independence of the three n()′
dimensional vectors given by ϕ̃k,ij = ϕ̃′k,ij (s1 ), . . . , ϕ̃′k,ij (sn ) . Let us denote thefour disjoint equal-area sub-regions in Dk,ij generated from the division of Dk,ij byl
the lines x1 = (i + 0.5)/2k and x2 = (j + 0.5)/2k as Ek,ij, l = 1, 2, 3, 4.
Proposition 1. For observation stations of s1 , . . . , sn scattered irregularly over[0, 1]2 , let ϕ̃k,ij (s), k = 0, 1, . . . , i, j = 0. . . . , 2k − 1 be the modified Haar wavelets.
A suﬃcient condition for the 3 × 3 matrix Lk,ij to be positive definite is that eachlof the four disjoint sub-regions Ek,ij, l = 1, 2, 3, 4 in Dk,ij contains at least one
observation station.
The proof is given in Section 6. Let P be the set of (k, i, j)s that satisfy thesuﬃcient condition for the positive definiteness of Lk,ij . In the following subsections, the block orthogonal basis ϕ̃k,ij for (k, i, j) ∈ P plays a role in spatial andspatio-temporal extension of the frequency domain expression in (4) for time series.
3.2. Wavelet model for spatio-temporal data. We shall extend the frequency
domain expression in (4) for time series to that for spatial and moreover spatiotemporal models by employing the modified Haar wavelets in (8) that constitutethe block orthogonal basis under the empirical measure in (7). Let us start from themodelling of spatial data where observation stations are located in sp , p = 1, . . . , nscattered irregularly over [0, 1]2 , for which we define P as the set of (k, i, j)s thatguarantee the positive definiteness of Lk,ij in (10) by Proposition 1. In analogy
with the frequency domain expression in (4) for time series, we extend it to spatial
data model by∑(11)
Z(sp ) =ϕ̃k,ij (sp )βk,ij + εp ,(k,i,j)∈P

where εp is independent and identically distributed random variables with mean
0 and variance σ 2 , which we call nugget eﬀect, and βk,ij is the three dimensional

WAVELET ANALYSIS

7

independent random vector with mean 0 and variance matrix Ak,ij I3 for the 3 × 3
identity matrix I3 , which is a wavelet version of spectral density in analogy withthe time series expression in (4).
Unlike the stationary covariances given by (6) for (4), our proposed model in(11) has nonstationary covariances. In fact, by simple calculations, the covariancefunction is evaluated as∑
Cov (Z(sp ), Z(sq )) =ϕ̃k,ij (sp )ϕ̃′k,ij (sq )Ak,ij + σ 2 δpq ,(k,i,j)∈P

which shows nonstationarity by the nature of the Haar wavelets.
Next we shall conduct a spatio-temporal extension from the spatial model in(11) by state space models. Suppose that , for time points t = 1, . . . , T , stations
are located in st,p , p = 1, . . . , nt that are scattered irregularly over [0, 1]2 , on which
we observe Z(t, st,p ). We allow st,p to depend on time to account for an existenceof NA. Let ϕ̃k,ij (t, s) be the modified Haar wavelets defined for the stations attime t, and Pt for the set of (k, i, j) that satisfy the suﬃcient condition for the
positive definiteness of Lk,ij in (10). Let P0 = ∩Tt=1 Pt . For each t, we propose thespatial model in (11) to Z(t, st,p ) and conduct temporal extension by regarding theregression coeﬃcients βk,ij as the state vector that follows autoregressive models.
Let βk,ij (1) be the initial state vectors that are three dimensional independent
random vectors with mean 0 and variance matrix given by Ak,ij I3 , we propose thestate space models for spatio-temporal data by∑(12)
ϕ̃k,ij (t, st,p )βk,ij (t) + εp (t),Z(t, st,p ) =(k,i,j)∈P0

βk,ij (t + 1) = ρk,ij βk,ij (t) + uk,ij (t),for t = 1, . . . , T , where uk,ij (t) is the three dimensional independent random vectors
with mean 0 and variance matrix given by(1 − ρ2k,ij )Ak,ij I3 ,and we assume that |ρk,ij | < 1.
The model in (12) can express non-separable space time covariance functions.
In fact, by simple calculations, we have(13)

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.
3.3. Parametric modelling of sapce time covariance functions. We shallpropose a specific parametric form for the functions Ak,ij and ρk,ij in the wavelet
model (12). Let Dk,ij be the support for the modified Haar wavelet ϕ̃k,ij . Since
Ak,ij and ρk,ij are the quantities that describe the strength for the wavelet ϕ̃k,ij
defined over Dk,ij , they can be regarded as mappings from Dk,ij to R. The point
for the modelling is how we find the parametric mappings to describe the spacetime covariance function in (13).

8

YASUMASA MATSUDA

A suﬃcient condition for (13) to be well defined by the convergent series is that∞∑k=0

max |Ak,ij | < ∞,i,jand |ρk,ij | < 1. As a model that satisfies the suﬃcient condition, we propose, for
0 < τk,ij < 1 and Bk,ij that is positive and uniformly bounded,kAk,ij = Bk,ij (1 − τk,ij )τk,ij.

From the functional form in (13) for h = 0, it can be seen that Bk,ij and τk,ij specify
a magnitude and a smoothness for spatial covariance structures, respectively, while
ρk,ij describes temporal correlation structures.
We shall identify Bk,ij , τk,ij and ρk,ij with the functions f1 , f2 and f3 definedover [0, 1]2 through the integral averages∫(14)|Dk,ij |−1
fm (x)dx,
Dk,ij

for m = 1, 2 and 3, respectively, where |Dk,ij | is the area of the domain Dk,ij . The
underlying function fm (x) is given by a Fourier expansion
∑{
}
(15)
fm (x) = a0 +
am,p cos(ωp′ x) + bm,p sin(ωp′ x) ,p∈Q

where ωp is a frequency (2πip /C1 , 2πjp /C2 )′ for the set Q of mesh points (ip , jp )except for the origin over the upper half plane of R2 and positive constants C1
and C2 . It follows that the Fourier expansion for each of the underlying functions
provides each of Bk,ij , τk,ij and ρk,ij with the 2|Q|+1 dimensional parametric form,
which will give a specific parametric form for the space time covariance function asa result.
3.4. Whittle likelihood estimate. This section considers estimation for a parameter that describes the scalar functions Ak,ij and ρk,ij in the spatio-temporalmodel in (12), when we observe in the stations st,p , t = 1, : . . . , T, p = 1, . . . , nt .
To show explicitly the dependency on the parameter, we express the functions as
Ak,ij (θ) and ρk,ij (θ) for the parameter θ ∈ Θ.
In analogy with the discrete Fourier transform for time series, let us define thewavelet transform for the spatial data Z(t) = Z(t, st,1 ), . . . , Z(t, st,nt ))′ at time t.
()′
Let Φ̃k,ij (t) be the nt × 3 matrix given by ϕ̃′k,ij (t, st,1 ), . . . , ϕ′k,ij (t, st,nt ) . Thenthe wavelet transform (WT) is defined by the ordinary least square (OLS) estimator
for the coeﬃcient on the regressor Φ̃k,ij (t) for (k, i, j) ∈ P0 , namely by
(
)−1
(16)
wk,ij (t) = Φ̃′k,ij (t)Φ̃k,ij (t)Φ̃′k,ij (t)Z(t),which is the 3 dimensional vector that corresponds with the OLS estimator for
βk,ij (t) in (12).
By the block orthogonality of Φ̃k,ij (t) shown in (10), the state space model forZ(t, st,p ) in (12) reduces to, for (k, i, j) ∈ P0 ,(17)

wk,ij (t) = βk,ij (t) + fk,ij (t),βk,ij (t + 1) = ρk,ij (θ)βk,ij (t) + uk,ij (t),WAVELET ANALYSIS

9

where fk,ij (t) is the independent three dimensional random vector with mean 0 andvariance matrix σ 2 (Φ̃′k,ij (t)Φ̃k,ij (t))−1 .
Let us evaluate the likelihood function for the WT that follows the state space
expression in (17) by Kalman recursion. Putµk,ij (t) = E(βk,ij (t)|wk,ij (t), . . . , wk,ij (1)),
Qk,ij (t) = V ar(βk,ij (t)|wk,ij (t), . . . , wk,ij (1)),
mk,ij (t + 1) = E(βk,ij (t + 1)|wk,ij (t), . . . , wk,ij (1)),
Pk,ij (t + 1) = V ar(βk,ij (t + 1)|wk,ij (t), . . . , wk,ij (1)),
which are evaluated by the Kalman recursion when we initialize by mk,ij (1) = 0
and Pk,ij (1) = Ak,ij I3 , namely byvk,ij (t) = wk,ij (t) − mk,ij (t),()−1
Fk,ij (t) = Pk,ij (t) + σ 2 Φ̃′k,ij (t)Φ̃k,ij (t)
,
(18)

−1
µk,ij (t) = mk,ij (t) + Pk,ij (t)Fk,ij(t)vk,ij (t),−1
Qk,ij (t) = Pk,ij (t) − Pk,ij (t)Fk,ij(t)Pk,ij (t),mk,ij (t + 1) = ρk,ij µk,ij (t),Pk,ij (t + 1) = ρ2k,ij Qk,ij (t) + (1 − ρ2k,ij )Ak,ij I3 ,for t = 1, . . . , T and (k, i, j) ∈ P0 .
Define the 3 × 3 periodogram matrix in analogy with (3) for time series by′Ik,ij (t) := (wk,ij (t) − mk,ij (t, θ)) (wk,ij (t) − mk,ij (t, θ)) .
Then the likelihood function for the WT is evaluated asT)}1∑ ∑ { (−1(19) log Lww (θ) := −tr Ik,ij (t)Fk,ij(t, θ) + log |Fk,ij (t, θ)| ,2 t=1(k,i,j)∈P0

which we call the wavelet version of Whittle likelihood function in analogy with
(2) for time series. The parameter θ that describes the functions Ak,ij and ρk,ij in(12) is estimated by maximizing the Whittle likelihood function, which provides an
eﬃcient algorithm conducted just by the three dimensional matrix operations.
3.5. Kriging and forecasting. Following the spatio-temporal data analysis literatures, we define estimation for the values on unknown points at 1 ≤ t ≤ T as
kriging and that for the values on any points at t > T as forecasting. Both kriging
and forecasting can be conducted eﬃciently as a result of the Kalman filtering in(18). The kriging for a spatial point u at time 1 ≤ t ≤ T is given by∑
ϕ̃k,ij (t, u)µk,ij (t)(k,i,j)∈P0

for which the mean squared error is evaluated as∑ϕ̃k,ij (t, u)Qk,ij (t)ϕ̃′k,ij (t, u) + σ 2 ,(k,i,j)∈P0

while the forecasting for a spatial point u at time T + 1 is given by∑
ϕ̃k,ij (T, u)mk,ij (T + 1)
(k,i,j)∈P0

10YASUMASA MATSUDA

for which the mean squared error is evaluated as∑ϕ̃k,ij (T, u)Pk,ij (T + 1)ϕ̃′k,ij (T, u) + σ 2 .
(k,i,j)∈P0

3.6. Mean function estimation. We have been considering the case when themean function is assumed to be 0 for spatio-temporal data. Here we will consider
the case when the mean function is not ignored. Let g(t, s) be a mean function attime t and spatial point s. We will consider the case when the mean function given
by a regression form is inserted in the spatio-temporal model in (12), which is given
by
∑
(20)
Z(t, st,p ) = g(t, st,p ) +ϕ̃k,ij (t, st,p )β(t) + εt,p ,(k,i,j)∈P0g(t, st,p ) =q∑

Xi (t, st,p )αi (t),i=1for independent variables Xi (t, s) and temporally dependent regression coeﬃcients
αi (t), i = 1, . . . , q.
First, we assume that the case when Ak,ij and σ 2 are known in (20). Suppose
we aim to estimate the coeﬃcients α(t0 ) = (α1 (t0 ), . . . , αq (t0 ))′ at a temporal point
t0 , 1 ≤ t0 ≤ T . Let α̂(t0 ) be the initial estimate for α(t0 ), which is typically the
OLS
α̂(t0 ) = {X ′ (t0 )X(t0 )}−1

X ′ (t0 )Z(t0 ).

Follow the recursion:(1) Calculate the wavelet transform in (16) at t0 for the residual data given by(21)

Z̃(t0 , st0 ,p ) = Z(t0 , st0 ,p ) −q∑

Xi (t0 , st0 ,p )α̂i (t0 ),i=1

for p = 1, . . . , nt0 .
(2) Evaluate µk,ij (t0 ) in the recursion (18) when we initialize the Kalman filter
as mk,ij (t0 − 1) = 0 and Pk,ij (t0 − 1) = Ak,ij .
(3) Calculate the residual process by∑Ẑ(t0 , st0 ,p ) = Z(t0 , st0 ,p ) −ϕ̃k,ij (t0 , st0 ,p )µk,ij (t0 )(k,i,j)∈P0

for p = 1, . . . , nt0 .
(4) Estimate α(t0 ) by the OLS for the residualα̂(t0 ) = {X ′ (t0 )X(t0 )}−1

X ′ (t0 )Ẑ(t0 ).

(5) Return to (1).
Conducting the recursion until it converges, we have the best linear unbiased estimator for α(t0 ) given by the converged one under the covariance structures identified
by Ak,ij and σ 2 , which we call the wavelet version of generalized least square (GLS)
estimator.
In the practical situation when Ak,ij and σ 2 are unknown, we need to estimate
them in prior to estimate the mean function. For the residual with the OLS forα(t0 ), apply the Whittle likelihood estimation to the OLS residual to obtain theestimators for θ and σ 2 . And replacing the parameters with the estimators in theWAVELET ANALYSIS

11

recursive procedures, we estimate the mean function by following the recursion withthe estimated parameters until it converges.
4. Empirical examples

latitude

25

30

35

40

45

This section examines the empirical performances of the wavelet modelling proposed in Section 3 by applying to real and simulated data. First, we apply it to
US precipitation data, which are observations over about 70,000 space time points
that shows crucial nonstationarity in the spatial covariances, to catch the empiricalproperties. Next, we conduct simulation studies to confirm whether they still holdgenerally for simulated samples generated by the exact model that we identify inthe real analysis.

−120

−110

−100

−90

−80

−70

longitude

Figure 2. The locations of 11918 stations for US precipitation data.
4.1. Applications to real data. Let us start from real data analysis by the
wavelet methods. We focus on US precipitation data, which are observations of
the total monthly precipitation collected at 11918 stations scattered irregularly allover US, and will examine the empirical performances of the wavelet modelling byusing specifically the yearly observations in August from 1987 to 1997.
The data are available in the web site:
http://www.image.ucar.edu/Data/US.monthly.met/USmonthlyMet.shtml.
The locations of the observation stations are shown in Figure 2, in which it is foundthat they are irregularly scattered all over US continent. All of the stations do notnecessarily conduct observation every month. Table 1 shows the numbers of thestations that have the observations in August from 1987 to 1997, from which wefind that NAs happened frequently.
We fit the wavelet model in (20) to the yearly precipitation in August for theten years from 1987 to 1996, where we take constant and height of the stations as

12

YASUMASA MATSUDA

year 87
88
89
90
91
92
93
94
95
96
no. 7042 7118 7176 7079 7009 7055 7043 6927 6823 6716
Table 1. The numbers of stations that have the observations for
US monthly precipitation in August from 1987 to 1997.

97
6747

regressors for the mean function. Based on the fitted model, we conduct kriging
and forecasting for some spatial points in 1997 to examine the goodness of fit.
Parametric model in (14) for Ak,ij and ρk,ij in the wavelet model (12) must be
specified for application to the yearly data in August. We will employ the specificcase when we identify Ak,ij by
Ak,ij = Bk,ij (1 − τ )τ k
for a constant τ that does not depend on k, i, j. It means that the smoothness forthe spatial covariances is assumed to be constant all over the US continent, while
the magnitudes of them and temporal correlations may be spatially dependent by
spatial dependencies of B and ρ, respectively And for the specific modelling of thespatial dependent Bk,ij and ρk,ij , we set C1 = C2 = 1.2 with Q = {(i, j)||i| ≤ 2, 0 ≤j ≤ 2} ∩ {(i, j)|i > 0 or j > 0} in (15), which means that we employ the modelwith the 50 parameters for the data size of 69,988 as a result.
We conduct Whittle estimation in (19) to estimate the parameters in the model
when the mean function estimated by OLS is deleted from the original data in eachyear from 1987 to 1996. And with the estimated parameters, we conduct kriging
and forecasting for some points in 1997 when we estimate the mean function by therecursive procedure in section 3.6. For the locations of the points for kriging andforecasting, we chose locations with high temporal correlations as well as randomly
scattered points to see the forecasting performances in relations with the temporalcorrelations.
As a result of Whittle estimation, τ and σ 2 are estimated as 0.83 and 6.34,
respectively. Since the 24 parameters to describe each of Bk,ij and ρk,ij are too
many to show the estimated values here, we summarise them in Figure 3, whichshows the figures of the underlying functions in (14) for Bk,ij and ρk,ij identified
with the estimated parameters. Table 2 describes the averages of the squared errorsof kriging and forecasting for 100 points randomly chosen from 6747 stations in 1997
and those of forecasting for 125 points in 1997 that have the temporal correlations
identified as being larger than 0.95 by the underlying function for ρk,ij . Note thatkriging and forecasting were calculated by the model estimated with the data until
1996, as if the selected points in 1997 for kriging and forecasting were unknown.
As the benchmarks for comparisons, the weighted averages of the data in 1997 and1996 for kriging and forecasting, respectively, are calculated by the normal kernel
with the bandwidths given by 0.1, 1, 2 and ∞, where the one with the bandwidth∞ means just the sample average of the data.
Figure 3, which are the identified figures of the underlying functions for Bk,ijand ρk,ij , detects well the spatial dependencies of the spatial covariances and thetemporal correlations. And Table 2 shows that the wavelet model has the best
performances in kriging and forecasting except for the one case, which means theWAVELET ANALYSIS

13

25
45

20
40

15
35
10
30
5

25
−120

−110

−100

−90

−80

−70

1.0

45

0.8

0.6
40

0.4
35
0.2

30
0.0

−0.2

25
−120

−110

−100

−90

−80

−70

Figure 3. The identified figures of the underlying functions in(14) for the magnitudes Bk,ij and temporal correlations ρk,ij in
the upper and lower, respectively.

waveletweighted average
bandwidth
0.1
12∞
kriging
6.61 6.75
7.07
8.27 18.08
forecast 1
20.23 67.91 42.53 36.41 18.22
forecast 26.65 14.38 13.73 12.52
9.20
Table 2. MSEs of the kriging and forecasting 1 for 100 randomly
chosen points in 1997 and those of the forecasting 2 for 125 points
that have the temporal correlations identified as being larger than
0.95.

reasonable goodness of fit for the identified wavelet model. The temporal correlation at the location for forecasting has crucial eﬀects for the forecasting performance. Precisely, the forecasting performances of the weighted averages dependon the choice of the bandwidth in relations with the temporal correlations, while
the wavelet method adjusts automatically to the temporal correlation to provide a
better forecast.

14

YASUMASA MATSUDA

waveletweighted average
0.1
1
2
∞
8.15
14.59
16.23
27.27
32.51
(1.10)(2.04)(2.49)(7.16)(9.05)
forecast 1
29.58
47.60
36.03
33.94
32.50
(7.68) (10.56) (9.41)(9.52)(9.05)
forecast 2
15.58
21.15
19.21
21.83
22.90
(9.91) (12.84) (11.93) (16.50) (19.29)
Table 3. Averages and standard deviations in parenthesis of
MSEs of kriging and forecasting for the data points in the 11th
year, where the locations of the points for forecasting and kriging
were selected to correspond exactly with those used in the real
analysis in Table 2.
bandwidth
kriging

4.2. Simulation studies. In order to confirm the observations obtained throughthe applications to US precipitation data, we will conduct simulation studies by
the exact model identified in the real analysis.
We simulate 100 sets of spatio-temporal data for 11 years by the estimated modelin the real analysis when the locations of the stations are designed to be exactly
the same as those of original data. We did not include the mean function in thesimulation model for simplicity. By the first 10 years data set, we conduct Whittle
likelihood estimation in (19) for the parameters as if they were unknown and compare the estimators with the true values to evaluate the estimation performances.
And for the last 11th year spatial data set, we conduct kriging and forecasting
based on the estimated model to see their performances in comparisons with thebenchmarks by the weighted averages.
We will summarize the results of the estimation, kriging and forecasting for the
100 sets of simulated data to show the performances as simple as possible. Figure
4 shows the figures of median, 5% and 95% points of the estimated underlying
functions for B and ρ with the true values as functions of latitudes when longitudeis fixed as 36.5. Table 3 shows the averages and standard deviations of the meansquared errors of kriging and forecasting for the 100 points and those of forecastingfor the 125 points, where the locations of the points for forecasting and kriging wereselected to correspond exactly with those used in the real analysis.
Figure 4 demonstrates the reasonable performances of Whittle estimators. The
estimators have almost no bias and the true values are included within the 90%
intervals empirically evaluated by the 100 estimators. Table 3 confirms the observations that we stated through the results in Table 2. The averages of MSEs for
all the cases in Table 3 are similar to those of Table 2, which are all included in theintervals within one standard deviations. The wavelet method automatically adjusts to the temporal correlation in the point for forecasting to provide reasonablekriging and forecasting, while the weighted averages require careful choice of the
bandwidth in relations with the temporal correlation that is usually unknown.
5. discussion
This paper propose a model for spatio-temporal data by an interpretation of
Whittle likelihood function in stationary time series. The striking feature of theWAVELET ANALYSIS

15

10

15

20

25

magnitude

−120

−110

−100

−90

−80

longitude

−0.2

0.0

0.2

0.4

0.6

0.8

temporal correlation

−120

−110

−100

−90

−80

longitude

Figure 4. The true values and 5%, 95% and 50% points of theestimators as functions of longitudes with the latitude being 36.5.
The upper and lower figures show those of the underlying functions
for the magnitude and temporal correlation, respectively.

proposed model is that it opens a way for analysis of nonstationary huge spatiotemporal data set larger than several hundred thousand by the Kalman recursion
for Whittle likelihood estimation, kriging and forecasting.
We state some points that are to be studied in the futures. First one is thenonstationary extension of the stationary temporal covariances that depend just on
the time lag h in the nonseparable space time covariances in (13). In order to allow
nonstationary also in the temporal dimension, there are two possible ways. One is to
apply random walk model instead of the autoregressive model to the state vectors in(12). This approach is expected to be eﬀective for nonstationary data in both space
and time when the sampling is conducted regularly in time but irregularly in spacesuch as yearly collection of land price data at many locations, which are popular insocial science areas. The other one is the use of three dimensional Haar wavelets
to construct an orthogonal basis under irregular sampling over three dimensional
space. Then the model by the modified three dimensional Haar wavelets directly
describe nonstationary covariances not only in space but also in time not via thestate space temporal extension. This approach is expected to be eﬀective for spatiotemporal data collected irregularly in both space and time, which are frequently
observed in natural science areas.

16

YASUMASA MATSUDA

Second one is the use of continuous wavelets other than the Haar wavelets.
Since the Haar wavelets are used as a basis for our modelling, kriged and forecast
values model are necessarily piecewise constants. An extension of Haar wavelets tocontinuous ones that can keep the orthogonality under irregular sampling must be
conducted to have continuous estimates. Presently the extension is possible only
in mesh data cases.
Third one is a multivariate extension of our wavelet model, which we expect ispossible by the method of coregionalization (Banerjee et. al. (2004)). It requires
intensive considerations in the space time points where a part of components ofobservations are missing,
Finally, it is a challenging problem to construct asymptotic theories for the
wavelet models in the fixed domain asymptotics (see Stein (1999)), when the lengthof time is fixed. The problem reduces to the very simple one to estimate theparameters c and δ in independent samples Xk , k = 1, . . . , n with mean 0 andvariance cδ k for 0 < δ < 1. The diﬃculty lies in the asymptotic decay of thevariance as n tends to infinity.
6. proofl
In this section, we will prove Proposition 1. Put N (Ek,ij) = nl for l = 1, 2, 3, 4.
Let us prove by contradiction. Assume that Lk,ij is singular. Since ϕ̃k,ij is linearly
dependent, there exist the non-zero constants k1 , k2 and k3 that satisfy
1
1
1
k1
+ k2
+ k3
=0
n1 + n2
n1 + n3
n1 + n4
1
1
1
k1
− k2
− k3
=0
n1 + n2
n2 + n4
n2 + n3
1
1
1
−k1
+ k2
− k3
=0
n3 + n4
n1 + n3
n2 + n3
1
1
1
−k1
− k2
+ k3
= 0.
n3 + n4
n2 + n4
n1 + n4
It follows that the matrix given by


1
1
1

B=

n1 +n21
n1 +n2
1− n3 +n4

n1 +n31− n2 +n
4
1
n1 +n3n1 +n41
− n2 +n3
1− n2 +n3

has the 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



n1 +n2n1 +n3n1 +n40+n2 +n3 +n4− (nn21 +n
4 )(n1 +n3 )+n2 +n3 +n4− (nn11 +n
4 )(n2 +n3 )00

4 (n1 +n2 +n3 +n4 )− 2n(n1 +n4 )(n2 +n3 ).

n4 must be 0 to let the matrix B be singular, which contradicts the assumptionthat n1 , n2 , n3 , n4 must be all larger than 1.
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.

WAVELET ANALYSIS

17[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, 2nd edn. SpringerVerlag.
[5] Gneiting, T. (2002) Nonseparable, stationary covariance functions for space-time data. J. Am.
Statist. Assoc., 97, 590-600.
[6] 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-1555.
[7] Nason, G., von Sachs, R. and Kroisandt, G. (2000) Wavelet processes and adaptive estimationof the evolutionary wavelet spectrum. J. R. Statist. Soc. B., 62, 281-292.
[8] Ogden, R. T. (1997) Essential Wavelets for Statistical Applications and Data Analysis.
BirkHäuser, Boston.
[9] Stein, M. L. (1999) Interpolation of Spatial Data. New York: Springer-Verlag.
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