Discussion Paper No. 124
Convolutional regression
for big spatial data
Yasumasa Matsuda and Xin Yuan

February, 2022

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Convolutional regression
for big spatial data
Yasumasa Matsuda and Xin Yuan*

Abstract
Recently it is common to collect big spatial data on a nationalor continental scale at discrete time points. This paper aims at aregression model when both dependent and independent variables are
big spatial data. Regarding spatial data as functions over a region,
we propose a functional regression by a parametric convolution kernel
together with the least squares estimation on the frequency domain
by applying Fourier transform. It can handle massive datasets withasymptotic validations under the mixed asymptotics. The regression
is applied to Covid-19 weekly new cases and human mobility collectedin city levels all over Japan to find that an increase of human mobility
is followed by an increase of Covid-19 new cases in time lag of two
weeks.

keywords: CARMA kernels, Convolution, Fourier transform, Irregularly
spaced data, Least squares estimation, NTT DoCoMo spatial statistics, Periodogram, Spurious regression.

1

Introduction

With rapid technological innovations, data structures become increasingly
complex and big data analysis becomes especially important. Recently, it
is common to collect big spatial data on a national or continental scale atdiscrete time points. For example, NTT DoCoMo, Inc., a Japanese mobile
phone company, provides big spatial dataset recording number of people
staying at 500 meter meshes all over Japan every one hour, resulting interabytes scale spatial data.
Regarding a big spatial data recorded at a time point as one sample ofdependent or independent variable in regression analysis, we aim at a regression model in this paper that identifies a non-trivial relationship among big
spatial data. For example, it is an interesting problem how much increaseof human mobility promote Covid-19 infections, which will be carefully examined in the later section.
*

Graduate School of Economics and Management, Tohoku University

1Let yt (s) and xt1 (s), . . . , xtp (s), t = 1, . . . , T , be discrete time observations of big spatial data for s ∈ S ⊂ R2 , a large scale region such as nationor continent. Regarding them at each time point as a functional data, weconsider a regression model,yt (s) = Φ1 (xt1 )(s) + · · · + Φp (xtp )(s) + εt (s), t = 1, . . . , T,
where Φj is not a constant in usual regression but a linear operator on L2 (S).
This paper specifies it with a convolution, which is given byZφj (s − u)xt (u)du,Φj (xt )(s) =S

where φj (·) is a function on S called a convolution kernel. We call a regression for functional data with a convolution kernel as convolutional regression.
Convolution kernel is a special case of Hilbert Schmidt kernel defined byZΨj (xt )(s) =ψj (s, u)xt (u)du,S

with the relation of φ(s − u) = ψ(s, u). A convolution kernel at s puts moreweights on closer points to s. Convolutional regression exploits this relation.
Existing studies in the fields of regression for functional data are foundin the literatures of functional data analysis. Ramsay and Silverman (2005)
introduces functional regression together with the estimation approach offunctional principal component analysis (fPCA). It is regarded as a nonparametric approach to identify Hilbert Schmidt kernels with empiricallyorthogonal functions. There have been plenty of studies on functional regression under varieties of settings. See, for functional regression, Yao et al.
(2005), Horvath and Kokoszka (2012), Crambes and Mas (2013), Benatia etal. (2017) and so on, while see, for functional auto-regression, Bosq (2000).
Liu et al. (2016) proposes the unique approach for functional auto-regression
with a convolutional kernel by a sieve estimation, which is another nonparametric method different from fPCA. Baseline of them supposes a regression
or auto-regression with Hilbert-Schmidt kernels, and identifies the kernels
by fPCA under the condition of stationarity across time without restrictionsof stationarity across space. They basically suppose an application to functional data on [a, b], a fixed compact interval on R, although it is claimed
that they can be extended to S ⊂ R2 .
This paper tries an approach for spatial regression on huge scales as
an alternative to existing ones on fPCA from the two viewpoints. One isthat our approach employs parametric identification for convolution kernels rather than nonparametric one for Hilbert Schmidt kernels. The otherone is that Fourier transform works to estimate convolution kernels in our
approach, while empirically orthogonal functions by fPCA does to identify
Hilbert Schmidt kernels in existing approaches. It means that our approach2requires no singular value decompositions but Fourier basis functions instead. Fourier transform restricts the applications to spatial data that must
be stationary across space but need not be across time, while fPCR restricts
spatial data that must be stationary across time but need not be acrossspace. As a result, our approach can handle spatio-temporal data nonstationary across time from short to large time periods including pure spatial
data without temporal points. Besides our approach does not suffer spurious
regression problems indicated by Granger and Newbold (1974), which will
be demonstrated in our simulation studies.
One feature of our approach to be stated finally is the asymptotic validation of our estimation procedures on the frequency domain by Fouriertransform, which is obtained by extending the frequency domain techniquesdeveloped by time series researchers such as Brockwell and Davis (1991,
chap. 10), Robinson (1995), Hosoya (1997) and so on, from time series tospatial data. Under the asymptotic scheme called the mixed domain aymptotics (Stein, 1999, pp. 62), the consistency and asymptotic normality are
obtained under the condition of stationarity across space. The consistentestimation of the asymptotic covariance matrix, easily evaluated in practice, makes it possible to conduct statistical inference by t tests in the same
way with usual regression analysis. It is actually demonstrated in the empirical section that the statistical inference is successfully applied to detect
two week lagged effects between human mobility and Covid-19 new cases
collected on national scale in Japan.

2

Convolutional regression for spatial data

We introduce a regression model for spatial data on national or continental
scales. Convolutional kernels are employed to define a regression between
spatial data.

2.1

Parametric convolution kernels

Let L2 (S) be the set of square integrable functions on a region S ⊂ R2 .
We regard spatial data on S ⊂ R2 at discrete time of t = 1, . . . , T as realizations of L2 (S) valued random variables, denoted as yt (s) and xt (s) =(xt1 (s), . . . , xtp (s))0 , dependent and independent variables, respectively.
For f (s), a L2 (S) function, a convolution is an operator that transforms
f (s) to L2 (S) function byZ
φ(s − u)f (u)du,S

where φ(s) is a convolution kernel. A convolution is a special case of Hilbert3Schmidt operator given byZ
ψ(s, u)f (u)du,S

where ψ(s, u) is a Hilbert Schmidt kernel.
Applying φj (s; θ), a parametric convolution kernel, to xtj (s), we definea convolutional regression given byyt (s) =

p Z
X
j=1

φj (s − u; θ)xtj (u)du + εt (s), s ∈ S,

(1)

S

where εt (s) is an unobserved error term regarded as L2 (S) valued randomvariables.
Parametric convolutional kernels are employed in our approach rather
than nonparametric Hilbert Schmidt kernels in traditional approaches offunctional regression. Traditional approaches mainly target a regression
for L2 [a, b], one dimensional functions on a compact region. Although it isclaimed that they are extended to that for L2 (S), SsubsetR2 , the extension is
not straightforward especially for big spatial data on national or continental
scale.
The spatial regression in (1) allows yt−1 (s), . . . , yt−k (s), temporally lagged
terms of dependent variable, to be included as independent variables. Specifically, convolutional auto-regression is defined in the same way as
yt (s) =

k Z
X
j=1

ψj (s − u; θ)yt−j (u)du +

S

p Z
X
j=1

φj (s − u; θ)xtj (u)du + εt (s).

S

Conditions to confirm stationary across time can be derived by the arguments in Bosq (2000). We omit details here as stationarity across time is notnecessary to validate the estimation in Section 3. We include convolutional
auto-regression in (1) throughout this paper. Estimation and asymptotic
validations in the later sections work for both convolutional regression andauto-regression.
Let us comment on intercept terms in the convolutional regression. We
can include a parametric intercept term c(s; θ) in (1). However, each of independent variables should be adjusted to be zero mean to avoid identification
issues indicated by Kokoszka and Reimherr (2017, pp 75).

2.2

Least squares estimation on the frequency domain

This section considers estimation of parameters that describe the parametricconvolutional kernels in (1) in practical situations when sampling points
are irregularly spaced and may not be identical for each time point. We

4will propose the least squares estimation on the frequency domain after
introducing the one on the spatial domain as the equivalent alternative.
The usual least squares estimation of θ is obtained by minimizing
Qsp (θ) =T Z
Xyt (s) −

t=1

p Z
X

2

ψa (s − u; θ)xta (u)du ds,

(2)

a=1

which we call is the least squares estimation on the spatial domain.
Let us consider least squares estimation on the frequency domain as analternative to the one on the spatial domain, both of which are equivalentmathematically. Fourier transform of f (s) ∈ L2 (R2 ) is defined byZ0˜f (s)e−is ω ds, ω ∈ R2 .
f (ω) =
R2

It is well known that Fourier transform of a convolutionZg(s) =
ψa (s − u)f (u)duR2

is given byg̃(ω) = ψ̃a (ω)f˜(ω).
Hence, applying the Fourier transform to (1), we obtainỹt (ω) =pXψ̃j (ω; θ)x̃tj (ω) + ε̃t (ω), ω ∈ R2 ,

(3)

j=1

and we estimate the parameter θ by minimizing

Qf req (θ) =T ZXỹt (ω) −

t=1

p
X

2

ψ̃j (ω; θ)x̃tj (ω) dω,

(4)

j=1

which we call the least squares estimation on the frequency domain.
Since Qsp (θ) = Qf req (θ) by Parseval’s theorem, the least squares estimation on both spatial and frequency domains are completely equivalent.
Practical situations when we evaluate them by several approximations, however, makes the two least squares estimation be different. This paper will
focus the least squares on the frequency domain by two reasons. One is lesscomputational costs of the convolution caused by Fourier transform. The
other one is the asymptotic validations in Section 3 that make it possible toconduct statistical inference.
Suppose we observe yt (s) and xta (s) on randomly distributed points
inside a rectangular A = [0, A1 ] × [0, A2 ], where observation points may
or may not be identical for t = 1, . . . , T . Let us denote them as s0tj , j =5

1, . . . , n0t for yt (s) and satj , j = 1, . . . , nat for xta (s). Then let us define thediscrete Fourier transform (DFT) of yt (s) and xt (s) by, for t = 1, . . . , T ,
n

0t
1 X
ŷt (ω) =yt (s0tj ) exp(−iω 0 s0tj ),
n0t

j=1

x̂ta (ω) =

nat
1 X

nat

xta (satj ) exp(−iω 0 satj ), a = 1, . . . , p.

j=1

Then replacing ỹt (ω), x̃ta (ω) in (3) with the DFTs on mesh points of Fourier
frequencies,
2πf1 2πf2
ωf =,, (f1 , f2 ) ∈ Z2 ,A1
A2
in a compact region D on R2 , we estimate the parameter θ by minimizingQ̃f req (θ) =T X
Xŷt (ωf ) −t=1 ωf ∈DpX2ψ̃a (ωf ; θ)x̂ta (ωf ) ,

(5)

a=1

which we denote as θ̂. Q̃f req (θ) is regarded as an empirical evaluation of
Qf req (θ) = Qsp (θ). We propose the one minimizing Q̃f req (θ) as the estimator in the practical situations when sampling points x and y are irregularly
spaced and may or may not be identical for each t with possibly time varying
sample sizes.
Remark 1. Let us introduce the empirically adjusted version of Qsp (θ)
in (2). Approximating the convolution with Riemannian summation, we
estimate θ by minimizing, for observation points satj , j = 1, . . . , Nat , a =0, 1, . . . , p,
Q̃sp (θ) =

N0t
T X
X
t=1 j=1

)2
p
Nat
X
1 X
ψa (s0tj − satk ; θ)xat (satk ) ,yt (s0tj ) −
Katj

(

a=1

k=1

(6)
with
Katj =NatX

φ0 (s0tj − satk ; θ).

k=1

Although Qsp (θ) = Qf req (θ) theoretically, their empirical counterparts ofQ̃sp (θ) and Q̃f req (θ) are not equal because of several differences in the approximations. It follows that the two least squares estimation on the spatial
and frequency domains are not the same. We claim advantages of the estimation on the frequency domain for big spatial data on huge scales that6

come from lower computational costs and the asymptotic validations in thefollowing section.
Remark 2. As in usual regression analysis, we often need to constructpredicted values. Approximating the convolution in (1) by Riemannian summation, we evaluate the forecast at v and t by
p
Nat
X
1 X
ψ(v − satk ; θ̂)xat (satk ),ŷt (v) =
Kat
a=1

Kat =

Nat
X

(7)

k=1

ψ(v − satk ; θ̂).

k=1

2.3

Parametric family of CARMA kernels

Suitable parametric family of kernels are necessary to apply the convolutional regression in (1) to real data in practice. They need to satisfy two
requirements. First, they cover reasonably broad class of convolutions by aparametric family. Second, they have the closed form expressions of Fouriertransforms necessary to estimate the parameters by the least squares on thefrequency domain in (5). CARMA kernels, which were proposed in Brockwell and Matsuda (2017) to define random fields as an extension of Ornstein
Uhlenbeck processes
Q (Uhlenbeck and Ornstein, 1930), are introduced.
Q
Let a(z) = pi=1 (z 2 − λ2i ) with Re(λi ) < 0 and b(z) = qi=1 (z 2 − ξi2 )with 0 ≤ q < p and λ2i 6= ξj2 for all i and j. CARMA(p, q) kernel is definedbyφ(s) =pXb(λi )i=1

a0 (λi )exp(λi ||s||), s ∈ R2 ,where a0 denotes the derivative of the polynomial a and ||s|| is the Euclidean
norm of he vector s. CARMA kernel is an isotropic function of ||s|| and has
the closed form expressions of Fourier transform given byφ̃(ω) = −pX

λi b(λi ).
0a (λi )(||ω||2 + λ2i )3/2
i=1

The simplest kernel of CAR(1), obtained by substituting p = 1, q = 0
in the CARMA(p, q) kernel, and the Fourier transform are given by, after
transforming λ1 with −1/δ for later use,
β −||s||/δe, δ > 0, s ∈ R2 ,2πδ 2
= βφ0 (s; δ), say,βφ̃(ω; β, δ) = 3, ω ∈ R2 ,
δ (||ω||2 + δ −2 )3/2
φ(s; β, δ) == β φ̃0 (ω; δ), say.
7

(8)

(9)

The CAR(1) kernel here is normalized to let the integral over R2 be β andmake easier interpretation for the parameter β as a result. The convolutional
regression with the CAR(1) kernels is expressed aspXZ

1 −||s−u||/δj
xt (u)du + εt (s), δ > 0, s ∈ S,2e
S 2πδjj=1
Z
pX
φ0 (s − u; δj )xt (u)du + εt (s).
=βjyt (s) =

j=1

βj

(10)

S

The parameter β measures a degree of spatial association between x and y,while δ describes spatial effects that measure how far the neighbours of x(s)affect y(s). As δj → 0, it reduces toyt (s) =pXβj xtj (s) + εt (s), s ∈ S,

j=1

the usual regression model without spatial effects. One point to be noticed
here is that β = 0 in (8) makes δ unidentifiable. Testing spatial independence
between x(s) and y(s) by β = 0 should be tried for a fixed value of δ. In thelater simulation studies, we tried the test of independence by testing β = 0
for δ = 0.

3

Asymptotic properties

This section considers asymptotic properties of the least squares estimator in(5) when we have samples of yt (s), xta (s), a = 1, . . . , p, on irregularly spaced
points on [0, A1 ] × [0, A2 ] for t = 1, . . . , T . Three kinds of asymptotics are
known in spatial statistics literatures (Stein, 1999). This paper employs
the mixed asymptotic among the three when both region of [0, A1 ] × [0, A2 ]and number of sampling points diverge jointly in the spatial dimension,
while the temporal size of T is finite and fixed in the temporal dimension.
The techniques to prove the asymptotic results are based on Matsuda andYajima (2009) and Matsuda and Yajima (2018), which have developed themethod to manage irregularly spaced data under the mixed asymptotics.
Let us denote the least squares estimator minimizing (5) as θ̂ throughout
this section.

3.1

Consistency and asymptotic normality

Let us clarify the conditions under which consistency and asymptotic normality of our estimator is validated.
C1. The sample sizes nat and the sampling region A = [0, A1 ] × [0, A2 ]diverge jointly such that A1 → ∞, A2 → ∞, A1 /A2 , A2 /A1 = O(1),8nat /nbt = 0(1) and |A|/nat → 0, a, b = 0, 1, . . . , p, t = 1, . . . , T , for thearea |A| = A1 × A2 . We shall employ a suffix k such as A = Ak , nat =natk when we indicate explicitly that they diverge as k tends to infinity.
C2. Let Sat be the set of sampling points in A = [0, A1 ]×[0, A2 ]. We assume
that elements in Sat are written as, for a = 0, 1, . . . , p, t = 1, . . . , T ,satj = (A1 u1,atj , A2 u2,atj ), j = 1, . . . , nat ,where uatj = (u1,atj , u2,atj )0 is a sequence of independent and identically distributed random vectors with a probability density function
g(x) supported on [0, 1]2 which has continuous first derivatives. Sampling points s and observed variables xt , εt at s are independent for allt.
C3. xt (s) = (xt,1 (s), . . . , xt,p (s))0 has finite moments of all orders for any
t = 1, . . . , T , s ∈ R2 , and is stationary across s ∈ R2 with all higher
order spectral density functions that are all bounded and integrable.
The 2nd order spectral density function,Z
0−2
fx,t1 t2 ,a1 a2 (ω) = (2π)
Cov(Xt1 ,a1 (s), Xt2 ,a2 (s − u))e−iu ω du,R2

is positive definite and twice differentiable.
C4. xt (s) is strictly exogenous in the sense that xt1 (s1 ) is independent of
εt2 (s2 ) for any t1 , t2 = 1, . . . , T, s1 , s2 ∈ R2 . εt (s) has finite momentsof all orders, and is stationary across s ∈ R2 with all higher order
spectral density functions that are all bounded and integrable. The
2nd order spectral density,Z0fε,t1 t2 (ω) = (2π)−2Cov(εt1 (s), εt2 (s − u))e−iu ω du,R2

is positive and twice differentiable.
C4’. xt (s) is contemporaneously exogenous in the sense that xt (s1 ) is independent of εt (s2 ) for any t = 1, . . . , T, s1 , s2 ∈ R2 . εt (s) has finite
moments of all orders, is serially independent across t and stationary
across s ∈ R2 with all higher order spectral density functions. The
2nd order spectral density function,Z
0−2fε,t (ω) = (2π)
Cov(εt (s), εt (s − u))e−iu ω du,R2

is positive definite and twice differentiable.

9C5. Let Θ be a compact parameter space in Rq for θ = (θ1 , . . . , θq ), and D
be a fixed symmetric compact region on R2 . ψ̃(ω; θ) has a continuous
first derivative with respect to θ on Θ × D. θ1 6= θ2 implies thatψ̃(ω; θ1 ) 6= ψ̃(ω; θ2 ) on a subset of D with positive Lebesgue measure.
C1 clarifies the details of the mixed asymptotics under which the asymptotic
properties of the estimator is validated. Notice that the asymptotics in spaceis employed with fixed and finite temporal size T rather than joint asymptotics of space and time. C2 specifies irregularity of sampling points over thediverging region A = [0, A1 ] × [0, A2 ] with the identical density function of|A|−1 g(s/A). The distribution of sampling points needs to be independentof dependent and independent variables. It should be emphasized that C3and C4 (C4’) confine zt (s) = (yt (s), xt (s)) to be stationary across s without
stationary restrictions across t. Namely, zt (s) can be nonstationary acrosst. C4 assumes the strict exogeneity of xt (s) to allow serial correlation of theerror terms, while C4’ assumes serially independency of the error terms torelax the strict exogeneity with the contemporaneous exogeneity. C5 is theusual assumption of identifiability of parametric convolutional kernels.
Now we are in a position to state the consistency and asymptotic normality. The proofs are shown in the final section of appendix.
Theorem 1. Under Assumptions C2, C3, C4 and C5,θ̂ →p θ0in the asymptotic scheme defined by C1. Assumption C4 can be replacedwith C4’ to have the consistency.
Theorem 2. Under Assumptions C2, C3, C4 and C5, if |A|3/2 /nat =o(1), a = 0, 1, . . . , p,p

|A|(θ̂ − θ0 ) → N 0, (2π)2 bg Ω−1 ΣΩ−1in the asymptotic scheme defined by C1, where
(Z
) (Z
)−2
|g(u)|4 dubg =|g(u)|2 du[0,1]2ZΩ=Ψ̇(ω; θ0 )DZΣ=Ψ̇(ω; θ0 )D,[0,1]2( TX)
fx,tt (ω) Ψ̇0 (ω; θ0 )dω,

t=1
T
1 X
2

fε,t1 t2 (ω)fx,t1 t2 (ω)

Ψ̇0 (ω; θ0 )dω,

t1 ,t2 =1

and where fx,tt (ω) is the spectral density matrix whose (a, b)th element is
fx,tt,ab (ω) in Assumption C3 and Ψ̇(ω; θ) is the q by p matrix whose (i, j)th10

element is∂ ψ̃j (ω; θ), i = 1, . . . , q, j = 1, . . . , p.
∂θi
If Assumption C4’ is assigned for Assumption C4, Σ in the asymptotic variance matrix reduces to
( T
)
Z
1X
Σ=
Ψ̇(ω; θ0 )fε,tt (ω)fx,tt (ω) Ψ̇0 (ω; θ0 )dω.
2
D
t=1

We will give a few remarks on the asymptotic results. First, the asymptotic results are validated under the asymptotics when observation regionand samples inside together diverge for a fixed time period, which makes
it possible to relax stationarity across time under stationarity only acrossspace. Secondly, bg is the term that clarifies the effects of distributions ofsampling points to the asymptotic variance. Sampling points uniformly distributed over a rectangular make the least squares estimation most efficient.
Non-uniform distributions over a region that deviates from a rectangularworsen the efficiency. Thirdly, the asymptotic variance corresponds withthe heteroskedasticity- robust one in the traditional regression analysis. In
other words, the error term in (3) is regarded as the one with heteroskedastic variance on the frequency domain. Finally, |A|, the area of the samplingregion, is the normalizing factor in the asymptotic distribution. This is thenotifying feature in the asymptotics that contrasts with the sample size inthe traditional central limit theorems for iid sequences.

3.2

consistent asymptotic variance estimation

Let us consider the consistent estimation of the asymptotic variance matrix
in Theorem 2, which shall be used to conduct statistical inference for theparameters in the kernels ψa (u, θ), a = 1, . . . , p. First, we estimate thedensity g(x) in condition C2 that appears in the asymptotic variance matrix.
Let K be a kernel function that is positive and continuous on [0, 1]2 , and−1
2
let Kh (s) = h−1
1 h2 K(s1 /h, s2 /h), s = (s1 , s2 ) ∈ [0, 1] for h = (h1 , h2 ), a
bandwdith. Estimate g(x), x ∈ [0, 1]2 bypĝ(x) =TXX
1ĝat (x),(p + 1)T
a=0 t=1ĝat (x) =

1
nat

nat
X

Kh (sati /A − x),

i=1

11

where sati /A = (sati,1 /A1 , sati,2 /A2 ). Then for xi ∈ [0, 1]2 , i = 1, . . . , N 2 ,mesh points equally dividing [0, 1]2 into N 2 squares, let
21N2b̂g = NX

ĝ(xi )4

i=1
2

1

N
X

 N2

i=1

ĝ(xi )22 .

Then the asymptotic variance matrix evaluated in Theorem 2 as
|A|−1 (2π)2 bg Ω−1 ΣΩ−1 ,is consistently estimated byb̂g Ω̂−1 (θ̂)Σ̂(θ̂)Ω̂(θ̂)−1 ,whereΩ̂(θ) =XΨ̇(ωf ; θ)

( T
X

ωf ∈D

Σ̂(θ) =x̂t (ωf )x̂t (ωf ))Ψ̇(ωf ; θ)0 ,

(11)

t=1


X

0

Ψ̇(ωf ; θ) 

ωf ∈D

T
X


{Lt1 t2 (ωf ; θ) + Mt1 t2 (ωf , θ)} Ψ̇(ωf ; θ)0 ,t1 ,t2 =1

Lt1 t2 (ωf ; θ) = Re{ε̂t1 (ωf ; θ)}Re{ε̂t2 (ωf ; θ)}Re{x̂t1 (ωf )}Re{x̂t2 (ωf )0 },Mt1 t2 (ωf ; θ) = Im{ε̂t1 (ωf ; θ)}Im{ε̂t2 (ωf ; θ)}Im{x̂t1 (ωf )}Im{x̂t2 (ωf )0 },ε̂t (ωf ; θ) = ŷt (ωf ) −pXψ̃a (ωf ; θ)x̂ta (ωf ).

a=1

The consistent estimator for the asymptotic variance makes it possible toconduct the t test for the least squares estimators. The term of (2π)2 |A|−1
in the asymptotic variance cancels in the consistent estimation with the
Riemannian approximation for the integral.

4

Simulation studies

This section examines empirical performances of the least squares estimationon the frequency domain proposed in Section 2.2. We validated the asymptotic consistency and normality under Assumption C2-C4 in the asymptotic
scheme in C1, where C3 assumes stationary only across space, does not
across time. Namely, our estimation is validated even for temporally unit
root processes. We check the empirical performances under C3 when temporal stationary is not necessarily satisfied.
12

First, we study the spatial extension of spurious regression problem foundby Granger and Newbold (1974). They found that, for yt and xt , which aremutually independent unit root processes, the t test for the simple regression
analysis rejects the independence more often than significance level, which isthe so called spurious regression. Now we shall check it for our least squares
estimation.
We simulate yt (s), xt (s), which are mutually independent spatio-temporal
data on A = [0, 20]2 by spatial moving averages of unit root processes. Let
uk , vk , k = 1, . . . , 2000 be knots uniformly distributed on A = [0, 40]2 andsimulate byyt (s) =

2000
X

(1)

ηtk e−||s−uk || , xt (s) =

k=1

2000
X

(2)

ηtk e−||s−vk || ,

(12)

k=1

(i)

where ηt,k , i = 1, 2 are mutually independent unit root processes for each k,simulated by
(i)

(i)

(i)

ηtk = ηt−1,k + εtk , t = 1, . . . , T,(i)where εtk are independent iid standard normal variables. We simulated
xt (s) and yt (s) on two different sets of uniformly distributed 2000 points
on A for t = 1, . . . , T . They are mutually independent processes by thesimulating models.
For the simulated dataset of xt (s) and yt (s) that are mutually independent, we fitted the convolutional regression model in (1) with the CAR(1)kernel in (8). The true kernel function in this case is attained by (8) with
β = 0 for any values of δ > 0, namely δ is not an identifiable parameter in
the independent case. Hence we fix δ with 0.5 and 2.0 in the kernel as examples, and estimate β by the least squares in (5) for the fixed δ = 0.5, 2.0,
where D is the half circle with the centre of origin on the upper half plane
adjusted to make the number of elements to 2000/2 = 1000. Figure 1 depicts the histograms of the t values of β for T = 30 and 200, evaluated by
100 simulations.
Figure 1 demonstrates that the asymptotic theory in section 3.1 works fornon-stationary unit root processes across time. Namely, the t values follow
distributions well approximated by standard normal one. The sizes of thettest are 3, 6, 6 and 6 for the cases of δ = 0.5, T = 30. δ = 0.5, T = 200,
δ = 2, T = 30 and δ = 2, T = 200, under the 5% significance level. Hence it isclarified empirically as well as asymptotically that convolutional regressions
do not suffer from the spurious regression problems indicated by Granger
and Newbold (1974).
Next, let us move on to dependent cases when yt (s) is given by a convo-

13

0.6

delta=0.5,T=200

−2

0

2

0.4
0.2

4

−4

−2

0
t−valuedelta=2,T=30

delta=2,T=200
0.6t−value

4

2

4

0.4

mean=−0.029, std=1.05

0.0

0.0

0.2

Density0.4

mean=−0.14, std=1.00

2

0.2

0.6

−4

Densitymean=−0.049, std=0.94

0.0

0.2

Density0.4

mean=−0.15, std=0.90

0.0

Density

0.6

delta=0.5,T=30

−4

−2

0

2

4

−4

−2

t−value0t−value

Figure 1: The histograms of t values for the mutually independent spatiotemporal dataset simulated by (12) evaluated by 100 simulations.
lution of xt (s). We simulated yt (s) by, for A = [0, 20]2 ,Zyt (s) = βφ0 (s − u; δ)xt (u)du + εt (u), t = 1, . . . , T,A

where φ0 is the CAR(1) kernel in (8), andxt (s) =

6000
X

(1)
ξtk e−||s−uk || , εt (s) =

k=1
(i)
(i)
(i)
ξtk = ρξt−1,k + ftk , i = 1, 2,

6000
X

(2)

ξtk e−||s−vk || ,

k=1

(i)

and where uk , vk are knots uniformly distributed on A and ftk , i = 1, 2,
are independent standard normal variables. The convolution operator wasapproximated with the Riemannian summation when we simulated in practice. Notice that the simulated dataset satisfy Assumption C4 for any valueof ρ. We have an interest in the estimation performances of β and δ when
ρ is designed as 0 and 1 that simulate stationary and unit root across time,
respectively.
14

mean
RMSE
ave(se)

ρ=0
β
δ
1.01
0.49
0.041 0.023
0.061 0.052

ρ=1βδ
1.04 0.51
0.18 0.10
0.16 0.11

Table 1: The means and RMSEs of LSE minimizing (5) for β = 1, δ = 0.5 in
the CAR(1) kernel in the cases of ρ = 0 and 1, which simulate stationary andunit root across time, respectively, evaluated by 100 simulations. ave(se) isthe average of estimated standard errors in (11).
For 100 sets of yt (s), xt (s) simulated for β = 1, δ = 0.5 on 6000 uniformly
scattered points on A with T = 30 in cases 1 and 2 of ρ = 0 and 1, respectively. We estimated β, δ in the CAR(1) kernel by minimizing (5). Table 1
lists the average and root mean squared error of β̂ and δ̂.
Table 1 demonstrates that the least squares by minimizing (5) works
efficiently for both cases of ρ = 0, error term is stationary across time, andρ = 0, unit root across time. They have small biases and root mean squaresthat are well approximated by the standard errors in (11) based on theasymptotic variance in Theorem 2. The simulation results validate that theasymptotic distribution approximates the empirical distribution accurately
for both cases of stationary and unit root non-stationary of error terms.

5Real data analysis

In this section, we demonstrate the applications of the convolutional regression model in (1) to two real spatio-temporal data. They are monthly
precipitation recorded at US weather stations and weekly number of newcases of COVID-19 in every city all over Japan. Observation points in both
examples are irregularly spaced in the way we supposed for the estimationin section 2.2.
Here let us denote several remarks for the demonstrations. First, weshall employ CAR(1) kernels in (8) for all the convolutional regression models in the section. Secondly, all the spatial data we considered for dependentand independent variables at all temporal points were detrended in order
to avoid the intercept issues in Section 2.1 as well as as well as to satisfy
stationary conditions in Assumption C4 in Section 3. Specifically we detrended separately for each temporal points by fitting CAR(1) random fieldby Brockwell and Matsuda (2017),ψj (s; α) = e−α||s−uj || , α > 0, s ∈ R2 ,(13)

where uj , j = 1, . . . , N , are knots scattered randomly over the observation
region. Finally, we took two benchmarks for comparisons with the forecast15

Figure 2: Weather stations in United States
ing performances. For our forecast with the least squares estimation on thefrequency domain in (5), the first benchmark is the predicted value with theestimation on the spatial domain in (6), and the other one is the trend that
we identified for the removal of the trend function by CAR(1) random field.

5.1

US weather data

Here we apply the functional regression model in (1) to US precipitation
data, which is monthly total precipitation, maximum and minimum temperatures observed at weather stations all over the US from 1895 through 1997,
available on the Institute for Mathematics Applied Geosciences (IMAG)
website. All the observations are recorded monthly as spatio-temporal datawith their longitudes and latitudes. For more details, visit
http://www.image.ucar.edu/Data/US.monthly.met/USmonthlyMet.shtml.
Around 6,000 weather stations are shown in Figure 2.
Regarding the monthly precipitation, and maximum and minimum temperatures as dependent and independent variables, respectively, we fitted thefunctional spatial regression model in (1) given by, for the CAR(1) kernel in(10),
Z
pptt (s) = β1 φ0 (s − u; δ1 )tmint (u)duZ+ β2 φ0 (s − u; δ2 )tmaxt (u)du + εt (s).
(14)
We took the observations from 1986 to 1990 as in-samples and those in16

variable
min temp.
max temp.

param.
β1
δ1
β2
δ2

LSE by (5)
1.97
0.62
−1.64
0.38

se
(0.15)(0.053)(0.10)(0.023)

Table 2: The estimated parameters for (14) with their standard errors,
when we took A = [0, 50] × [0, 30] to construct the periodogram on the meshpoints (2πi/A1 , 2πj/A2 ), i, j 6= 0 in D = {ω ∈ R2 , ||ω|| < K}, where K was
adjusted to make the number of the elements be 4000.

model
CAR(1) by (5)CAR(1) by (6)trend fit by (13)forecasting MSE
in-sample out-of-sample
23.7
28.8
23.4
28.6
28.7
36.3

Table 3: Forecasting MSE for in-samples (60 months in 1986-1990) andout-of-samples (12 months in 1991) in US precipitation over 6000 weather
stations.
1991 as out-of-samples. The precipitation, maximum, and minimum temperature do not share the same observation points that differ from month tomonth. Thus the numbers of the observations depend on month. They arearound 6000 for precipitation and 4000 for maximum and minimum temperature on average. The total period is 72 months, 60 months for in-samples
and 12 months for out-of-samples.
We transformed the coordinate of latitudes and longitudes to record theobservation locations into Cartesian coordinate based on global distancewith an unit of 100 km. Applying the discrete Fourier transform to thein-sample observations, we estimated the parameters, β1 , β2 , δ1 , δ2 by minimizing the equation in (5) with (9), where we took A = [0, 50] × [0, 30] toconstruct the periodogram on the mesh points (2πi/A1 , 2πj/A2 ), i, j 6= 0 in
D = {ω ∈ R2 , ||ω|| < K}, where K was adjusted to make the number of theelements be 4000. The estimation results are shown in Table 2.
Then we constructed the forecasts by (14) in (7) with the estimatedparameters to compare the mean squared errors with those of the benchmarks. The MSEs evaluated by, for vkt , k = 1, . . . , nt , the out-of-samples at
t = 61, . . . , 72,
72

n

t
XX

1
ˆ t (vkt ) 2
pptt (vkt ) − ppt
n61 + · · · + n72

t=61 k=1

are shown in Table 3.
17

We find from Table 2 that we identified the non-trivial relations. Namely,
precipitation is negatively correlated with maximum temperatures, while itis positively correlated with minimum temperatures. The evaluated standard errors for the estimators indicate that they are strongly significant
with the p-values less than 1%. Table 3 demonstrates that our estimationon the frequency domain works as well as that on the spatial domain, while
it outperforms the benchmark of the fitted trend. The two least squares
estimation on the frequency and spatial domains shares almost equivalent
results, although the one on the spatial domain is slightly better. It may
come from the approximation to evaluate the out-of-sample forecats in (7)that gives advantage to the estimation on the spatial domain in (6) over theone on the frequency domain in (5).

5.2

COVID-19 outbreak in Japan

This section applies the convolutional regression in (1) with CAR(1) kernels
to analysis of COVID-19 outbreak in Japan in relation with human mobilities. We aim to identify the impact of human mobility to outbreak of
COVID-19.
We collected the weekly new cases of COVID-19 infections at all of 1,896
cities in Japan as dependent variables, while we evaluated weekly humanmobility based on dataset provided by NTT DOCOMO, Inc., one of thelargest Japanese mobile phone companies as an independent variable. Moreprecisely, the dependent variable is the covid19 infection rate, which is, fornti and Nti , weekly new cases and population at a city i and a week t, definedby,nti
covid19t (si ) =× 10, 000.
Nti
NTT Docomo, inc. provides spatio-temporal dataset that counts numbersof people in 500 meter meshes at every one hour period all over Japan. We
aggregate it into city level that counts the weekly number of people whoenter a city from outside the city. For mti , weekly number of people whoenter a city i at si from outside at a week t, and define the independent
variable bymti
hmobilt (si ) =.
Nti
In addition, we collected weekly total precipitation at 936 weather stationsin Japan as one more possible independent variable.
All the variables have longitudes and latitudes to indicate their locations.
We transformed the coordinate of latitudes and longitudes into Cartesian
coordinate based on global distance with an unit of 50 km. We took thedataset collected weekly for 86 weeks from Feb 12th, 2020 to Sep. 30th,
2021, and the last ten weeks were separated as the out-of -samples.
18

Figure 3: Locations of 1,896 Cities in Japan
We fitted the convolutional regression that includes lagged terms of dependent variables ascovid19t (s)ZZ= β1 φ0 (s − u; δ1 )covid19t−1 (u)du + β2 φ0 (s − u; δ2 )covid19t−2 (u)duZZ+ β3 φ0 (s − u; δ3 )hmobilt−1 (u)du + β4 φ0 (s − u; δ4 )hmobilt−2 (u)duZ+ β5 φ0 (s − u; δ5 )pptt−2 (u)du + εt (s),(15)
for the CAR(1) kernel φ0 (·; δ) in (10).
We estimated the parameters by minimizing (5) with (9), where wetook A = [0, 40] × [0, 40] to construct the periodogram on the mesh points(2πi/A1 , 2πj/A2 ), i, j 6= 0 in D = {ω ∈ R2 , ||ω|| < K}, where K was adjusted to make the number of the elements be 1800. The estimation results
are shown in Table 4.
Then we constructed the forecasts by (15) with (7) with the estimatedparameters to compare the mean squared errors with those of the benchmarks. The MSEs evaluated by, for vkt , k = 1, . . . , nt , t = 77, . . . , 86, theout-of-sample period,
86

n

t 
2
XX
1
ˆ t (vkt )covid19t (vkt ) − covid19
n77 + · · · + n86

t=77 k=1

are shown in Table 5.
We find from Table 4 that covid-19 is a strongly autocorrelated data by
the estimators for the first and second lagged autoregressive terms. Hmobil
19variable
covid−1
covid−2
hmobil−1
hmobil−2
ppt−2

param.
β1
δ1
β2
δ2
β3
δ3
β4
δ4
β5
δ5

LSE by (5)
0.99
0.14
−0.39
0.45
−0.50
0.11
0.59
0.12
−0.00310.13

se(0.021)(0.0078)(0.027)(0.027)(0.062)(0.052)(0.065)(0.043)(0.00092)(0.19)

Table 4: The estimated parameters for (15) with their standard errors,
when we took A = [0, 40] × [0, 40] to construct the periodogram on the meshpoints (2πi/A1 , 2πj/A2 ), i, j 6= 0 in D = {ω ∈ R2 , ||ω|| < K}, where K was
adjusted to make the number of the elements be 1800.
is negatively correlated with covid-19 at the first lag, while it is positively
related at the second lag, when we control the first and second lagged covid19. The negative correlation with the first lag means that human mobility isstrictly controlled by local or central government. Namely the government
predicts the tendency of increasing or decreasing new cases and fixes thepolicies to control the mobilities based on the prediction, which can produce
the negative correlation at the first lag. On the other hand, the second lagged
human mobility positively related with covid-19, which demonstrates thewell known result that it takes around two weeks for a man to be infectedwith the virus. Finally we find that the precipitation at the second lag
correlates covid-19 negatively. From the fitted models that are not shown
here, it is seen that the precipitation at any other lags is not significant.
Although it is not clear which causes the negative correlation in two weeks,
restricted human mobility by rain or effects of humidity, it is interesting todetect the negative relationship by the convolutional regression.
Table 5 indicates that the least square estimation on frequency domain
has larger MSEs than that on the spatial domain in terms of both in-sample
and out-of-sample comparisons, although they should be theoretically equivalent. The difference here is even greater than that in the US weather datacase in Table mse. The greater difference of MSE as well as the estimation
between the two estimation may come from the difference of the sample
sizes. There are around 6,000 weather stations in the first example, while
there are around 1,800 cities recorded for weekly new cases of infections.
Frequency domain estimation requires more sample sizes than that on thespatial domain to attain the equivalent approximation accuracy.

20

modelCAR(1) by (5)CAR(1) by (6)trend fit by (13)forecasting MSE
in-sample out-of-sample
1.82
12.7
1.71
11.8
2.45
31.3

Table 5: Forecasting MSE for in-samples (77 weeks in Feb. 2020-Jul. 2021)
and out-of-samples (8 weeks in Aug. -Sep. 2021) in Covid-19 new cases over
1,760 cities in Japan.

6

Discussion

This paper tries a parametric convolutional regression for big spatial dataon huge scales at discrete time points. Applying Fourier transform to theregression, we propose the least squares estimation on the frequency domain.
We have found several features of the convolutional regression for big spatial
data.
First, our approach can handle a regression among big spatial data on
national or continental scale. Second, our approach does not require stationarity across time but do across space, while existing approaches on functional
PCA need not stationarity across space but do across time. Non-stationary
behaviours including unit root across time or even one time point observation
suffices to estimate parametric convolution kernels consistently. Stationarity
across space necessary for our approach to work is the price for no restrictions across time. As a result, our regression is free from spurious regression
problem indicated by Granger and Newbold (1974), which is well demonstrated in the simulation studies. Thirdly, our approach allows irregularly
spaced data over space and irregularity can change with time varying sample sizes across time, which is common in collecting big spatial data. The
discrete Fourier transform works critically to make it possible to estimate
properly under irregular sampling across space. Finally, the leat squares
estimation on the frequency domain is validated asymptotically under theso called mixed asymptotics. Since the consistent asymptotic variance matrix estimator is easily constructed, statistical inference for convolutional
regression becomes possible. The t test in the empirical study of Covid-19
analysis in relation with human mobility data collected on national scalein Japan detects the significant relation statistically that two week lagged
human mobility associates positively current week Covid-19 new cases.

21

7

ProofsWe define L(θ) by normalizing Q(θ) in (5) as
L(θ) = (2π)2 Q(θ)To(2π)2 X X n=
It,yy (ωf ) − 2φ̃(ωf ; θ)It,xy (ωf ) + φ̃(ωf ; θ)It,xx (ωf )φ̃0 (ωf ; θ) ,|A|
ωf ∈D t=1

where φ̃(ω; θ) = (φ̃1 (ω; θ), . . . , φ̃p (ω, θ)) and0It,xx (ω) = |A|x̂t (ω)x̂t (ω) ,It,xy (ω) = |A|x̂t (ω)ŷt (ω),
It,yy (ω) = |A|ŷt (ω)ŷt (ω),
the periodograms.

7.1

Proof of Theorem 1

Let θ1 ∈ Θ be a parameter that is not equal to θ0 . By applying Lemmas 1
and 2 to the periodograms with the relations,ft,yy (ω) = φ̃(ω; θ0 )ft,xx (ω)φ̃0 (ω; θ0 ) + ft,εε (ω), t = 1, . . . , T,
ft,xy (ω) = ft,xx (ω)φ̃0 (ω; θ0 ),ft,xε (ω) = 0,
we haveL(θ1 ) →τg

T Z 
X
t=1

φ̃(ω; θ1 ) − φ̃(ω; θ0 ) ft,xx (ω) φ̃(ω; θ1 ) − φ̃(ω; θ0 )

0+ ft,εε (ω) dωD:= L∞ (θ1 ),say, in probability as k tends to be infinity in the asymptotic regime of C1.
Since ft,xx (ω) is positive definite,
L∞ (θ1 ) − L∞ (θ0 ) =
T Z 


0
X
φ̃(ω; θ1 ) − φ̃(ω; θ0 ) ft,xx (ω) φ̃(ω; θ1 ) − φ̃(ω; θ0 ) dω
τg
t=1

D

> 0.
It follows that, for any positive constant K(θ0 , θ1 ) that is smaller thanL∞ (θ1 ) − L∞ (θ0 ),lim P (L(θ0 ) − L(θ1 ) < −K(θ0 , θ1 )) = 1.

k→∞

22

Let Mt (ω) be the maximum eigenvalue of the periodogram matrix!
0It,yy (ω) It,xy (ω)
It (ω) =, ω ∈ D, t = 1, . . . , T.
It,xy (ω) It,xx (ω)
Then for any θ1 and θ2 that satisfy |θ2 − θ1 | < δ,|L(θ2 ) − L(θ1 )| < Cδ sup max Mt (ω) = Hδ ,ω∈Dtsay. It is seen that there exists a δ > 0 such thatlim P (Hδ < K(θ0 , θ1 )) = 1.

k→∞

Applying Lemma 2 in Walker (1964), we have the consistency.

7.2

Proof of Theorem 2

θ̂)
Applying Taylor series expansion to ∂L(∂θ at θ0 , we have0=

∂L(θ0 ) ∂ 2 L(θ∗ )∂L(θ̂)=+(θ̂ − θ0 ),
∂θ
∂θ
∂θ∂θ0

where θ∗ is the mean value between θ0 and θ̂. Hencep |A| θ̂ − θ0 =

 2
−1


p
∂ L(θ∗ )∂L(θ0 )
|A| −
.
∂θ∂θ0
∂θ

The first factor, which is the Hessian matrix, is evaluated as,()(2π)2 X ∂ φ̃(ωk ; θ∗ ) X∂ φ̃0 (ωk ; θ∗ )
2
I
(ω
)
t,xx
k
|A|
∂θ0
∂θ
t
ωk ∈D

p

o ∂ 2 φ̃ (ω ; θ∗ )(2π)2 X X X na fIt,yxa (ωf ) − It,xa x (ωf )φ̃(ωf ; θ∗ )
,
−2
|A| t
∂θ∂θ0
ωf ∈D a=1

which converges in probability by Lemmas 1 and 2 to2τg (2π)2 Ω,(16)

since θ∗ converges to θ0 under the consistency, and ft,yx (ω) = ft,xx (ω)φ̃0 (ω; θ0 ).
The second factor, which is the score vector, is evaluated aspTno|A|(2π)2 X X ∂ φ̃(ωf ; θ0 )−2
|A|x̂t (ωf ) ŷt (ωf ) − φ̃(ωf ; θ0 )x̂t (ωf )|A|∂θ
ωf ∈D t=1

23

which is, by Lemma 2, equal to( T)p|A|(2π)2 X ∂ φ̃(ωf ; θ0 ) X−2It,xε (ωf ) + J,
|A|
∂θ
ωf ∈D

(17)

t=1

where
It,xε (ω) = |A|x̂t (ω)ε̂t (ω),
n0t1 X
εt (s0tj ) exp(−iω 0 s0tj ),ε̂t (ω) =n0t
j=1

and E|J| is bounded by
C max max max
ω∈D 1≤t≤T 1≤a≤p
v
u
2
u
n0t
nat
1/2 X
u 1/2 |A|1/2 X
0 ω
|A|
0
−is
×t|A| Exta (s0tj − u)e−i(s0tj −u) ω −xta (satj )e atj
n0t
nat
j=1

j=1

=o(1),by Lemma 1. Hence it is equivalent to show the asymptotic distribution of)( T
p Z ∂ φ̃(ω; θ0 ) X
K = 2 |A|It,xε (ω) dω,
∂θ
D
t=1

by applying Lemma 3 to the first term.
We shall show in Lemmas 4, 5 and 6 that
E(Ka ) → 0, a = 1, . . . , q,E(Ka K b ) → 4(2π)6 τg2 Σab , a, b = 1, . . . , q,cum (Ka1 , . . . , Kar )) → 0, for 1 ≤ a1 , . . . , ar ≤ q, r ≥ 3,
respectively, as k tends to ∞, which proves the asymptotic normality in
Theorem 2 in combinations with (16).

8LemmasFor the unified treatment of the periodograms, we define, for t = 1, . . . , T ,zt0 (s) = yt (s),zta (s) = xta (s), a = 1, . . . , p,and, re-define the periodogram and spectral density matrix of z by, fora, b = 0, 1, . . . , p,It,ab (ω) = |A|ẑta (ω)ẑtb (ω),
ft,ab (ω) = ft,z,ab (ω), a, b = 0, 1, . . . , p,
24

Lemma 1.


2

EIt,ab (ω) = (2π) τg ft,ab (ω) + O

ntab |A|−2−2+ A1 + A2, a, b = 0, 1, . . . , p,
nta ntb

where nta , ntb and ntab are the number of elements in Sta , Stb and Sta ∩ Stb ,respectively, andZ|g(x)|2 dx.
τg =[0,1]2

Proof. This is an immediate consequence of Lemma 3 in Matsuda and Yajima (2009).
Lemma 2. For a square integrable function ψ(ω), ω ∈ R2 ,
 (2π)2 X
varIt,ab (ωj )ψ(ωj ) = O |A|−1 , a, b = 0, 1, . . . , p.
 |A|

ωj ∈D

Proof. Letψ̂(s) =

0(2π)2 Xψ(ωj )e−iωj s , s ∈ A = [0, A1 ] × [0, A2 ],
|A|

ωj ∈D

which is extended periodically to [−A1 , A1 ] × [−A2 , A2 ]. Then the object
for the variance is evaluated as
n

n

ta X
tb
|A|3/2 X
zta (satc )zb (sbtd )ψ̂(satc − sbtd ).
nta ntb

c=1 d=1

The variance is given by|A|3 X X X X
E 2 2
cum (zta (satc1 ), ztb (sbtd1 ), zta (satc2 ), ztb (sbtd2 ))nta ntb c
c
1

d1

2

d2


+ γt,aa (satc1 − satc2 )γt,bb (sbtd1 − sbtd2 ) + γt,ab (satc1 − sbtd2 )γt,ba (sbtd1 − satc2 )× ψ̂(satc1 − sbtd1 )ψ̂(satc2 − sbtd2 )Z Z Z Z =cum(zta (u1 ), ztb (v1 ), zta (u2 ), ztb (v2 ))
A

A

A

A


+ γt,aa (u1 − u2 )γt,bb (v1 − v2 ) + γt,ab (u1 − v2 )γt,ba (v1 − u2 )× ψ̂(u1 − v1 )ψ̂(u2 − v2 )|A|−1 g(u1 /A)g(v1 /A)g(u2 /A)g(v2 /A)du1 dv1 du2 dv2 + o(1).

25The first term is, by expressing the cumulant term with the cumulant spectrum:mXft,abab (ω1 , ω2 , ω3 ) =κt,ef gh G̃tae (ω1 )G̃tbf (ω2 )G̃tag (ω3 )G̃tbh (ω1 + ω2 + ω3 ),e,f,g,h=1(18)
given by
Z Z Z

0

0

0

ft,abab (ω1 , ω2 , ω3 )eiω1 (u1 −v2 ) eiω2 (v1 −v2 ) eiω3 (u2 −v2 ) dω1 dω2 dω322
R2 R
Z ZR Z Z|A|−1ψ̂(u1 − v1 )ψ̂(u2 − v2 )g(u1 /A)g(v1 /A)g(u2 )/Ag(v2 /A)du1 dv1 du2 dv2 ,
A

A

A

A

which is, by Schwarz inequality, bounded bysZ Z Z
2
Y|ft,abab (ω1 , ω2 , ω3 )|Pj dω1 dω2 dω3 ,
R2

j=1

R2

(19)

R2

for
−1

Z Z

0

A

P2 = |A|−1

A

Z Z
A

2

0

ψ̂(u1 − v1 )g(u1 /A)g(v1 /A)eiω1 u1 eiω2 v1 du1 dv1 ,P1 = |A|002ψ̂(u2 − v2 )g(u2 /A)g(v2 /A)e−i(ω1 +ω2 ) v2 eiω3 (u2 −v2 ) du2 dv2 .

A

By applying Perseval’s equality to both terms, (19) is bounded by
Z A1 Z A2C2

ψ̂(s) ds,
−A1

(20)

−A2

which is evaluated as
4C|A|(2π)4 X|ψ(ωj )|2 < C 0
|A|2
ωj ∈D

Z

|ψ(ω)|2 dω = O(1).

D

Also the second and third terms in the variance are bounded by a constantwith the same argument, which completes the proof.
Lemma 3. For a square integrable function ψ(ω), ω ∈ R2 ,pp Z
|A|(2π)2 X
It,xa ε (ωj )ψ(ωj ) − |A|
It,xa ε (ω)ψ(ω)dω
|A|
D
j∈JD

= op (1), a = 1, . . . , p, t = 1, . . . , T.

26

Proof. Let0(2π)2 Xψ(ωj )e−iωj s , s ∈ A = [0, A1 ] × [0, A2 ],
|A|
ωj ∈D
Z
0
ψ̃(s) =
ψ(ω)e−iω s dω, s ∈ R2 ,ψ̂(s) =

D

the first one of which is extended periodically to [−A1 , A1 ] × [−A2 , A2 ]. For
δ(s) = ψ̂(s) − ψ̃(s), the difference is evaluated as
n

n

L=

ta X
t0
|A|3/2 X
xa (stac )ε(st0d )δ(stac − st0d ).
nta nt0

c=1 d=1

E|L|2 , which is similarly evaluated till (20) in Lemma 2, is bounded by
Z A1 Z A2(21)|δ(s)|2 ds.
C
−A1

−A2

Notice that ψ̂(s), ψ̃(s) are square integrable, since
Z A1 Z A2
Z
X
2
4
−2
2
|ψ̂(s)| ds = |A|(2π) |A||ψ(ωj )| < C|ψ(ω)|2 dω < ∞,
0

0

D

ωj ∈D

Z

|ψ̃(s)|2 ds = (2π)2R2Z|ψ(ω)|2 dω < ∞.

D

It follows that, for any ε > 0, there exists a compact set BM = [−M1 , M1 ] ×[−M2 , M2 ] ⊂ [−A1 , A1 ] × [−A2 , A2 ] such that
Z A1 Z A22
ψ̂(s) − ψ̂(s)IBM (s) ds < ε,−A1 −A2Z
2ψ̃(s) − ψ̃(s)IBM (s) ds < ε.
R2Then (21) is bounded by
Z A 1 Z A 2
Z22
Cψ̂(s) − ψ̂(s)IBM (s) ds +ψ̂(s) − ψ̃(s) ds
−A
−A2
BM

Z 1
2
+
ψ̃(s)IBM (s) − ψ̃(s) ds
R20< C ε,
which completes the proof.
Lemma 4.
E(Ka ) = O−3/2 1/2
1/2 −3/2
A1 A2 + A1 A2+

T
X|A|3/2
t=1

27

nta

!
, a = 1, . . . , q.

Proof. It follows immediately from Lemma 3 in Matsuda and Yajima (2009)
and ft,εx (ω) = 0.
Lemma 5.
E(Ka K b ) → 4(2π)6 τg2 Σab , a, b = 1, . . . , q,
Proof. Re-express Ka by
p
n0t
nct X
T X
X
|A|3/2 X
t=1 c=1n0t nct

Xtc (sctj )εt (s0tl )ψ̃ac (sctj − s0tl ),(22)

j=1 l=1

where∂ φ̃c (ω; θ0 )
,
∂θa
Z
0
ψ̃ac (s) =ψac (ω)e−iω s dω,ψac (ω) =D

First, ψ̃ac (s) in Ka defined in (22) may be replaced withMψ̃ac(s) = ψ̃ac (s)IBM (s),for a sufficiently large compact set BM = [−M1 , M1 ] × [−M2 , M2 ], since thevariance of the difference between them may be made arbitrary small byfollowing the argument till (20) in Lemma 2. Then EKa K̄b replaced withM (s) and ψ̃ M (s) is evaluated as
ψ̃ac
bc
T
X

p
X

−1

Z Z Z Z|A|A

t1 ,t2 =1 c1 ,c2 =1

A

A

A

M
M
γx,t1 t2 ,c1 c2 (u1 − u2 )γε,t1 t2 (v1 − v2 )ψ̃ac(u1 − v1 )ψ̃bc(u2 − v2 )1
2×g(u1 /A)g(v1 /A)g(u2 /A)g(v2/A)du1 dv1 du2 dv2 + o(1).
The first term is evaluated as
Z Z
p
T
X
X
fx,t1 t2 ,,c1 c2 (ω1 )fε,t1 t2 (ω2 )|A|−1 P dω1 dω2
2
t1 ,t2 =1 c1 ,c2 =1 R

R2for
Z Z Z Z
P =
A

A

A

A

0

0

M
M
ψ̃ac1
(u1 − v1 )ψ̃bc(u2 − v2 )eiω1 (u1 −u2 ) eiω2 (v1 −v2 )2×g(u1 /A)g(v1 /A)g(u2 /A)g(v2 /A)du1 dv1 du2 dv2 .
By change of variables by u1 − v1 = l1 , u2 − v2 = l2 and the compactness ofthe supports of φ̃M , |A|−1 P is evaluated asMM(2π)4 ψac(ω1 )ψbc(ω1 )Fg2 (ω1 + ω2 ) + o(1),1228whereFg2 (ω) = |A|−1

Z

M
ψac
(ω) = (2π)−2120g 2 (u/A)eiω u du ,

ZA

0

BM

M
ψ̃ac
(s)eiω s ds.
1

It follows by Lemma 1(c) in Matsuda and Yajima (2009) that the first term
converges to
p
X

T
X

6

(2π) τg2

t1 ,t2 =1 c1 ,c2 =1

Z
R2

M
M
(ω)dω.
(ω)ψbcfx.t1 t2 ,c1 c2 (ω)fε,t1 t2 (ω)ψac
1
2

We have the result arbitrary close to (2π)6 τg2 Σab by taking BM large.
Lemma 6. For r ≥ 3,
cum(Ka1 , . . . , Kar ) = O |A|−r/2+1 .
Proof. First, ψ̃ac (s) in Ka defined in (22) may be replaced withMψ̃ac(s) = ψ̃ac (s)IBM (s), a = 1, . . . , q,

for a sufficiently large compact set BM = [−M1 , M1 ] × [−M2 , M2 ], since thevariance of the difference between them may be made arbitrary small byfollowing the argument till (20) in Lemma 2.
M (s) is evaluated as
The cumulant for the one replaced with ψ̃accum(Ka1 , . . . , Kar ) =

q
X

T
X

···

t1 ,...,tr =1 c1 ,d1 =1

Z|A|−r/2

cr ,dr =1

Z

×

···
A
r
Y

×

q
X

cum(xt1 ,c1 (u1 )εt1 (v1 ), . . . , xtr ,cr (ur )εtr (vr ))Aψ̃aMj cj (uj − vj )g(uj /A)g(vj /A)duj dvj + o(1),j=1in which the cumulant is evaluated ascum(xt1 ,c1 (u1 ), . . . , xtr ,cr (ur )) × cum(εt1 (v1 ), . . . , εtr (vr ))ZZ
r−1
Y iω0 (u −u )=···fx,t1 ...tr ,c1 ···cr (ω1 , . . . , ωr−1 )e j j r dω1 · · · dωr−1
R2

R2

Z
×

···
R2

j=1

Z
R2

fε,t1 ...tr (λ1 , . . . , λr−1 )

r−1
Y
j=1

29

0

eiωj (vj −vr ) dλ1 · · · dλr−1 .

By replacing the cumulant term with the spectrum expression, the summandof the corresponding term is given byZZ−r/2
fx,t1 ...tr ,c1 ···cr (ω1 , . . . , ωr−1 )dω1 · · · ωr−1
···
|A|
2
R2
ZR
Z
fε,t1 ...tr (λ1 , . . . , λr−1 )dλ1 · · · λr−1
···
×
×

R2
r−1
YZ

j=1 A

R2

Z

0

A

0

ψ̃aMj cj (uj − vj )eiωj uj eiλj vj g(uj /A)g(vj /A)duj dvj

Z Z

0

0

ψ̃aMr cr (ur − vr )e−i(ω1 +···+ωr−1 ) ur e−i(λ1 +···+λr−1 ) vr g(ur /A)g(vr /A)dur dvr
A A
Z
Z
−r/2
fx,t1 ...tr ,c1 ···cr (ω1 , . . . , ωr−1 )dω1 · · · ωr−1
···
=C|A|

×

R2

Z
···

×
R2

R2

Z
R2

fε,t1 ...tr (λ1 , . . . , λr−1 )dλ1 · · · λr−1 ×

r−1Y

D2 (ωj + λj )ψaMj cj (ωj )

j=1

×D2 (−ω1 − λ1 − · · · − ωr−1 − λr−1 )ψaMj cj (−ω1 − · · · − ωr−1 ) + o(1),whereZ
D2 (ω) =0g 2 (u/A)e−ω u du.

A

We find it to be O(|A|−r/2+1 ) by Lemma 2 in Matsuda and Yajima (2009).

References[1] Benatia, D., Carrasco, M. and Florens, J. (2017) Functional linear regression with functional response, Journal of Econometrics, 201, 269-291.
[2] Bosq, D. (2000). Linear Processes in Function Spaces. Springer-Verlag,
New York.
[3] Brockwell, P. J. and Davis, R. A. (1991) Time Series: Theory and Methods, 2nd edn. Springer-Verlag.
[4] Brockwell, P. J. and Matsuda, Y. (2017). Continuous auto-regressive
moving average random fields on Rn . Journal of the Royal StatisticalSociety: Series B, 79(3), 833-857.
[5] Crambes, C. and Mas, A. (2013). Asymptotics of prediction in functional
linear regression with functional outputs. Bernoulli 19, 2627–2651.

30

[6] Ferraty, F., Vieu, P. (2006). Nonparametric Functional Data Analysis:
Methods, Theory, Applications and Implementations. Springer-Verlag,
London.
[7] Granger, C. W. J. and Newbold, P. (1974). Spurious regressions in econometrics. Journal of Econometrics, 2, 111-120.
[8] Horvath, L., Kokoszka, P. (2012). Inference for Functional Data with
Applications. Springer, New-York.
[9] Hosoya, Y. (1997). A limit theory for long-range dependence and statistical inference on related models Annals of Statistics, .25, 1, 105-37.
[10] Kokoszka, P. and Reimherr, M. (2017). Introduction to Functional Data
Analysis. CRC Press, Boca Raton.
[11] Liu, X., Xiao, H. and Chen, R. (2016). Convolutional autoregressivemodels for functional time series. Journal of Econometrics, 194, 263282.
[12] Matsuda, Y. and Yajima, Y. (2009). Fourier analysis of irregularly
spaced data on Rd . Journal of the Royal Statistical Society: Series B,
71(1), 191-217.
[13] Matsuda, Y. and Yajima, Y. (2018). Locally stationary spatio-temporalprocesses. Jpn. J. Stat. Data Sci., 1, 41-57.
[14] Ramsay, J.O., Silverman, B.W.(2005). Functional Data Analysis, second ed. Springer, New York.
[15] Robinson, P. M. (1995). Gaussian semiparametric estimation of longrange dependence. Annals of Statistics, 23, 1630-1661.
[16] Stein, M. L. (1999). Interpolation of Spatial Data. Springer, New York.
[17] Uhlenbeck, G. E. and Ornstein, L. S. (1930). On the theory of Brownian
Motion, Physical Review, 36, 823-841.
[18] Walker, A. M. (1964). Asymptotic properties of least-squares estimatesof parameters of the spectrum of a stationary non-deterministic timeseries. Journal of the Australian Mathematical Society, 4, 363-384.
[19] Yao, F., Müller, H.G., Wang, J.L. (2005). Functional linear regression
analysis for longitudinal data. Ann. Statist. 33 (6), 2873–2903.

31