Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
81,779 characters · 6 sections · 66 citation commands
Nonparametric prediction with spatial data
Random models for spatial or spatio-temporal data play an important role in many disciplines of economics, such as environmental, urban, development or agricultural economics as well as economic geography, among others. When data is collected over time such models are termed `noncausal' and have drawn interest in economics, see for instance Breidt2001 among others for some early examples. Other studies may be found in the special volume by Baltagi2007 or Cressie. Classic treatments include the work by Mercer1911 on wheat crop yield data (see also Gao2006) or Batchelor1918 which was employed as an example and analysed in the celebrated paper by Whittle1954. Other illustrations are given in Cressie1999, see also Fernandez-Casal2003. With a view towards applications in environmental and agricultural economics, Mitchell2005 employed a model of the type studied in this paper to analyse the effect of carbon dioxide on crops, whereas Genton2008 examine the yield of barley in UK. The latter manuscripts shed light on how these models can be useful when there is evidence of spatial movement, such as that of pollutants, due to winds or ocean currents.
Doubtless one of the main aims when analysing data is to provide predicted values of realizations of the process. More specifically, assume that we have a realization $\mathcal{X}_{n}=\left\{ x_{t_{i}}\right\} _{i=1}^{n}$ at locations $t_{1},...,t_{n}$ of a process $\left\{ x_{t}\right\} _{t\in \mathcal{D}}$, where $\mathcal{D}$ is a subset of $\mathbb{R}^{d}$. We wish then to predict the value of $x_{t}$ at some unobserved location $t_{0}$, say $x_{t_{0}}$. For instance in a time series context, we wish to predict the value $x_{n+1}$ at the unobserved location (future time) $n+1$ given a stretch of data $x_{1},..,x_{n}$. It is often the case that the predictor of $x_{t_{0}}$ is based on a weighted average of the data $\mathcal{X}_{n}$, that is
where the weights $\beta _{1},...,\beta _{n}$ are chosen to minimize the $ \mathcal{L}_{2}$-risk function
with respect to $b_{1},...,b_{n}$. With spatial data, the solution in $ \left( \ref{pred1}\right) $ is referred as the Kriging predictor, see Stein1999, which is also the best linear predictor for $x_{t_{0}}$. Notice that under Gaussianity or our Condition $C1$ below, the best linear predictor is also the best predictor. It is important to bear in mind that with spatial data prediction is also associated with both interpolation as well as extrapolation.
The optimal weights $\left\{ \beta _{i}\right\} _{i=1}^{n}$ in $\left( \ref {pred1}\right) $ depend on the covariogram (or variogram) structure of $ \left\{ x_{t_{1}},...,x_{t_{n}};x_{t_{0}}\right\} =:\left\{ \mathcal{X} _{n};x_{t_{0}}\right\} $, see among others Stein1999 or Cressie . That is, denoting the covariogram by $Cov\left( x_{t_{i}},x_{t_{j}}\right) =:C\left( t_{i},t_{j}\right) $ and assuming stationarity, so that $C\left( t_{i},t_{j}\right) =:C\left( \left\vert t_{i}-t_{j}\right\vert \right) $, we have that the best linear predictor $\left( \ref{pred1}\right) $ becomes
where
When the data is regularly observed, the unknown covariogram function $ C\left( h\right) $ is replaced by its sample analogue
where $n\left( h\right) =\left\{ \left( t_{i},t_{j}\right) :\left\vert t_{i}-t_{j}\right\vert =h\right\} $ and $\left\vert n\left( h\right) \right\vert $ denotes the cardinality of the set $n\left( h\right) $. When the data is not regularly spaced some modifications of $\widehat{C}\left( h\right) $ have been suggested, see Cressie $\left( 1993,p.70\right) $ for details. One problem with the above estimator $\widehat{C}\left( h\right) $ is that it can only be employed for lags $h$ which are found in the data, and hence the Kriging predictor $\left( \ref{pred2}\right) $ cannot be computed if $\left\vert t_{i}-t_{0}\right\vert \not=h$ for any $h$ such that $n\left( h\right) $ is not an empty set. To avoid this problem a typical solution is to assume some specific parametric function $C\left( h\right) =:C\left( h;\theta \right) $, so that one computes $\left( \ref{pred2} \right) $ with $C\left( h;\widehat{\theta }\right) $ replacing $C\left( h\right) $ therein, where $\widehat{\theta }$ is some estimator of $\theta $.
In this paper, we shall consider the situation when the spatial data is collected regularly, that is on a lattice. This may occur as a consequence of some planned experiment or due to a systematic sampling scheme, or when we can regard the (possibly non-gridded) observations as the result of aggregation over a set of covering regions rather than values at a particular site, see e.g. Conley1999, Conley2007a, Bester2011, Wang2013, Nychka2015, Bester2014. As a result of this ability to map locations to a regular grid, lattice data are frequently studied in the econometrics literature, see e.g. Roknossadati2010, Robinson2011 and Jenish2016. Nonsystematic patterns may occur, although these might arise as a consequence of missing observations, see Jenish2012 for a study that covers irregular spatial data.
However contrary to the solution given in $\left( \ref{pred2}\right) $, our aim is to provide an estimator of $\left( \ref{pred1}\right) $ without assuming any particular parameterization of the dynamic or covariogram structure of the data a priori, for instance without assuming any particular functional form for the covariogram $C\left( h\right) $. The latter might be of interest as we avoid the risk that misspecification might induce on the predictor. In this sense, this paper may be seen as a spatial analog of contributions in a standard time series context such as Bhansali1974 and Hidalgo2002.
The remainder of the paper is organized as follows. In the next section, we describe the multilateral and unilateral representation of the data and their links with a Wold-type decomposition. We also describe the canonical factorization of the spectral density function, which plays an important role in our prediction methodology described in Section (ref), wherein we examine its statistical properties. Section (ref) describes a small Monte-Carlo experiment to gain some information regarding the finite sample properties of the algorithm, and compares our frequency domain predictor to a potential `space-domain' competitor. Because land value and real-estate prices comprise classical applications of spatial methods, see e.g. IversenJr2001, Banerjee2004, Majumdar2006, we apply the procedures to prediction of house prices in Los Angeles in Section (ref). Finally, Section (ref) gives a summary of the paper whereas the proofs are confined to the mathematical appendix.
Before we describe how to predict the value of the process $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{d}}$ at unobserved locations, for $d\geq 1$ , it is worth discussing what do we understand by multilateral and unilateral representations of the process and, more importantly, the link with the Wold-type decomposition. Recall that in the prediction theory of stationary time series, i.e. when $d=1$, the Wold decomposition plays a key role. For that purpose, and using the notation that for any $a\in \mathbb{Z} ^{d}$, $a=\left( a\left[ 1\right] ,...,a\left[ d\right] \right) $, so that $ t-j$ stands for $\left( t\left[ 1\right] -j\left[ 1\right] ,....,t\left[ d \right] -j\left[ d\right] \right) $, we shall assume that the (spatial) process $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{d}}$ admits a representation given by
where the $\varepsilon _{t}$ are independent and identically distributed random variables with zero mean, unit variance and finite fourth moments. The model in $\left( \ref{a1}\right) $ denotes the dynamics of $x_{t}$ and it is known as the multilateral representation of $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{d}}$. It is worth pointing that a consequence of the latter representation is that the sequence $\left\{ \varepsilon _{t}\right\} _{t\in \mathbb{Z}^{d}}$ loses its interpretation as being the \textquotedblleft prediction\textquotedblright\ error of the model, and thus they can no longer be regarded as innovations, as was first noticed by Whittle1954. When $d=1$, this multilateral representation gives rise to so-called noncausal models or, in Whittle1954's terminology, linear transect models. These models can be regarded as forward looking and have gained some consideration in economics, see for instance Lanne2011, Davis2013, Lanne2013 or Cavaliere2018.
It is worth remarking that, contrary to $d=1$, it is not sufficient for the coefficients $\psi _{j}$ in $\left( \ref{a1}\right) $ to be $O\left( \left\vert j\right\vert ^{-3-\eta }\right) $ for any $\eta >0$ as our next example illustrates. Indeed, suppose that $\psi _{j}=\left( j\left[ 1\right] +j\left[ 2\right] \right) ^{-4}=O\left( \left\Vert j\right\Vert ^{-4}\right) $. However it is known that the sequence $\left\{ \sum_{\ell =1}^{d}j^{2} \left[ \ell \right] \right\} \left\vert \psi _{j}\right\vert $ is not summable. That is, see for instance Limaye2009,
One classical parameterization of $\left( \ref{a1}\right) $ is the $ARMA$ field model
where $\mathbb{Z}_{1}^{d}$ and $\mathbb{Z}_{2}^{d}$ are finite subsets of $ \mathbb{Z}^{d}$ and henceforth $z^{j}=\prod\nolimits_{\ell =1}^{d}z\left[ \ell \right] ^{j\left[ \ell \right] }$ with the convention that $0^{0}=1$. As an example, we have the $ARMA\left( -k_{1},k_{2};-\ell _{1},\ell _{2}\right) $ field
As mentioned above, the Wold decomposition, and hence the concept of past and future, plays a key role in the theory of prediction when $d=1$. However, contrary to the situation when $d=1$, an intrinsic problem with spatial or lattice data is that we cannot assign a unique meaning to the concept of \textquotedblleft past\textquotedblright\ and/or \textquotedblleft future\textquotedblright . One immediate consequence is then that different definitions of what might be considered as past (or future) will yield different Wold-type decompositions. More specifically, denote a \textquotedblleft half-plane\textquotedblright\ of $\mathbb{Z}^{2}$ according to the lexicographical (dictionary) ordering \textquotedblleft $ \prec $\textquotedblright\ defined as
where herewith we shall consider the case when $d=2$, often encountered with real data. The half-plane defined by \textquotedblleft $\prec $ \textquotedblright\ is illustrated in Figure (ref). Following earlier work by Helson1958,Helson1961, there exists then a Wold-type representation of the (spatial) process $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ given by
where $\left\{ \vartheta _{t}\right\} _{t\in \mathbb{Z}^{2}}$ is a zero mean white noise sequence with finite second moments $\sigma _{\vartheta }^{2}$. It is worth recalling that $\vartheta _{t}$ once again has the interpretation of being the \textquotedblleft one-step\textquotedblright\ prediction error. Often $\left( \ref{uni_1}\right) $ is called a unilateral representation of $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ as opposed to the multilateral representation in $\left( \ref{a1}\right) $. See also Whittle1954 for some earlier work on multilateral versus unilateral representations. As an example, $\left( \ref{arma}\right) $ becomes a unilateral or causal model when $\ell _{1}=k_{1}=0$. $\left( \ref{uni_1} \right) $ might be regarded as a particular way to model the dependence of $ x_{t}$ induced by the lexicographic ordering in $\left( \ref{lex_1}\right) $ . Of course, the choice of the \textquotedblleft half-plane\textquotedblright\ of $\mathbb{Z}^{2}$ according to the associated chosen lexicographic ordering is not the only possible one. That is, a different choice of \textquotedblleft half-plane\textquotedblright\ of $\mathbb{Z}^{2}$, induced by the lexicographic ordering, will yield a \textquotedblleft similar\textquotedblright\ but different representation of $x_{t}$ to that given in $\left( \ref{uni_1}\right) $. As it will become clear in the next section, the choice of a specific lexicographic ordering, or its associated half-plane, will depend very much on practical purposes. For instance, the choice of $\left( \ref{lex_1}\right) $ will depend on the location where we wish to predict $x_{t}$. Last but not least it is worth, and important, mentioning that the sequences $\left\{ \varepsilon _{t}\right\} _{t\in \mathbb{Z}^{2}}$ and $\left\{ \vartheta _{t}\right\} _{t\in \mathbb{Z}^{2}}$ are not the same. Recall that a similar phenomenon occurs when $d=1$ and the practitioner allows for noncausal/bilateral representations of the sequence $x_{t}$. When this is the case, the \textquotedblleft bilateral or noncausal\textquotedblright\ representation has errors which are independent and identically distributed, whereas for its \textquotedblleft unilateral or causal\textquotedblright representation, the corresponding errors are only a white noise sequence.
It is clear from the introduction that to provide accurate and valid (linear) predictions (or interpolations), a key component is to obtain the covariogram function of the sequence $\left\{ x_{t}\right\} _{t\in \mathbb{Z} ^{2}}$, that is $C\left( h\right) =Cov\left( x_{t},x_{t+h}\right) $, which is related to the spectral density function $f\left( \lambda \right) $ via the expression
where $\Pi =\left( -\pi ,\pi \right] $. Henceforth the notation \textquotedblleft $h\cdot \lambda $\textquotedblright\ means the inner product of the vectors $h$ and $\lambda $. It is worth observing that we can factorize $f\left( \lambda \right) $ as
where $\sigma _{\varepsilon }^{2}=E\varepsilon _{t}^{2}$ and $\sigma _{\vartheta }^{2}=E\vartheta _{t}^{2}$, and
The latter displayed expressions indicate that either $\Psi \left( \lambda \right) $ or $\Upsilon \left( \lambda \right) $ summarize the covariogram structure of $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$.
When $d=1$ and the sequence $\left\{ x_{t}\right\} _{t\in \mathbb{Z}}$ is purely nondeterministic we know, see Whittle $\left( 1961\text{, } p.26\right) $ or Brillinger $\left( 1981\text{, Theorem }3.8.4\right) $, that the spectral density $f\left( \lambda \right) $ admits a representation
where by definition $A\left( \lambda \right) =:\exp \left\{ -\sum_{k=1}^{\infty }\alpha _{k}e^{ik\cdot \lambda }\right\} $. The latter expression is referred to as the canonical factorization of the spectral density function and is also known as Bloomfield's model. One important consequence of the canonical factorization is that the sequence $\left\{ x_{t}\right\} _{t\in \mathbb{Z}}$ can be written as
where $\vartheta _{t}$ is a zero mean white noise sequence with finite second moments and $a_{j}$ are the Fourier coefficients of $A\left( \lambda \right) $, that is
with $2\pi \exp \left( \alpha _{0}\right) =\sigma _{\vartheta }^{2}$, i.e. the one-step prediction error. However, more importantly, denoting
we have that its Fourier coefficients equal the coefficients $\zeta _{j}$ in $\left( \ref{uni_1}\right) $.
Whittle1954, Section 6, signalled that a similar argument can be used when $d>1$. However a formal and theoretical justification for a canonical factorization of $f\left( \lambda \right) $ when $d>1$ was discussed in Korezlioglu1986, see also Solo1986. More specifically, they show that the spectral density function of $\left\{ x_{t}\right\} _{t\in \mathbb{Z }^{2}}$ might be characterized using the representation
where
which is sometimes known as the Cepstrum model by Solo1986, who notes that if $0<f\left( \lambda \right) <M$ then the representation of the spectral density in $\left( \ref{bloom_1}\right) $ or in $\left( \ref {arbloompa}\right) $ exists, see also mcelroy2014. Note that the coefficients $\alpha _{k}$ in $\left( \ref{arbloompa}\right) $ are the Fourier coefficients of $\log \left( f\left( \lambda \right) \right) $, that is
where $\widetilde{\Pi }^{2}=\left[ 0,\pi \right] \times \Pi $, that is $ \lambda \in \widetilde{\Pi }^{2}$ if $\lambda \left[ 1\right] \in \left[ 0,\pi \right] $ and $\lambda \left[ 2\right] \in \Pi $.
As it is the case when $d=1$, there is a relationship between the representation in $\left( \ref{uni_1}\right) $ and $\left( \ref{bloom_1} \right) /\left( \ref{arbloompa}\right) $, i.e. between the coefficients $ \zeta _{j}$ and $\alpha _{k}$. So, it will be convenient to discuss the relationship between the representations of the sequence $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ in the \textquotedblleft frequency\textquotedblright\ and \textquotedblleft space\textquotedblright\ domains. The link among these coefficients turns out to play a crucial role in our prediction algorithm. For that purpose, consider the lexicographic ordering given in $\left( \ref{lex_1}\right) $. Then, denoting the Fourier coefficients of $A\left( \lambda \right) $ by
and $a_{0}=1$, the sequence $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ has a unilateral representation given by
where $\left\{ \vartheta _{t}\right\} _{t\in \mathbb{Z}^{2}}$ is the sequence given in $\left( \ref{uni_1}\right) $. But also we have that the coefficients $\zeta _{j}$ in $\left( \ref{uni_1}\right) $ are the Fourier coefficients of $B\left( \lambda \right) =:A^{-1}\left( \lambda \right) =\exp \left\{ \sum_{0\prec k}\alpha _{k}e^{-ik\cdot \lambda }\right\} $. That is,
see Section 1.2 of Korezlioglu1986. The latter might be considered as an extension of the canonical factorization given in Brillinger1981 to the case $d>1$. However, one key aspect is that there is a direct link between $\alpha _{k}$ and the coefficients of the Wold-type decomposition of its autoregressive representation, that is $a_{j}/\zeta _{j}$ and $\alpha _{k}$. This observation will be important for our prediction methodology in the next section.
The purpose of the section is to present and examine a prediction algorithm, extending the methodology in Bhansali1974 or Hidalgo2002, to the case when $d=2$. Similar to the aforementioned work, a key component of the methodology will be based on the canonical factorization of the spectral density in $\left( \ref{bloom_1}\right) $. Due to the rather unusual notation in this paper, we have decided to collate it at this stage for convenience. Given two vectors $a$ and $b$, $a\geq \left( \leq \right) b$ means that $a\left[ \ell \right] \geq \left( \leq \right) b\left[ \ell \right] $ for all $\ell =1,2$. Denote
where $\lambda _{k}=\left( \lambda _{k\left[ 1\right] },\lambda _{k\left[ 2 \right] }\right) $ are the Fourier frequencies and $\widetilde{\Pi } _{n}^{2}=\left\{ \lambda _{k}\in \Pi _{n}^{2}:\lambda _{k\left[ 1\right] }>0\right\} $. Finally, we denote
Similarly, we denote
where we are using the convention that for any $k\in \mathbb{Z}^{2}$, we write $d_{k}$ as
Observe that $\sum_{j\preceq J}^{+}+\sum_{j\preceq J}^{-}=\sum_{-J<j\leq J}$ , and likewise $\int_{\lambda \preceq \pi }^{+}+\int_{\lambda \preceq \pi }^{-}=\int_{\lambda \in \Pi ^{2}}$.
Before we describe our prediction algorithm, we shall introduce our set of regularity conditions.
We now comment on Conditions $C1$ and $C2$. First, Condition $C2$ can be generalized to allow for different rates of convergence to zero of $n^{-1} \left[ \ell \right] $, $\ell =1,2$. However, for notational simplicity, we prefer to keep it as it stands. Condition $C1$ could have been written in terms of the multilateral representation in $\left( \ref{a1}\right) $. However since the prediction employs the representation in $\left( \ref{SAR} \right) $ or $\left( \ref{uni_1}\right) $, we have opted to write $C1$ as it stands. Part $\left( a\right) $ of Condition $C1$ seems to be a minimal condition for our results below to hold true. Sufficient regularity conditions required for the validity of the expansion in $\left( \ref{SAR} \right) $ is $\Upsilon \left( z\right) $ be nonzero for any $z\left[ \ell \right] $, $\ell =1,2$. The latter condition guarantees that $f\left( \lambda \right) >0$ for all $\lambda \in \widetilde{\Pi }^{2}$. Part $\left( \mathbf{c}\right) $ entails that the spectral density $f\left( \lambda \right) $ is $4$ times continuously differentiable. This is needed if one wants to achieve a similar rate of approximation of sums by their integrals when $d=1$ and the function is twice continuously differentiable. Indeed whereas when $d=1$, we have that
with two continuous derivatives for $g\left( x\right) $, to have a \textquotedblleft similar\textquotedblright\ result when $d=2$ one needs $ g(x)$ to be $4$ times continuously differentiable. See Lemma (ref) in the appendix for some extra insight.
We now discuss the methodology to predict the\ value of $x_{t}$ at an unobserved location without imposing any specific parametric model for $ f\left( \lambda \right) $. In addition, as a by-product, we provide a simple estimator of the coefficients $\zeta _{j}$ or $a_{j}$. First, $A\left( \lambda \right) $ and expression $\left( \ref{alpha_1}\right) $ suggest that to compute an estimator of the coefficients $\alpha _{j}$ and/or $a_{j}$, it suffices to obtain an estimator of $f\left( \lambda \right) $. To that end, for a generic sequence $\left\{ v_{t}\right\} _{t=1}^{n}$, we shall define the discrete Fourier transform, $DFT$, as
and the periodogram as
where, in what follows, we use the notation that for any $g=\left( g\left[ 1 \right] ,g\left[ 2\right] \right) $,
In real applications, in order to make use of the fast Fourier transform, the periodogram will be evaluated at the Fourier frequencies $\lambda _{k}$.
However as noted by Guyon1982, due to non-negligible end effects (the edge effect), the bias of the periodogram does not converge to zero fast enough when $d>1$. We therefore proceed as in Dahlhaus1987, and employ the tapered periodogram defined as
where $w_{v}^{T}\left( \lambda _{j}\right) $ denotes the taper discrete Fourier transform, $DFT$. One common taper is the cosine-bell (or Hanning) function, which is defined as
see Brillinger1981. It is worth observing the cosine-bell taper DFT is related to $w_{v}\left( \lambda \right) $ by the equality
In this paper we shall explicitly consider the cosine-bell, although the same results follow employing other taper functions such as Parzen or Kolmogorov tapers Brillinger1981. \ This is formalized in the next condition.
Using notation in $\left( \ref{gblack}\right) $, we shall estimate $f\left( \lambda \right) $ by the average tapered periodogram
where $m\left[ \ell \right] /n\left[ \ell \right] +m\left[ \ell \right] ^{-1}=o\left( 1\right) $,$\ $for $\ell =1,2$. Next, we denote $\widetilde{ \lambda }_{k}=\left( \widetilde{\lambda }_{k\left[ 1\right] },\widetilde{ \lambda }_{k\left[ 2\right] }\right) ^{\prime }$, for $k\left[ 1\right] =0,1,...,M\left[ 1\right] =:\tilde{n}\left[ 1\right] /m\left[ 1\right] $ and $k\left[ 2\right] =0,\pm 1,...,\pm M\left[ 2\right] =:\tilde{n}\left[ 2 \right] /m\left[ 2\right] $, where
Bearing in mind $\left( \ref{notd}\right) $, denoting $\mathcal{M=}\left\{ j:~\left( 0\prec j;j=0\right) \text{ }\wedge \left( -M<j\leq M\right) \right\} $ and abbreviating $\phi \left( \widetilde{\lambda }_{k}\right) $ by $\phi _{k}$ for a generic function $\phi \left( \lambda \right) $, we estimate the coefficients $a_{j}$, $j=1,...,M$, as
It is also worth defining the quantities $\left( \ref{ahat_j}\right) $ and $ \left( \ref{cr_1}\right) $ when $\widehat{f}\left( \lambda \right) $ is replaced by $f\left( \lambda \right) $, that is
That is,
and also we denote
We shall now begin describing how we can predict a value $x_{t}$ at the location $s=\left( s\left[ 1\right] ,s\left[ 2\right] \right) $ such that $ 1\leq s\left[ 1\right] \leq n\left[ 1\right] $ and $1\leq s\left[ 2\right] \leq n\left[ 2\right] $. For instance, we wish to predict the unobserved value\ $x_{s}$
Now, the location of $s$ suggests that a convenient unilateral representation of $x_{t}$ appears to be
which comes from the lexicographic ordering in $\left( \ref{lex_1}\right) $. Since we need to estimate the coefficients $a_{k}$, the prediction will then become
where $\widehat{a}_{k}~x_{s-k}=:\widehat{a}_{k\left[ 1\right] ,k\left[ 2 \right] }~x_{s\left[ 1\right] -k\left[ 1\right] ,s\left[ 2\right] -k\left[ 2 \right] }$. However, it may be very plausible that the value$\ $of $M$ is such that we may not observe the process at some of the locations employed to compute $\left( \ref{1}\right) $. That is, consider the situation where we want to predict $x_{s}$
In this case we observe that to compute $\left( \ref{1}\right) $, we first need to obtain a predictor of values of $x_{s-k}$ when say $k\left[ 1\right] =1$ and $k\left[ 2\right] <0$, since $x_{s-k}$ is not observed at those locations, which in its computation needs predictors of the relevant values themselves. See $\left( \ref{pre_1}\right) $ for more exact details. However, in this case one can avoid this extra computational burden. Indeed, this is so as the relative location $\left( s\left[ 1\right] ,s\left[ 2 \right] \right) $ suggests that the practitioner might have employed the Wold-type representation
which can be regarded as induced by the lexicographic ordering\
Note that the lexicographic ordering $\left( \ref{lex_2}\right) $ is as that in $\left( \ref{lex_1}\right) $ but swapping $j\left[ 2\right] $ for $j\left[ 1\right] $. From here, we proceed as with $\left( \ref{1}\right) $ but with the \textquotedblleft coordinates\textquotedblright\ $\left[ 2\right] $ and $ \left[ 1\right] $ changing their roles.
Finally, consider the case where location we wish to predict $x_{s}$ is $ \left( n\left[ 1\right] +1,s\left[ 2\right] \right) $. That is,
Now, the location of $s=:\left( n\left[ 1\right] +1,s\left[ 2\right] \right) $ suggests that the more convenient representation of $x_{s}$ appears to be that in $\left( \ref{unil_1}\right) $ which comes from the lexicographic ordering in $\left( \ref{lex_1}\right) $, and hence our prediction is given in $\left( \ref{1}\right) $. That is, since we need to estimate the coefficients $a_{k}$, the prediction will then become
However to compute the prediction we also need to replace the unobserved $ x_{s}$ by its prediction. As with \textquotedblleft standard\textquotedblright\ time series when we wish to predict beyond $1$ period ahead, this is done by recursion, that is we make use of formula $ \left( \ref{1}\right) $ starting say from the value $x_{n\left[ 1\right] +1,s \left[ 2\right] -M\left[ 2\right] }$. Once we have \textquotedblleft predicted\textquotedblright\ the value for this observation, we then predict $x_{n\left[ 1\right] +1,s\left[ 2\right] -M\left[ 2\right] +1}$ and so on. For instance, for any $r\left[ 1\right] =0,...,r$ and $r\left[ 2\right] =0,...,r=\min \left\{ n\left[ 2\right] /8;M\left[ 2\right] \right\} $,
where we take the convention that $\widehat{x}_{s}=x_{s}$ if the location were observed and $=:0$ when $s\left[ 2\right] <-r$ or $\left\{ s\left[ 2 \right] <0\text{ }\wedge s\left[ 1\right] <n\left[ 1\right] -r\right\} $. Finally, if we were interested to predict $x_{t}$ at the unobserved location $\left( t\left[ 1\right] ,n\left[ 2\right] +1\right) $, then it suggests to employ the lexicographic ordering in $\left( \ref{lex_2}\right) $ and hence the representation given in $\left( \ref{unil_2}\right) $, and then we would proceed as above but again with the \textquotedblleft coordinates\textquotedblright\ $\left[ 2\right] $ and $\left[ 1\right] $ changing their roles.
Before we examine the statistical properties of $\widehat{x}_{t}$ in $\left( \ref{1}\right) $ or $\left( \ref{1a}\right) $, we shall look at those of $ \widehat{\alpha }_{j}$ or $\widehat{A}_{j}$. For that purpose, denote
Also, denote $\left\{ \boldsymbol{\xi }_{j}\right\} _{j}$ the Fourier coefficients of $g\left( \lambda \right) $ \ given by
Notice that Condition $C1$ implies that $g\left( \lambda \right) $ is twice continuous differentiable, so that $\left\{ \boldsymbol{\xi }_{j}\right\} _{j}$ is summable.
We introduce one extra condition relating the rate of increase of $m\left[ \ell \right] $ with respect of $n\left[ \ell \right] $.
We shall now denote $a_{\upsilon }=0$ if $\upsilon \prec 0$.
Once we have obtained the asymptotic properties of the estimators of $a_{j}$ , for $0\prec j$ and $j\leq M$, we are in a position to examine the asymptotic properties of the predictor $\widehat{x}_{s}$ in $\left( \ref{1} \right) $ or $\left( \ref{1a}\right) $. To that end, denote by $\left\{ x_{t}^{\ast }\right\} _{t\in \mathbb{Z}^{2}}$ a new independent replicate sequence with the same statistical properties of the original sequence $ \left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ not used in the estimation of the spectral density function. Then let $\widehat{x}_{s}^{\ast }$ be as in $ \left( \ref{1}\right) $ but with $\widehat{x}_{t}$ replaced by $x_{t}^{\ast } $ there, that is
or $\left( \ref{1a}\right) $ but with $\widehat{x}_{t}$ being replaced by $ x_{t}^{\ast }$ there, that is
We examine the finite-sample behaviour of our algorithm in a set of Monte Carlo simulations. As in Robinson2006 and robinson2007nonparametric we used the model
similar to one considered in Haining1978. Then
with $\nu \left( \lambda \right) =\prod_{j=1}^{2}\left( 1+2\cos \lambda _{j}\right) -1$. Robinson2006 show that a sufficient condition for invertibility of ((ref)) is
We first generated a $40\times 41$ lattice using ((ref)), with $ \tau =0.05, 0.075, 0.10$ and the $\epsilon _{t}$ drawn independently from three different distributions for each $\tau$: $U(-5,5)$, $N(0,1)$ and $\chi _{9}^{2}-9$. The aim of this section is to examine the performance of both prediction algorithms in predicting the $20,20$-th element of this lattice. We did this by assuming a situation in which the practitioner has available data sets of various sizes, generated from ((ref)). To permit a clear like-for-like comparison of improvement in performance as sample size increases, we construct the prediction coefficients using the samples generated in each replication and then use these to construct predictions for the 20,20-th element of the $40\times 41$ lattice.
We took $n[1]=n^{\ast }+1$ and $n[2]=2n^{\ast }+1$, for some positive integer $n^{\ast }$, implying $\mathbf{n}=\left( 2n^{\ast }+1\right) (n^{\ast }+1)$, and generated iid $\epsilon_t$ from each of the three distributions mentioned in the previous paragraph. In each of the 1000 replications we experimented with $\tau =0.05,0.075,0.10$ and $n^{\ast }=5,10,20$ and $40$. The choices of $\tau $ satisfy ((ref)).
Given the different sample sizes in each dimension, we can experiment with more values of $m[1], m[2]$ and $p_1,p_2$ as $n^*$ increases. We make the following choices:
The flexible exponential approach requires a nonparametric estimate of $ f(\lambda)$. Two such estimates are available to use: the first one based on the tapered periodogram described in ((ref)), which we denote $\hat{f} (\lambda)$, and the second based on the autoregressive approach in Gupta2016. The latter also provides a rival prediction methodology based on a nonparametric algorithm using AR model fitting, extending well established results for $d=1$, see Bhansali1978 and Lewis1985. The idea is first to obtain a least squares predictor based on a truncated autoregression of order $p=\left( p_{L_{1}},p_{U_{1}};p_{L_{2}},p_{U_{2}}\right) $, for non-negative integers $ p_{L_{\ell }},\;p_{U_{\ell }}$, $\ell =1,2$, with the truncation allowed to diverge as $N\rightarrow \infty $. That is, we approximate the infinite unilateral representation in ((ref)) by one of increasing order.
In view of the half-plane representation we can a priori set, say, $ p_{L_{2}}=0$ when considering $\preccurlyeq $. If we could observe the AR prediction coefficients $a_{k}$, say, a prediction of $x_{s}$ based on $ \preccurlyeq $ could be constructed as
where $S\left[ -p_{L},p_{U}\right]$ is the intersection of the set $\left\{ t\in \mathbb{L}:-p_{L_{\ell }}\leq t_{\ell }\ \leq p_{U_{\ell }},\;\ell =1,2\right\}$ with the prediction half-plane. This is the spatial version of one-step prediction and again we follow the convention that $\check{x}_s=x_s$ if $x_s$ is observed. However ((ref)) is not feasible and needs to be replaced by an approximate version, as described below.
Writing $p_{\ell }=p_{L_{\ell }}+p_{U_{\ell }}$, we assume throughout that $ n[\ell ]>p_{\ell }$ for $\ell =1,2$, and denote $n_{p}=\prod_{\ell =1}^{2}\left( n[\ell ]-p_{\ell }\right) $, $\mathfrak{h}(p)=p_{U_{2}}+\left( p_{1}+1\right) p_{U_{2}}$ , i.e. the cardinality of $S\left[ -p_{L},p_{U} \right] $. Suppose that the data are observed on $\left\{ \left( t_{1},t_{2}\right) :n_{L_{1}}\leq t_{1}\leq n_{U_{1}},-n_{L_{2}}\leq t_{2}\leq n_{U_{2}}\right\} $. Define a least squares predictor of order $ \mathfrak{h}(p)$ by
where $\sum_{j(p,n)}^{\prime \prime }$ runs over $\left\{ \left( j_{1},j_{2}\right) :p_{1}-n_{L_{1}}<j_{1}\leq n_{U_{1}}+1,p_{2}-n_{L_{2}}<j_{2}\leq n_{U_{2}}+1\right\} $. We denote the elements of $\check{d}_{p}$ by $\check{d}_{p}(k)$, $k\in S\left[ -p_{L},p_{U} \right] $, and the minimum value by $\check{\sigma}_{p}^{2}$. A feasible half-plane prediction based on a fitted autoregression of order $p$ is given by
The autoregressive nonparametric spectrum estimate is defined as
A predictor of $x_{s}$ based on $\left( \ref{1}\right) $ using $\hat{f} (\lambda )$ (respectively $\check{f}(\lambda )$) is denoted $\hat{x}_{s}$ (respectively $\tilde{x}_{s}$), while a predictor based on ((ref) ) is denoted $\check{x}_{s}$ as mentioned above.
Let $\vec{x}_{r,s}$ be a generic predictor of $x_{s}$ in replication $r$, $ r=1,\ldots ,1000$. We report a statistic called the root mean squared error (RMSE) of prediction, defined as
The results are reported in Tables (ref)-(ref). We observe an improvement in prediction performance as $ n^*$ increases, and also as the bandwidths ($(m[1],m[2])$ and $p^*$) increase as function of $n^*$. This is as expected in the theory. Nevertheless, even for rather small sample sizes the RMSE is acceptable. For example, for $\epsilon_t\sim U(-5,5)$ and $\epsilon_t\sim N(0,1)$ with $ n^*=5 $, we can obtain predictions with RMSE that are not radically different from the $n^*=10$ case, even though this change in $n^*$ entails a sample that is nearly four times larger (231 against 66). In comparison the RMSE with the smaller sample size can be quite close to those obtained with more data in some cases, cf. $\check x _{20,20}$ for any error distribution.
For the smallest sample size $\check x_{20,20}$ can outperform $\hat x_{ 20,20}$ and $\tilde x _{20,20}$, but with increasing $n^*$ the latter two clearly begin to dominate. An inspection of Tables (ref)-(ref) reveals that the use of the flexible exponential algorithm proposed in this paper together with either the tapered periodogram or the AR spectral estimator of Gupta2016 outperforms autoregressive prediction in moderate to large sample sizes. There is little to choose from between the two best performing algorithms, and a practitioner might choose to use either one. However the AR prediction is clearly dominated by our algorithm.
In this section we show how the techniques established in the paper can be used to predict house prices. This can be of interest in real estate and urban economics, as well as for property developers. Indeed, spatial methods are frequently used in these fields, as studied for instance by IversenJr2001, Banerjee2004 and Majumdar2006. We use median house price data for census blocks in California from the 1990 census from Pace1997a, available at www.spatial-statistics.com. We confine our analysis to the city of Los Angeles. The data is gridded as follows: a $ 14\times 23$ grid of square cells is superimposed on Los Angeles, from $ 33.75^{\circ}$N to $34.17^{\circ}$N and $117.75^{\circ}$W to $118.44^{\circ}$ W. The grid covers a total of 5259 observations. The average of the median house values for each cell is calculated and the 322 such observations form our sample. The gridding is shown in Figure (ref), in which the 8 empty cells are filled and marked with a cross. We wish to predict the house price for these cells. House price data is not a zero mean process, so we subtract the sample mean using the whole sample from each cell.
We proceed in the following way: to obtain the coefficients $\hat{a}_{\ell }$ , and $\check{d}(\ell )$ in $\left( \ref{1}\right) $ and ((ref)) we use the $14\times 19$ sublattice formed of the first 19 columns of cells. This sublattice contains no missing observations. Once the coefficients are obtained we construct predictions using the remaining $4\times 19$ sublattice, in a step-by-step manner. The shaded-and-crossed cell (8,20) is predicted first, followed by (8,21) and (8,22). We then predict (4,21), followed by (7,22), (9,23), (6,23) and (1,23).
The predicted values are tabulated for various values of $\left( m[1],m[2]\right) $ and $\left( p_{1},p_{2}\right) $ in Tables (ref) and (ref). The predicted values are quite stable across the choices $\left( m[1],m[2]\right) =(2,2),(2,3)$ using either the periodogram or AR spectral estimate. They most closely match those obtained when $\left( p_{1},p_{2}\right) =(2,2),(2,3)$ in ((ref)). In the latter case we compare in Table (ref) the order selection criteria proposed by Gupta2016, which include the usual FPE and BIC (denoted with a $\widehat{}$ ) as well as corrected version that account for the spatial case (denoted with $\widetilde{}$ and $\bar{}$ ). The FPE tends to favour longer lag lengths no matter which version is used, as do $\widehat{\text{BIC}}$ and $ \overline{\text{BIC}}$. However the latter as well as $\widetilde{\text{FPE}} $ are not monotonically decreasing in lag length, unlike $\widehat{BIC}$, $ \widehat{FPE}$ and $\overline{FPE}$. Thus the latter three are likely to overfit and seem undesirable. If we impose a selection rule that picks the desirable lag order as the first instance when the selection criteria shows an increase with lag length, then we get $\left( p_{1},p_{2}\right) =(2,2)$ using $\widetilde{FPE}$ and $\overline{BIC} $. $\widetilde{BIC}$ indicates a choice of $\left( p_{1},p_{2}\right) =(1,2)$, on the other hand. All considered, it seems that $\left( p_{1},p_{2}\right) =(2,2)$ is a reasonable choice.
In this paper we have dealt with the problem of prediction when the data $ \left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ is collected on a lattice. To do so, we considered unilateral representations of $\left\{ x_{t}\right\} _{t\in \mathbb{Z}^{2}}$ and in particular the canonical factorization of the spectral density function, the latter being possible as observed by Whittle1954. Our approach does not need any parameterization of the model (i.e. the covariogram structure of the data), so we avoid the consequences that a wrong parameterization can have in the predictor. We have also compared our methodology to one based on the space domain by using a finite approximation of the unilateral autoregressive model in $\left( \ref{SAR} \right) $.
However, it might be interesting to examine how our proposed methodology compares with one based on the conditional autoregressive ($CAR$) representation of Besag1974. That is, let $x_{t}$ be given by
Note that our definition in $\left( \ref{CAR}\right) $ implies that $x_{t}$ is, among other characteristics, homogeneous. The representation of $x_{t}$ given in $\left( \ref{CAR}\right) $ suggests to predict a value $x_{t}$ at a location $s=\left( s\left[ 1\right] ,s\left[ 2\right] \right) $, $1\leq s \left[ 1\right] \leq n\left[ 1\right] $ and $1\leq s\left[ 2\right] \leq n \left[ 2\right] $, by
where $\widehat{\mu }$ and $\widehat{\zeta }_{\left\vert r-s\right\vert }$ are respectively the least squares estimator of $\mu $ and $\zeta _{\left\vert r-t\right\vert }$, and with the convention that $x_{r}=0$ if it were not observed. This is in the same spirit as we did with our predictor in $\left( \ref{1}\right) $. On the other hand, if we were interesting to predict a value $x_{t}$ at a location $s=\left( n\left[ 1\right] +1,s\left[ 2 \right] \right) $, we might then use
However to compute the prediction we would also need to replace the unobserved $x_{r}$ by its prediction as in $\left( \ref{1a} \right) $. The latter might be done in an iterative fashion similar to what we did in $\left( \ref{pre_1}\right) $.
.