# Dynamic Spatiotemporal ARCH Models

Philipp Otto
Leibniz University Hannover, Germany
Osman Doğan
Department of Economics, University of Illinois at Urbana-Champaign, U.S.A.
Süleyman Taşpınar
Department of Economics, Queens College CUNY, New York, U.S.A.

###### Abstract

Geo-referenced data are characterized by an inherent spatial dependence due to
the geographical proximity. In this paper, we introduce a dynamic
spatiotemporal autoregressive conditional heteroscedasticity (ARCH) process to
describe the effects of (i) the log-squared time-lagged outcome variable,
i.e., the temporal effect, (ii) the spatial lag of the log-squared outcome
variable, i.e., the spatial effect, and (iii) the spatial lag of the log-
squared time-lagged outcome variable, i.e., the spatiotemporal effect, on the
volatility of an outcome variable. Furthermore, our suggested process allows
for the fixed effects over time and space to account for the unobserved
heterogeneity. For this dynamic spatiotemporal ARCH model, we derive a
generalized method of moments (GMM) estimator based on the linear and
quadratic moment conditions of a specific transformation. We show the
consistency and asymptotic normality of the GMM estimator, and determine the
best set of moment functions. We investigate the finite-sample properties of
the proposed GMM estimator in a series of Monte-Carlo simulations with
different model specifications and error distributions. Our simulation results
show that our suggested GMM estimator has good finite sample properties. In an
empirical application, we use monthly log-returns of the average condominium
prices of each postcode of Berlin from 1995 to 2015 (190 spatial units, 240
time points) to demonstrate the use of our suggested model. Our estimation
results show that the temporal, spatial and spatiotemporal lags of the log-
squared returns have statistically significant effects on the volatility of
the log-returns.

Keywords: Spatial ARCH, GMM, volatility clustering, volatility, house price
returns, local real-estate market

##  1 Introduction

In a standard autoregressive conditional heteroskedasticity (ARCH) model, the
volatility is modeled as a linear function of the lagged squared outcome
variable in order to account for the volatility clustering patterns observed
in the outcome variable (Engle,, 1982; Bollerslev et al.,, 1992; Engle and
Bollerslev,, 1986). However, when analyzing geo-referenced time series, a
further phenomenom occurs – observations close in space are more similar than
observations farther away – known as Tobler’s first law of geography (Tobler,,
1970). From a statistical perspective, this spatial dependence may occur in
the means and the volatility of a random process (Otto et al.,, 2018). Thus,
in this paper, we extend the standard ARCH model to spatiotemporal data by
using some tools from spatial econometrics. In our suggested specification,
the log-volatility term may depend on (i) the log-squared time-lagged outcome
variable, (ii) the higher-order spatial lags of the log-squared outcome
variable, (iii) the higher-order spatial lags of the log-squared time-lagged
outcome variable, (iv) exogenous variables, and (v) the unobserved
heterogeneity across regions and time. The presence of higher-order spatial
lags in our specification indicates that the log-volatility term of a region
may depend on the current and time-lagged outcome variables in the neighboring
locations in differing ways, depending on the specifications of the associated
spatial weights matrices (e.g., different influences from different directions
or directional dependence, cf. Gupta and Robinson, 2015; Merk and Otto, 2021).
We refer to this extended model as the dynamic spatiotemporal ARCH model.

To introduce an estimation approach for our model, we transform the outcome
equation so that it is in the form of log-squared terms. We then substitute
the log-volatility equation into the the transformed outcome equation to
obtain an estimation equation for the log-squared outcome variable. The
resulting specification is in the form of a higher-order spatial dynamic panel
data model with disturbance terms that may not have a zero mean. We use an
orthonormal and a deviation from group-mean operator to wipe out the regional
and time fixed effects from the estimation specification. For the estimation
of the transformed model that is free of the regional and time fixed effects,
we propose a generalized method of moments (GMM) estimator formulated with a
set of linear and quadratic moment functions (Lee,, 2007; Lee and Liu,, 2010;
Lee and Yu,, 2014). We show that the resulting GMM estimator has the standard
large sample properties irrespective of whether the number of time periods is
large or finite. When the number of time periods is large, the precision
matrix of our GMM estimator simplifies significantly, allowing us to determine
a set of linear and quadratic moment functions that can lead to an efficient
estimator. We provide such a set of best linear and quadratic moment
functions, and establish the asymptotic properties of the resulting best GMM
estimator. In a Monte Carlo simulation study, we show that the proposed GMM
estimator performs well in finite samples.

In the literature, Robinson, (2009) introduces the log-square transformation
to a cross-sectional spatial stochastic volatility model, and consider a quasi
maximum likelihood (QML) estimation approach for the estimation of the
transformed model.111In the literature, the log-square transformation approach
is also used for the estimation of (i) the cross-sectional spatial stochastic
volatility models (Taşpınar et al.,, 2021), and (ii) the cross-sectional
spatial ARCH/GARCH models (Sato and Matsuda,, 2017; Otto,, 2019; Otto and
Schmid,, 2020; Sato and Matsuda,, 2021). In a similar manner, we may
alternatively consider the QML estimation approach rather than the GMM
approach for our transformed model (Yu et al.,, 2008; Lee and Yu,, 2010;
Hølleland and Karlsen,, 2020). Compared to the QML estimation approach, our
GMM estimation approach has the following advantages. First, the GMM approach
has the computational advantage over the QML approach since the QML estimation
involves calculation of the determinant of a Jacobian term at each iteration
during the estimation. The computational cost can be especially high when the
number of the cross-sectional units is large222LeSage and Pace, (2009) provide
some solutions based on various approximation methods to reduce the
computation cost significantly. Second, it is well known that the QML
estimator has an asymptotic bias, and thus requires a bias correction approach
even when the number of time period is large (Yu et al.,, 2008; Lee and Yu,,
2010). Finally, the QML estimator may have poor finite sample properties since
the distribution of the log-squared disturbance terms in the transformed model
is approximated by a normal distribution. In the time series literature on the
volatility models, it has been documented that the QML estimator obtained in
this way has poor finite sample properties (Jacquier et al.,, 1994; Shephard,,
1994; Kim et al.,, 1998; Sandmann and Koopman,, 1998).

In an empirical application, we use a monthly dataset of the real house price
returns in Berlin at the zip-code level over the period, January 1995 to
December 2015 to test the effect of temporal, spatial and spatiotemporal lags
of the log-squared returns on the volatility of the log-returns. That is, we
analyze the volatility of the house price returns in a real-estate market on a
small geographic area of around 900 $km^{2}$ (see also McMillen, 2014; Billé
et al., 2017; Zhang and Yi, 2017). In Section 2, we show that our dynamic
spatiotemporal ARCH model implies a spatial dynamic panel data model for the
log-squared returns. Therefore, the presence of spatial, temporal and
spatiotemporal effects in the log-squared returns will provide the empirical
evidence for our suggested specification. To motivate the presence of these
effects on the log-squared returns, Figure 1 displays the average log-squared
returns over Berlin’s zip-codes (the left figure), the estimated temporal
autocorrelation of the log-squared returns as a series of boxplots (the center
figure), and the estimated spatiotemporal autocorrelation in terms of Moran’s
$I$ across the time horizon (the right figure). The first figure shows a
clustering pattern in the log-squared returns, indicating the presence of
spatial dependence. From the ACF estimates, we can observe a clear temporal
volatility clustering, while the spatiotemporal dependence is of minor degree,
irregularly fluctuating around zero.

By using a first-order version of our dynamic spatiotemporal ARCH process for
the local house price returns, we separately identify temporal, spatial and
spatiotemporal interaction effects in the log-squared returns. Our estimation
results show that the temporal, spatial and spatiotemporal lags of the log-
squared returns have statistically significant effects on the log-volatility,
and for that there is significant variation in the log-volatility of the real
house price returns in Berlin over its zip-codes. This finding is not
surprising, because it has been documented in the literature that the spatial
dependence in house price variations might arise due to several factors such
as migration, equity transfer, spatial arbitrage and spatial patterns in the
determinants of house prices (Meen,, 1999). These patterns can change with the
local infrastructure (Chang and Diao,, 2021). Recently, Holmes et al., (2017)
and Bashar, (2021) particularly focus on intra-city house prices and show
significant temporal and spatial dependence in the growth rates. In contrast
to these studies, we focus on the analysis of the log-volatility as a measure
of the market risk.

Figure 1: Indication of spatial, temporal and spatiotemporal volatility
clustering. Left: Average log-squared house price return for each of the 190
zip-code areas over the period from February 1995 to December 2015. The map
shows a weak clustering effect especially for the outer regions (indicating a
positive spatial dependence in the volatility). Center: Temporal ACF depicted
as a series of boxplots showing the estimated temporal autocorrelation of all
190 locations. Right: Spatiotemporal correlation in terms of the slope of a
regression line between the log squared returns and their temporally lagged
neighbors (Moran’s $I$ of first spatiotemporal lag). There is no clear pattern
with varying coefficients around zero, which may indicate a weak
spatiotemporal dependence.

The rest of the paper proceeds in the following way. In Section 2, we state
our model specification and discuss its properties. In Section 3, we provide
the details of the GMM estimation approach for our model, and formally
establish its large sample properties. In Section 4, we investigate the finite
sample properties of our suggested algorithm through an extensive Monte Carlo
study. In Section 5, we provide the details of our empirical application on
Berlin’s house price returns. In Section 6, we offer our concluding comments
with some directions for future studies. Some technical results are relegated
to Appendix.

##  2 Model Specification

The outcome variable $y_{it}$ of region $i$ at time $t$ is modeled according
to

$\displaystyle y_{it}$
$\displaystyle=h_{it}^{1/2}\operatorname{\varepsilon}_{it},$ (1)
$\displaystyle\log h_{it}$
$\displaystyle=\sum_{l=1}^{p}\sum_{j=1}^{n}\rho_{l0}m_{l,ij}\log
y^{2}_{jt}+\gamma_{0}\log
y^{2}_{j,t-1}+\sum_{l=1}^{p}\sum_{j=1}^{n}\delta_{l0}m_{l,ij}\log
y^{2}_{j,t-1}$
$\displaystyle\quad+\mathbf{x}^{{}^{\prime}}_{it}\boldsymbol{\beta}_{0}+\mu_{i0}+\alpha_{t0},$
(2)

for $i=1,2,\ldots,n$ and $t=1,\ldots T$. The spatial locations indexed by
$i=1,\ldots,n$ are supposed to be on discrete regular (e.g., for image
processes) or irregular lattice, also known as spatial polygons. The latter
case is typically present in economics, e.g., regional legal units, districts,
countries, etc. Here, $h_{it}$ is considered as the volatility term in region
$i$ at time $t$, and $\operatorname{\varepsilon}_{it}$ are independent and
identically distributed random variables that has mean zero and unit variance.
The log-volatility terms follow the process in (2), where
$\\{m_{l,ij}\\}_{l=1}^{p}$, for $i,j=1,\ldots,n$, are the non-stochastic
spatial weights. Here, $p$ is a finite positive integers, and
$\\{m_{l,ii}\\}_{l=1}^{p}$ are zero for $i=1,\ldots,n$. The spatial, temporal
and spatiotemporal effects of the log-squared outcome variable on the log-
volatility are measured by the unknown parameters $\gamma_{0}$,
$\\{\rho_{l0}\\}_{l=1}^{p}$, and $\\{\delta_{l0}\\}_{l=1}^{q}$, respectively.
In (2), $\mathbf{x}_{i}$ is a $k\times 1$ vector of exogenous variables with
the associated parameter vector $\boldsymbol{\beta}_{0}$, and the regional and
time fixed effects are denoted by
$\boldsymbol{\mu}_{0}=(\mu_{10},\ldots,\mu_{n0})^{{}^{\prime}}$ and
$\boldsymbol{\alpha}_{0}=(\alpha_{10},\ldots,\alpha_{T0})^{{}^{\prime}}$. Both
$\boldsymbol{\mu}_{0}$ and $\boldsymbol{\alpha}_{0}$ can be correlated with
the exogenous variables in an arbitrary manner. We assume that the initial
value vector $\mathbf{Y}_{0}=(y_{10},\ldots,y_{n0})^{{}^{\prime}}$ is
observable.

Squaring both sides of (1) and then taking the natural logarithm yield

$\displaystyle y^{*}_{it}=h^{*}_{it}+\operatorname{\varepsilon}^{*}_{it},$ (3)

where $y^{*}_{it}=\log y^{2}_{it}$, $h^{*}_{it}=\log h_{it}$ and
$\operatorname{\varepsilon}^{*}_{it}=\log\operatorname{\varepsilon}^{2}_{it}$.
In vector form, we can express (3) and (2) as

$\displaystyle\mathbf{Y}^{*}_{t}=\mathbf{h}^{*}_{t}+\boldsymbol{\operatorname{\varepsilon}}^{*}_{t},$
(4)
$\displaystyle\mathbf{h}^{*}_{t}=\sum_{l=1}^{p}\rho_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t}+\gamma_{0}\mathbf{Y}^{*}_{t-1}+\sum_{l=1}^{p}\delta_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t-1}+\mathbf{X}_{t}\boldsymbol{\beta}_{0}+\boldsymbol{\mu}_{0}+\alpha_{t0}\mathbf{1}_{n},$
(5)

where $\mathbf{M}_{l}=(m_{l,ij})$ is the $n\times n$ spatial weight matrices,
$\mathbf{Y}^{*}_{t}=(y^{*}_{1t},\ldots,y^{*}_{nt})^{{}^{\prime}}$,
$\mathbf{h}^{*}_{t}=(h^{*}_{1t},\ldots,h^{*}_{nt})^{{}^{\prime}}$,
$\boldsymbol{\operatorname{\varepsilon}}^{*}_{t}=(\operatorname{\varepsilon}_{1t},\ldots,\operatorname{\varepsilon}^{*}_{nt})^{{}^{\prime}}$,
$\mathbf{X}_{t}=(\mathbf{x}_{1t},\ldots,\mathbf{x}_{nt})^{{}^{\prime}}$, and
$\mathbf{1}_{n}$ is the $n\times 1$ vector of ones. The process in (5)
indicates that $\mathbf{h}^{*}_{t}$ depends on the high order spatial lags of
$\mathbf{Y}^{*}_{t}$ and $\mathbf{Y}^{*}_{t-1}$. Substituting (5) into (4), we
obtain

$\displaystyle\mathbf{Y}^{*}_{t}=\sum_{l=1}^{p}\rho_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t}+\gamma_{0}\mathbf{Y}^{*}_{t-1}+\sum_{l=1}^{p}\delta_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t-1}+\mathbf{X}_{t}\boldsymbol{\beta}_{0}+\boldsymbol{\mu}_{0}+\alpha_{t0}\mathbf{1}_{n}+\boldsymbol{\operatorname{\varepsilon}}^{*}_{t},$
(6)

for $t=1,\ldots,T$. This transformed model indicates that our specification in
(2) implies a high-order spatial dynamic panel data model for the log-squared
outcome variable. In next section, we show how (6) can be used to estimate the
parameters of the model.

##  3 The Estimation Approach

The elements of $\boldsymbol{\operatorname{\varepsilon}}^{*}_{t}$ in (6) are
i.i.d across $i$ and $t$ but their mean may not be zero. Therefore, we add and
subtract
$\mathbb{E}\left(\boldsymbol{\operatorname{\varepsilon}}^{*}_{t}\right)$ to
obtain the following equation.

$\displaystyle\mathbf{Y}^{*}_{t}$
$\displaystyle=\sum_{l=1}^{p}\rho_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t}+\gamma_{0}\mathbf{Y}^{*}_{t-1}+\sum_{l=1}^{p}\delta_{l0}\mathbf{M}_{l}\mathbf{Y}^{*}_{t-1}+\mathbf{X}_{t}\boldsymbol{\beta}_{0}+\boldsymbol{\mu}_{0}+\alpha_{t0}\mathbf{1}_{n}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}+\mathbf{U}_{t},$
(7)

where
$\mathbf{U}_{t}=(u_{1t},\ldots,u_{nt})^{{}^{\prime}}=\boldsymbol{\operatorname{\varepsilon}}^{*}_{t}-\mathbb{E}\left(\boldsymbol{\operatorname{\varepsilon}}^{*}_{t}\right)$,
and
$\mu_{\operatorname{\varepsilon}}=\mathbb{E}\left(\operatorname{\varepsilon}^{*}_{it}\right)$.
Let $\sigma^{2}_{0}=\mathbb{E}(u^{2}_{it})$. Then, it follows that the
elements of $\mathbf{U}_{t}$ are i.i.d across $i$ and $t$ with mean zero and
variance $\sigma^{2}_{0}$. We need to eliminate both fixed effects terms from
the model in order to avoid the incidental parameter problem. To eliminate
$\boldsymbol{\mu}$ and $\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}$ from
the model, we consider an orthonormal transformation based on the matrix
decomposition of
$\mathbf{J}_{T}=\left(\mathbf{I}_{T}-\frac{1}{T}\mathbf{1}_{T}\mathbf{1}^{{}^{\prime}}_{T}\right)$,
where $\mathbf{I}_{T}$ is the $T\times T$ identity matrix. Let
$\left(\mathbf{F}_{T,T-1},\frac{1}{\sqrt{T}}\mathbf{1}_{T}\right)$ be the
orthonormal eigenvector matrix of $\mathbf{J}_{T}$, where $\mathbf{F}_{T,T-1}$
is the $T\times(T-1)$ sub-matrix containing eigenvectors corresponding to the
eiegenvalues of one. Let
$\mathbf{C}=\left(\mathbf{c}_{1},\ldots,\mathbf{c}_{T}\right)$ be an $n\times
T$ matrix, where $\mathbf{c}_{t}$ is an $n\times 1$ vector for $t=1,\ldots,T$.
Using $\mathbf{F}_{T,T-1}$, we can transform $\mathbf{C}$ into a
$n\times(T-1)$ matrix in the following way:
$\left(\mathbf{c}^{*}_{1},\ldots,\mathbf{c}^{*}_{T-1}\right)=\left(\mathbf{c}_{1},\ldots,\mathbf{c}_{T}\right)\mathbf{F}_{T,T-1}$,
where $\mathbf{c}^{*}_{j}$ is the $j$th column of
$\mathbf{C}\mathbf{F}_{T,T-1}$ for $j=1,\ldots,T-1$. If we apply
$\mathbf{F}_{T,T-1}$ to our model in (7) in a similar manner, we obtain

$\displaystyle\mathbf{Y}^{**}_{t}$
$\displaystyle=\sum_{l=1}^{q}\rho_{l0}\mathbf{M}_{l}\mathbf{Y}^{**}_{t}+\gamma_{0}\mathbf{Y}^{**,-1}_{t-1}+\sum_{l=1}^{q}\delta_{l0}\mathbf{M}_{l}\mathbf{Y}^{**,-1}_{t-1}+\mathbf{X}^{*}_{t}\boldsymbol{\beta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}+\mathbf{U}^{*}_{t},$
(8)

for $t=1,\ldots,T-1$, where
$\left(\mathbf{Y}^{**}_{1},\ldots,\mathbf{Y}^{**}_{T-1}\right)=\left(\mathbf{Y}^{*}_{1},\ldots,\mathbf{Y}^{*}_{T}\right)\mathbf{F}_{T,T-1}$,
$\left(\mathbf{Y}^{**,-1}_{0},\ldots,\mathbf{Y}^{**,-1}_{T-2}\right)=\left(\mathbf{Y}^{*}_{0},\ldots,\mathbf{Y}^{*}_{T-1}\right)\mathbf{F}_{T,T-1}$,
$\left(\mathbf{X}^{*}_{l1},\ldots,\mathbf{X}^{*}_{k,T-1}\right)=\left(\mathbf{X}_{l1},\ldots,\mathbf{X}_{lT}\right)\mathbf{F}_{T,T-1}$,
where $\mathbf{X}_{lt}$ is the $l$th column of $\mathbf{X}_{t}$ for
$l=1,\ldots,k$,
$\left(\alpha^{*}_{1},\ldots,\alpha^{*}_{T-1}\right)=\left(\alpha_{1},\ldots,\alpha_{T}\right)\mathbf{F}_{T,T-1}$
and
$\left(\mathbf{U}^{*}_{1},\ldots,\mathbf{U}^{*}_{T-1}\right)=\left(\mathbf{U}_{1},\ldots,\mathbf{U}_{T}\right)\mathbf{F}_{T,T-1}$.
Note that both $\boldsymbol{\mu}_{0}$ and
$\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}$ are dropped from the model
since

$\displaystyle\left(\boldsymbol{\mu}_{0},\ldots,\boldsymbol{\mu}_{0}\right)\mathbf{F}_{T,T-1}=\boldsymbol{\mu}_{0}\mathbf{1}^{{}^{\prime}}_{T}\mathbf{F}_{T,T-1}=\mathbf{0}_{n\times(T-1)},$
$\displaystyle\left(\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n},\ldots,\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}\right)\mathbf{F}_{T,T-1}=\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}\mathbf{1}^{{}^{\prime}}_{T}\mathbf{F}_{T,T-1}=\mathbf{0}_{n\times(T-1)}.$

Let $N=n(T-1)$ and
$\mathbf{U}_{N}=(\mathbf{U}^{*^{\prime}}_{1},\ldots,\mathbf{U}^{*^{\prime}}_{T-1})^{{}^{\prime}}$.
Note that we can express
$\left(\mathbf{U}^{*}_{1},\ldots,\mathbf{U}^{*}_{T-1}\right)=\left(\mathbf{U}_{1},\ldots,\mathbf{U}_{T}\right)\mathbf{F}_{T,T-1}$
as
$\left(\mathbf{U}^{*^{\prime}}_{1},\ldots,\mathbf{U}^{*^{\prime}}_{T-1}\right)^{{}^{\prime}}=\left(\mathbf{F}^{{}^{\prime}}_{T,T-1}\otimes\mathbf{I}_{n}\right)\left(\mathbf{U}^{{}^{\prime}}_{1},\ldots,\mathbf{U}^{{}^{\prime}}_{T}\right)^{{}^{\prime}}$.333Note
that the matrix equation $\mathbf{A}\mathbf{B}\mathbf{C}=\mathbf{D}$, where
$\mathbf{D}$, $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are suitable
matrices, can be expressed as
$\text{vec}(\mathbf{D})=(\mathbf{C}^{{}^{\prime}}\otimes\mathbf{A})\text{vec}(\mathbf{B})$,
where $\text{vec}(\mathbf{B})$ denotes the vectorization of the matrix
$\mathbf{B}$ (Abadir and Magnus,, 2005, p. 282). This property can be applied
to
$\left(\mathbf{U}^{*}_{1},\mathbf{U}^{*}_{2},\ldots,\mathbf{U}^{*}_{T-1}\right)=\left(\mathbf{U}_{1},\mathbf{U}_{2},\ldots,\mathbf{U}_{T}\right)\mathbf{F}_{T,T-1}$
by setting
$\mathbf{D}=\left(\mathbf{U}^{*}_{1},\mathbf{U}^{*}_{2},\ldots,\mathbf{U}^{*}_{T-1}\right)$,
$\mathbf{C}=\mathbf{F}_{T,T-1}$,
$\mathbf{B}=\left(\mathbf{U}_{1},\mathbf{U}_{2},\ldots,\mathbf{U}_{T}\right)$
and $\mathbf{A}=\mathbf{I}_{n}$.  Then, it follows that

$\displaystyle\mathbb{E}\left(\mathbf{U}_{N}\mathbf{U}^{{}^{\prime}}_{N}\right)=\mathbb{E}\left(\left(\mathbf{F}^{{}^{\prime}}_{T,T-1}\otimes\mathbf{I}_{n}\right)\left(\mathbf{U}^{{}^{\prime}}_{1},\ldots,\mathbf{U}^{{}^{\prime}}_{T}\right)^{{}^{\prime}}\left(\mathbf{U}^{{}^{\prime}}_{1},\ldots,\mathbf{U}^{{}^{\prime}}_{T}\right)\left(\mathbf{F}_{T,T-1}\otimes\mathbf{I}_{n}\right)\right)=\sigma^{2}_{0}\mathbf{I}_{N},$

indicating that the elements of $\mathbf{U}_{N}$ are uncorrelated. Among the
orthonormal transformations, Lee and Yu, (2014) show that the forward
orthogonal difference (the Helmert transformation) can be useful for the
spatial dynamic panel data models. Thus, we have the following explicit forms
for the transformed variables:
$\mathbf{Y}^{**}_{t}=\left(\frac{T-t}{T-t+1}\right)^{1/2}\left(\mathbf{Y}^{*}_{t}-\frac{1}{T-t}\sum_{h=t+1}^{T}\mathbf{Y}^{*}_{h}\right)$,
$\mathbf{Y}^{**,-1}_{t-1}=\left(\frac{T-t}{T-t+1}\right)^{1/2}\left(\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{h=t}^{T-1}\mathbf{Y}^{*}_{h}\right)$
and the other variables are expressed similarly.

The transformed model in (8) includes the transformed time fixed effects.
These terms can be eliminated by pre-multiplying the model with
$\mathbf{J}_{n}=\left(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}_{n}\mathbf{1}^{{}^{\prime}}_{n}\right)$
to get

$\displaystyle\mathbf{J}_{n}\mathbf{Y}^{**}_{t}$
$\displaystyle=\sum_{l=1}^{p}\rho_{l0}\mathbf{J}_{n}\mathbf{M}_{l}\mathbf{Y}^{**}_{t}+\gamma_{0}\mathbf{J}_{n}\mathbf{Y}^{**,-1}_{t-1}+\sum_{l=1}^{p}\delta_{l0}\mathbf{J}_{n}\mathbf{M}_{l}\mathbf{Y}^{**,-1}_{t-1}+\mathbf{J}_{n}\mathbf{X}^{*}_{t}\boldsymbol{\beta}_{0}+\mathbf{J}_{n}\mathbf{U}^{*}_{t},$
(9)

where we used the fact that $\mathbf{J}_{n}\mathbf{1}_{n}=\mathbf{0}_{n}$. Our
GMM estimation approach is based on (9). It is clear that we need to determine
IVs for the following terms:
$\\{\mathbf{M}_{j}\mathbf{Y}^{**}_{t}\\}_{j=1}^{p}$,
$\mathbf{Y}^{**,-1}_{t-1}$ and
$\\{\mathbf{M}_{j}\mathbf{Y}^{**,-1}_{t-1}\\}_{j=1}^{p}$ for $t=1,\ldots,T-1$.
That is, we need IVs for the following variables:

$\displaystyle\mathbf{J}_{n}\left(\mathbb{M}\mathbf{Y}^{**}_{t},\mathbf{Y}^{**,-1}_{t-1},\mathbb{M}\mathbf{Y}^{**,-1}_{t-1}\right),$
(10)

where
$\mathbb{M}\mathbf{Y}^{**}_{t}=\left(\mathbf{M}_{1}\mathbf{Y}^{**}_{t},\ldots,\mathbf{M}_{p}\mathbf{Y}^{**}_{t}\right)$
and
$\mathbb{M}\mathbf{Y}^{**,-1}_{t-1}=\left(\mathbf{M}_{1}\mathbf{Y}^{**,-1}_{t-1},\ldots,\mathbf{M}_{p}\mathbf{Y}^{**,-1}_{t-1}\right)$.
Let $\mathcal{F}_{t-1}$ be the $\sigma$-algebra generated by
$\left(\mathbf{Y}_{0},\ldots,\mathbf{Y}_{t-1}\right)$ conditional on
$\left(\mathbf{X}_{1},\ldots,\mathbf{X}_{T},\boldsymbol{\mu}_{0},\boldsymbol{\alpha}_{0}\right)$.
Then, we can formulate the theoretical linear IVs based on the expectation of
(10) conditional on $\mathcal{F}_{t-1}$.

We use
$\boldsymbol{\rho}_{0}=\left(\rho_{10},\ldots,\rho_{p0}\right)^{{}^{\prime}}$
and
$\boldsymbol{\delta}_{0}=\left(\delta_{10},\ldots,\delta_{p0}\right)^{{}^{\prime}}$
to denote the true parameter values, and
$\boldsymbol{\rho}=\left(\rho_{1},\ldots,\rho_{p}\right)^{{}^{\prime}}$ and
$\boldsymbol{\delta}=\left(\delta_{1},\ldots,\delta_{p}\right)^{{}^{\prime}}$
to denote arbitrary parameter values. Let
$\mathbf{S}(\boldsymbol{\rho})=\left(\mathbf{I}_{n}-\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\right)$,
$\mathbf{S}\equiv\mathbf{S}(\boldsymbol{\rho}_{0})$,
$\mathbf{A}(\boldsymbol{\rho},\boldsymbol{\delta},\gamma)=\mathbf{S}^{-1}(\boldsymbol{\rho})\left(\gamma\mathbf{I}_{n}+\sum_{l=1}^{p}\delta_{l}\mathbf{M}_{l}\right)$,
and
$\mathbf{A}\equiv\mathbf{A}(\boldsymbol{\rho}_{0},\boldsymbol{\delta}_{0},\gamma_{0})$.
Then, the reduced form of (8) can be expressed as

$\displaystyle\mathbf{Y}^{**}_{t}$
$\displaystyle=\mathbf{A}\mathbf{Y}^{**,-1}_{t-1}+\mathbf{S}^{-1}\left(\mathbf{X}^{*}_{t}\boldsymbol{\beta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}+\mathbf{U}^{*}_{t}\right).$
(11)

Let
$\mathbf{Z}^{**}_{t}=\left(\mathbf{Y}^{**,-1}_{t-1},\mathbb{M}\mathbf{Y}^{**,-1}_{t-1},\mathbf{X}^{*}_{t}\right)$
be the $n\times k_{z}$ matrix, where $k_{z}=p+k+1$, and
$\mathbf{Z}_{N}=\left(\mathbf{Z}^{**^{\prime}}_{1},\ldots,\mathbf{Z}^{**^{\prime}}_{T-1}\right)^{{}^{\prime}}$.
Then, using (8), we have

$\displaystyle\mathbf{M}_{r}\mathbf{Y}^{**}_{t}=\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right)+\mathbf{G}_{r}\mathbf{U}^{*}_{t},\quad
r=1,\ldots,p,$ (12)

where
$\boldsymbol{\eta}_{0}=(\gamma_{0},\boldsymbol{\delta}^{{}^{\prime}}_{0},\boldsymbol{\beta}^{{}^{\prime}}_{0})^{{}^{\prime}}$,
and $\mathbf{G}_{r}=\mathbf{M}_{r}\mathbf{S}^{-1}$. We can use (12) to
determine IVs for $\\{\mathbf{M}_{j}\mathbf{Y}^{**}_{t}\\}_{j=1}^{p}$. In the
case of $\mathbf{Y}^{**,-1}_{t-1}$, we can use all strictly exogenous
variables $\mathbf{X}^{*}_{s}$ for $s=1,\ldots,T-1$, and the time lag
variables $\mathbf{Y}^{*}_{0},\ldots,\mathbf{Y}^{*}_{t-1}$ as IVs. Similarly,
we can use $\mathbf{M}_{j}\mathbf{X}^{*}_{s}$ for $s=1,\ldots,T-1$, and
$\mathbf{M}_{j}\mathbf{Y}^{*}_{s}$ for $s=0,1,\ldots,t-1$ as IVs for
$\mathbf{M}_{j}\mathbf{Y}^{**,-1}_{t-1}$. Let $\mathbf{Q}_{t}$ be the $n\times
k_{q}$ matrix of IVs for $t=1,\ldots,T-1$, where $k_{q}\geq k+2p+1$. For
example, we may choose $\mathbf{Q}_{t}$ as

$\displaystyle\left(\mathbf{Y}^{*}_{t-1},\,\mathbb{M}\mathbf{Y}^{*}_{t-1},\mathbb{M}^{2}\mathbf{Y}^{*}_{t-1},\mathbf{X}^{*}_{t},\mathbb{M}\mathbf{X}^{*}_{t},\mathbb{M}^{2}\mathbf{X}^{*}_{t}\right),$
(13)

where
$\mathbb{M}^{2}\mathbf{Y}^{*}_{t-1}=\left(\mathbf{M}^{2}_{1}\mathbf{Y}^{*}_{t-1},\ldots,\mathbf{M}_{1}\mathbf{M}_{p}\mathbf{Y}^{*}_{t-1},\mathbf{M}_{2}\mathbf{M}_{1}\mathbf{Y}^{*}_{t-1},\ldots,\mathbf{M}_{1}\mathbf{M}_{p}\mathbf{Y}^{*}_{t-1},\ldots,\mathbf{M}^{2}_{p}\mathbf{Y}^{*}_{t-1}\right)$,
and $\mathbb{M}^{2}\mathbf{X}^{*}_{t}$ is defined similarly. Denote
$\mathbf{Q}_{N}=(\mathbf{Q}^{{}^{\prime}}_{1},\ldots,\mathbf{Q}^{{}^{\prime}}_{T-1})^{{}^{\prime}}$,
$\mathbf{J}_{N}=\mathbf{I}_{T-1}\otimes\mathbf{J}_{n}$ and
$\mathbf{S}_{N}(\boldsymbol{\rho})=\mathbf{I}_{T-1}\otimes\mathbf{S}(\boldsymbol{\rho})$.
Then, the linear moment conditions based on $\mathbf{Q}_{N}$ can be formulated
as

$\displaystyle\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta}),$
(14)

where
$\boldsymbol{\theta}=\left(\boldsymbol{\rho}^{{}^{\prime}},\boldsymbol{\eta}^{{}^{\prime}}\right)^{{}^{\prime}}$,
$\mathbf{U}_{N}(\boldsymbol{\theta})=\left(\mathbf{U}^{*^{\prime}}_{1}(\boldsymbol{\theta}),\ldots,\mathbf{U}^{*^{\prime}}_{T-1}(\boldsymbol{\theta})\right)^{{}^{\prime}}$
and
$\mathbf{U}^{*}_{t}(\boldsymbol{\theta})=\mathbf{S}(\boldsymbol{\rho})\mathbf{Y}^{**}_{t}-\mathbf{Z}^{**}_{t}\boldsymbol{\eta}-\alpha^{*}_{t}\mathbf{1}_{n}$.
Note that the transformed time fixed effects
$\boldsymbol{\alpha}^{*}=(\alpha^{*}_{1},\ldots,\alpha^{*}_{T-1})^{{}^{\prime}}$
will be eliminated in the moment function because
$\mathbf{U}_{N}(\boldsymbol{\theta})$ is pre-multiplied by $\mathbf{J}_{N}$.

Following Lee, (2007) and Lee and Yu, (2014), we also consider the quadratic
moment functions for estimation. The quadratic moment functions are based on
the idea that the vector $\mathbf{P}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}$ can
be uncorrelated with $\mathbf{J}_{n}\mathbf{U}^{*}_{t}$ for an $n\times n$
matrix $\mathbf{P}_{l}$ satisfying
$\mathrm{tr}\left(\mathbf{J}_{n}\mathbf{P}_{l}\mathbf{J}_{n}\right)=0$, while
it may be correlated with $\mathbf{G}_{r}\mathbf{U}^{*}_{t}$ in (12). Let
$\mathbf{P}_{lN}=\mathbf{I}_{T-1}\otimes\mathbf{P}_{l}$, and assume that there
are $m$ such quadratic moment matrices. Then, the quadratic moment functions
can be expressed as

$\displaystyle\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta}),$
(15)

for $l=1,2,\ldots,m$. Combining the linear and quadratic moment functions, we
obtain the following vector of moment functions,

$\displaystyle\mathbf{g}_{N}(\boldsymbol{\theta})=\begin{pmatrix}\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{1N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\vdots\\\
\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{mN}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\end{pmatrix}.$
(16)

Let $\text{vec}(\mathbf{P})$ be the vectorization of the square matrix
$\mathbf{P}$, $\text{vec}_{D}(\mathbf{P})$ be the column vector formed from
the diagonal elements of $\mathbf{P}$ and
$\mathbf{P}^{s}=\mathbf{P}+\mathbf{P}^{{}^{\prime}}$. Define
$\boldsymbol{\Omega}_{N}=\frac{1}{N}\mathbb{E}\left(\mathbf{g}_{N}(\boldsymbol{\theta}_{0})\mathbf{g}_{N}(\boldsymbol{\theta}_{0})\right)$.
Then, using Lemma 1, it can be shown that444In applying Lemma 1, we use the
fact that
$\mathrm{tr}\left(\mathbf{A}^{{}^{\prime}}\mathbf{B}\right)=\text{vec}^{{}^{\prime}}\left(\mathbf{A}\right)\text{vec}\left(\mathbf{B}\right)=\text{vec}^{{}^{\prime}}\left(\mathbf{B}\right)\text{vec}\left(\mathbf{A}\right)$,
where $\mathbf{A}$ and $\mathbf{B}$ are any two $N\times N$ matrices.

$\displaystyle\boldsymbol{\Omega}_{N}=\operatorname{plim}_{n\to\infty}\frac{\sigma^{4}_{0}}{N}\begin{pmatrix}\boldsymbol{\Delta}_{mN}&\mathbf{0}_{m\times
q}\\\ \mathbf{0}_{q\times
m}&\frac{1}{\sigma^{2}_{0}}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}_{N}\end{pmatrix}+\lim_{n\to\infty}\frac{\mu_{4}-3\sigma^{4}_{0}}{N}\begin{pmatrix}\boldsymbol{\omega}^{{}^{\prime}}_{mN}\boldsymbol{\omega}_{mN}&\mathbf{0}_{m\times
q}\\\ \mathbf{0}_{q\times m}&\mathbf{0}_{q\times q}\end{pmatrix},$ (17)

where $\mu_{4}$ is the fourth moment of $u_{it}$,
$\boldsymbol{\omega}_{mN}=\left(\text{vec}_{D}\left(\mathbf{J}_{N}\mathbf{P}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\right),\ldots,\text{vec}_{D}\left(\mathbf{J}_{N}\mathbf{P}^{{}^{\prime}}_{mN}\mathbf{J}_{N}\right)\right)$
and

$\displaystyle\boldsymbol{\Delta}_{mN}$
$\displaystyle=\left(\text{vec}\left(\mathbf{J}_{N}\mathbf{P}^{{}^{\prime}}_{1N}\right),\ldots,\text{vec}\left(\mathbf{J}_{N}\mathbf{P}^{{}^{\prime}}_{mN}\mathbf{J}_{N}\right)\right)^{{}^{\prime}}$
$\displaystyle\times\left(\text{vec}\left(\mathbf{J}_{N}\mathbf{P}^{s}_{1N}\mathbf{J}_{N}\right),\ldots,\text{vec}\left(\mathbf{J}_{N}\mathbf{P}^{s}_{mN}\mathbf{J}_{N}\right)\right).$

Let $\hat{\boldsymbol{\Omega}}_{N}$ be a consistent estimator of
$\boldsymbol{\Omega}_{N}$, i.e.,
$\hat{\boldsymbol{\Omega}}_{N}-\boldsymbol{\Omega}_{N}=o_{p}(1)$. Then, the
optimal GMM estimator is defined as

$\displaystyle\hat{\boldsymbol{\theta}}_{N}=\operatorname{argmin}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbf{g}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\hat{\boldsymbol{\Omega}}^{-1}_{N}\mathbf{g}_{N}(\boldsymbol{\theta}).$
(18)

To investigate the asymptotic properties of $\hat{\boldsymbol{\theta}}_{N}$,
we require the following assumptions.

######  Assumption 1.

The disturbance terms $u_{it}$, for $i=1,2,\ldots,n$, and $t=1,2,\ldots,T$,
are i.i.d. across $i$ and $t$ with mean zero, variance $\sigma^{2}_{0}$ and
$\mathbb{E}\left(|u_{it}|^{4+\kappa}\right)$ for some $\kappa>0$.

######  Assumption 2.

The spatial weights matrices $\\{\mathbf{M}_{l}\\}_{l=1}^{q}$ are uniformly
bounded in both row and column sums in absolute value.

######  Assumption 3.

(i) $\mathbf{S}(\boldsymbol{\rho})$ is invertible for all
$\boldsymbol{\rho}\in\boldsymbol{\Lambda}$, where $\boldsymbol{\Lambda}$ is a
compact parameter space, and $\boldsymbol{\rho}_{0}$ is in the interior of
$\boldsymbol{\Lambda}$. (ii) $\mathbf{S}^{-1}(\boldsymbol{\rho})$ is uniformly
bounded in both row and column sums in absolute value.

######  Assumption 4.

(i) $\mathbf{X}_{t}$ is non-stochastic with
$\sum_{t=1}^{T}\sum_{i=1}^{n}|x_{it,l}|^{2+\epsilon}<\infty$ for some
$\epsilon>0$, where $x_{it,l}$ is the $(i,t)$th element of the $l$th regressor
for $l=1,\ldots,k$. Moreover,
$\lim_{n\to\infty}\frac{1}{N}\mathbf{X}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{X}_{N}$
exists and is non-singular, where
$\mathbf{X}_{N}=\left(\mathbf{X}^{*^{\prime}}_{1},\ldots,\mathbf{X}^{*^{\prime}}_{T-1}\right)^{{}^{\prime}}$.
(ii) $\boldsymbol{\mu}_{0}$ and $\boldsymbol{\alpha}_{0}$ are non-stochastic
with $\sup_{n}\frac{1}{n}\sum_{i=1}^{n}|\mu_{i0}|^{2+\epsilon}<\infty$ and
$\sup_{T}\frac{1}{T}\sum_{t=1}^{T}|\alpha_{t0}|^{2+\epsilon}<\infty$.

######  Assumption 5.

(i)
$\mathbf{Y}^{*}_{0}=\sum_{h=0}^{\bar{h}}\mathbf{A}^{h}\mathbf{S}^{-1}\left(\mathbf{X}_{-h}\boldsymbol{\beta}_{0}+\boldsymbol{\mu}_{0}+\alpha_{-h0}\mathbf{1}_{n}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}+\mathbf{U}_{-h}\right)$,
where $\bar{h}$ can be finite or infinite. (ii)
$\sum_{h=0}^{\infty}\text{abs}(\mathbf{A}^{h})$ is uniformly bounded in both
row and column sums in absolute value, where the $(i,j)$th element of
$\text{abs}(\mathbf{A})$ is given by $|A_{ij}|$ and $A_{ij}$ is the $(i,j)$th
element of $\mathbf{A}$.

######  Assumption 6.

$\mathbb{E}\left(\mathbf{Q}_{t}|\mathcal{F}_{t-1}\right)=\mathbf{Q}_{t}$ and
$\mathbb{E}\left(|q_{it,l}|^{2+\epsilon}\right)<\infty$, where $q_{it,l}$ is
the $(i,t)$th element of the $l$th column of $\mathbf{Q}_{t}$. Moreover,
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{Z}_{N},\,\mathbf{L}_{N}\right)$
and
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}_{N}$
have full column ranks, where
$\mathbf{L}_{N}=\left(\mathbf{L}^{{}^{\prime}}_{1},\ldots,\mathbf{L}^{{}^{\prime}}_{T-1}\right)^{{}^{\prime}}$
with $\mathbf{L}_{t}=\left(\mathbf{L}_{1,t},\ldots,\mathbf{L}_{p,t}\right)$
and
$\mathbf{L}_{r,t}=\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}\right)$
for $r=1,2,\ldots,p$.

Assumption 1 specifies the distribution of the elements of $\mathbf{U}_{t}$
for $t=1,2,\ldots,T$. The moment condition in this assumption is required for
showing the asymptotic distribution of our set of moment functions.
Assumptions 2 and 3 are standard assumptions adopted in the literature for
limiting the degree of spatial correlation at a manageable degree, e.g., among
others, see Kelejian and Prucha, (2010); Lee, (2004). Assumption 4 provides
the regularity conditions for $\mathbf{X}_{N}$, $\boldsymbol{\mu}_{0}$ and
$\boldsymbol{\alpha}_{0}$. The first part of Assumption 5 specifies
$\mathbf{Y}^{*}_{0}$, and the remaining parts are required to limit dependence
over time and across cross section units (see Lee and Yu, (2014) for the
details). The sufficient conditions for the first part of Assumption 3, and
the second part of 5 can be determined. Let $\|\cdot\|$ be any matrix norm.
Then, the following respective conditions will be sufficient for ensuring
these parts: (i) $\|\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\|<1$ and (ii)
$\|\mathbf{A}(\boldsymbol{\rho},\boldsymbol{\delta},\gamma)\|<1$. Note that

$\displaystyle\left\|\sum_{j=1}^{p}\rho_{l}\mathbf{M}_{l}\right\|\leq|\rho_{1}|\cdot\|\mathbf{M}_{1}\|+\ldots+|\rho_{p}|\cdot\|\mathbf{M}_{p}\|\leq\left(\sum_{j=1}^{p}|\rho_{j}|\right)\times\max_{1\leq
j\leq p}\|\mathbf{M}_{{}_{j}}\|,$

Thus, a relatively restrictive condition for (i) is
$\left(\sum_{j=1}^{p}|\rho_{j}|\right)\times\max_{1\leq j\leq
p}\|\mathbf{M}_{{}_{j}}\|<1$. Similarly, we have

$\displaystyle\left\|\mathbf{A}(\boldsymbol{\rho},\boldsymbol{\delta},\gamma)\right\|\leq\left\|\mathbf{S}^{-1}(\boldsymbol{\rho})\right\|\times\left\|\gamma\mathbf{I}_{n}+\sum_{l=1}^{p}\delta_{l}\mathbf{M}_{l}\right\|$
$\displaystyle=\left\|\mathbf{I}_{n}+\left(\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\right)+\left(\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\right)^{2}+\ldots\right\|\times\left\|\gamma\mathbf{I}_{n}+\sum_{l=1}^{p}\delta_{l}\mathbf{M}_{l}\right\|$
$\displaystyle\leq\frac{1}{1-\tau_{1}}\times\left(|\gamma|+\left(\sum_{l=1}^{p}|\delta_{l}|\right)\times\max_{1\leq
l\leq p}\|\mathbf{M}_{l}\|\right),$

where $\tau_{1}=\left(\sum_{l=1}^{p}|\rho_{l}|\right)\times\max_{1\leq l\leq
p}\left\|\mathbf{M}_{l}\right\|<1$ is guaranteed by the first condition. This
result suggests that a relatively restrictive condition for (ii) is
$\frac{1}{1-\tau_{1}}\times\left(|\gamma|+\left(\sum_{l=1}^{p}|\delta_{l}|\right)\times\max_{1\leq
l\leq p}\|\mathbf{M}_{l}\|\right)<1$. When the spatial weights matrices are
row normalized these relatively restrictive conditions can be further
simplified. For example, if we use the matrix row sum norm, we will get the
following sufficient conditions: (i) $\left(\sum_{j=1}^{p}|\rho_{j}|\right)<1$
and (ii)
$\left(\sum_{j=1}^{p}|\rho_{j}|+|\gamma|+\sum_{l=1}^{p}|\delta_{l}|\right)<1$.

Assumption 6 provides the regularity conditions for the IV matrix
$\mathbf{Q}_{t}$. The first part states that $\mathbf{Q}_{t}$ is pre-
determined in the sense that
$\mathbb{E}\left(\mathbf{Q}_{t}|\mathcal{F}_{t-1}\right)=\mathbf{Q}_{t}$. The
moment condition in this assumption is required for the application of a CLT
to the set of our moment functions (see the CLT given in Lemma 3). The full
column rank condition in Assumption 6 gives the identification condition based
on the linear moment function in our setting. See Appendix B for the details
on the identification condition in our setting.

Let
$\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\left(\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\rho}^{{}^{\prime}}},\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}\right)$.
In Section C of Appendix, we show that
$\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\mathbf{D}_{1N}+\mathbf{D}_{2N}+O_{p}(N^{-1/2})$,
where $\mathbf{D}_{1N}=O(1)$ and $\mathbf{D}_{2N}=O(T^{-1})$.555The explicit
forms of $\mathbf{D}_{1N}$ and $\mathbf{D}_{2N}$ are given in Section C of
Appendix. The following result gives the limiting distribution of
$\hat{\boldsymbol{\theta}}_{N}$ under the large $T$ and finite $T$ cases.

######  Theorem 1.

Under Assumptions 1-6, we have the following results,

  1. 1.

When $T$ is finite and $n\to\infty$, we have

$\displaystyle\sqrt{n}\left(\hat{\boldsymbol{\theta}}_{N}-\boldsymbol{\theta}_{0}\right)\xrightarrow{d}N\left(\mathbf{0}_{k_{z}+p},\,\operatorname{plim}_{n\to\infty}\frac{1}{T-1}\left(\left(\mathbf{D}_{1N}+\mathbf{D}_{2N}\right)^{{}^{\prime}}_{N}\boldsymbol{\Omega}^{-1}_{N}\left(\mathbf{D}_{1N}+\mathbf{D}_{2N}\right)\right)^{-1}\right).$
(19)

  2. 2.

When $T\to\infty$ and $n\to\infty$, we have

$\displaystyle\sqrt{N}\left(\hat{\boldsymbol{\theta}}_{N}-\boldsymbol{\theta}_{0}\right)\xrightarrow{d}N\left(\mathbf{0}_{k_{z}+p},\,\operatorname{plim}_{n,T\to\infty}\left(\mathbf{D}^{{}^{\prime}}_{1N}\boldsymbol{\Omega}^{-1}_{N}\mathbf{D}_{1N}\right)^{-1}\right).$
(20)

###### Proof.

See Section C of Appendix. ∎

Our estimator defined in (18) requires a consistent estimator of $\Omega_{N}$.
We can use a plug-in estimator of $\Omega_{N}$ based on an initial GMM
estimator, or alternatively, we can formulate a 2SLS estimator based on
$\mathbf{Q}_{N}$. Let
$\mathbf{Y}_{N}=\left(\mathbf{Y}^{**^{\prime}}_{1},\ldots,\mathbf{Y}^{**^{\prime}}_{T-1}\right)^{{}^{\prime}}$
and
$\mathbf{Z}_{N}=\left(\mathbf{Z}^{**^{\prime}}_{1},\ldots,\mathbf{Z}^{**^{\prime}}_{T-1}\right)^{{}^{\prime}}$.
Then, the 2SLS estimator is

$\displaystyle\tilde{\boldsymbol{\theta}}_{N}=\left(\left(\mathbb{M}_{N}\mathbf{Y}_{N},\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\left(\mathbb{M}_{N}\mathbf{Y}_{N},\mathbf{Z}_{N}\right)\right)^{-1}\left(\mathbb{M}_{N}\mathbf{Y}_{N},\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\mathbf{Y}_{N},$
(21)

where
$\mathbf{M}_{\mathbf{Q}}=\mathbf{J}_{N}\mathbf{Q}_{N}\left(\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}_{N}\right)^{-1}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}$
and
$\mathbb{M}_{N}\mathbf{Y}_{N}=\left(\mathbf{M}_{1N}\mathbf{Y}_{N},\ldots,\mathbf{M}_{pN}\mathbf{Y}_{N}\right)$
with $\mathbf{M}_{jN}=\mathbf{I}_{T-1}\otimes\mathbf{M}_{j}$ for
$j=1,2,\ldots,p$. Our Theorem 1 suggests that

$\displaystyle\sqrt{n}\left(\tilde{\boldsymbol{\theta}}_{N}-\boldsymbol{\theta}_{0}\right)\xrightarrow{d}N\left(\mathbf{0},\,\sigma^{2}_{0}\operatorname{plim}_{n\to\infty}\left(\frac{1}{N}\left(\mathbf{L}_{N},\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\left(\mathbf{L}_{N},\mathbf{Z}_{N}\right)\right)^{-1}\right).$
(22)

In the case of both the initial GMM and 2SLS estimators, the linear IV matrix
can be
$\mathbf{Q}_{t}=\left(\mathbf{Y}^{*}_{t-1},\mathbb{M}\mathbf{Y}^{*}_{t-1},\mathbb{M}^{2}\mathbf{Y}^{*}_{t-1},\mathbf{X}^{*}_{t},\mathbb{M}\mathbf{X}^{*}_{t},\mathbb{M}^{2}\mathbf{X}^{*}_{t}\right)$
for $t=1,2,\ldots,T-1$. We may consider the following quadratic moment
matrices for the initial GMM estimator:
$\mathbf{P}_{j}=\left(\mathbf{M}_{j}-\frac{\mathrm{tr}(\mathbf{M}_{j}\mathbf{J}_{n})}{n-1}\mathbf{J}_{n}\right)$
and
$\mathbf{P}_{j+p}=\left(\mathbf{M}^{2}_{j}-\frac{\mathrm{tr}(\mathbf{M}^{2}_{j}\mathbf{J}_{n})}{n-1}\mathbf{J}_{n}\right)$
for $j=1,2,\ldots,p$. Then, the initial GMM estimator is given by
$\tilde{\boldsymbol{\theta}}_{N}=\operatorname{argmin}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbf{g}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{g}_{N}(\boldsymbol{\theta})$
where

$\displaystyle\mathbf{g}_{N}(\boldsymbol{\theta})=\begin{pmatrix}\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{1N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\vdots\\\
\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{2pN}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\end{pmatrix}.$
(23)

We can use $\tilde{\boldsymbol{\theta}}_{N}$ to formulate the plug-in
estimator of $\Omega_{N}$, which requires the estimators of $\sigma^{2}_{0}$
and $\mu_{4}$. Let
$\tilde{\mathbf{V}}_{t}=\mathbf{S}(\tilde{\boldsymbol{\rho}}_{N})\mathbf{Y}^{**}_{t}-\mathbf{Z}^{**}_{t}\tilde{\boldsymbol{\eta}}_{N}$.
Then, we can estimate $\sigma^{2}_{0}$ by
$\tilde{\sigma}^{2}_{N}=\frac{1}{N}\sum_{t=1}^{T-1}\tilde{\mathbf{V}}^{{}^{\prime}}_{t}\mathbf{J}_{n}\tilde{\mathbf{V}}_{t}$.
Let
$\Delta\tilde{\mathbf{V}}_{t}=\mathbf{S}(\tilde{\boldsymbol{\rho}}_{N})\Delta\mathbf{Y}^{**}_{t}-\Delta\mathbf{Z}^{**}_{t}\tilde{\boldsymbol{\eta}}_{N}$.
Then, following Lee and Yu, (2014), we can estimate $\mu_{4}$ by
$\tilde{\mu}_{4}=\frac{1}{2N}\sum_{i=1}^{n}\sum_{t=2}^{T}\left(\left[\mathbf{J}_{n}\Delta\tilde{\mathbf{V}}_{t}\right]_{i}\right)^{4}-3\tilde{\sigma}^{4}$,
where $\left[\mathbf{J}_{n}\Delta\tilde{\mathbf{V}}_{t}\right]_{i}$ is the
$i$th element of $\mathbf{J}_{n}\Delta\tilde{\mathbf{V}}_{t}$.

Our set of moment functions in (15) depends on the IV matrix $\mathbf{Q}_{N}$
and the quadratic moment matrices $\mathbf{P}_{lN}$ for $l=1,2,\ldots,m$. The
asymptotic efficiency of $\hat{\boldsymbol{\theta}}_{N}$ should be considered
in choosing the IV and quadratic moment matrices. The best set of IV and
quadratic moment matrices is the set that leads to the most efficient GMM
estimator. When $T$ is large, the precision matrix
$\hat{\boldsymbol{\theta}}_{N}$ takes a simple form allowing for determining
the best set of IV and quadratic moment matrices. When $T$ is large, the proof
of Theorem 1 indicates that
$\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\mathbf{D}_{1N}+O_{p}(N^{-1/2})$,
where

$\displaystyle\mathbf{D}_{1N}=-\frac{1}{N}\begin{pmatrix}\sigma^{2}_{0}\mathbf{C}_{N}&\mathbf{0}_{m\times
k_{z}}\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{L}_{N}&\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Z}_{N}\end{pmatrix}.$
(24)

Then, (20) shows that the precision matrix of
$\sqrt{N}\left(\hat{\boldsymbol{\theta}}_{N}-\boldsymbol{\theta}_{0}\right)$
is

$\displaystyle\mathbf{D}^{{}^{\prime}}_{1N}\boldsymbol{\Omega}^{-1}_{N}\mathbf{D}_{1N}$
$\displaystyle=\frac{1}{N}\begin{pmatrix}\mathbf{C}^{{}^{\prime}}_{N}\left(\boldsymbol{\Delta}_{mN}+\frac{\mu_{4}-3\sigma^{4}_{0}}{\sigma^{4}_{0}}\boldsymbol{\omega}^{{}^{\prime}}_{mN}\boldsymbol{\omega}_{mN}\right)^{-1}\mathbf{C}_{N}&\mathbf{0}_{p\times
k_{z}}\\\ \mathbf{0}_{k_{z}\times p}&\mathbf{0}_{k_{z}\times
k_{z}}\end{pmatrix}$
$\displaystyle+\frac{1}{N\sigma^{2}_{0}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right).$
(25)

Since the above precision matrix has the same form as the one given in Lee and
Yu, (2014), we use their approach to determine the best set of quadratic
moment matrices. We will choose the best quadratic matrices by maximizing
$\mathbf{C}^{{}^{\prime}}_{N}\left(\boldsymbol{\Delta}_{mN}+\frac{\mu_{4}-3\sigma^{4}_{0}}{\sigma^{4}_{0}}\boldsymbol{\omega}^{{}^{\prime}}_{mN}\boldsymbol{\omega}_{mN}\right)^{-1}\mathbf{C}_{N}$.
As shown in Lee and Yu, (2014), these matrices are

$\displaystyle\mathbf{P}^{*}_{j}=\left(\mathbf{G}_{j}-\frac{\mathrm{tr}\left(\mathbf{G}_{j}\mathbf{J}_{n}\right)}{n-1}\mathbf{J}_{n}\right)+c\left(\operatorname{Diag}\left(\mathbf{J}_{n}\mathbf{G}_{j}\mathbf{J}_{n}\right)-\frac{\mathrm{tr}\left(\mathbf{G}_{j}\mathbf{J}_{n}\right)}{n}\mathbf{I}_{n}\right),\quad
j=1,2,\ldots,p,$ (26)

where
$c=\left(\frac{n}{n-2}\right)^{2}\left(\frac{1}{n/(n-2)+(\eta_{4}-3)/2}-\frac{n-2}{n}\right)$
and $\eta_{4}=\mu_{4}/\sigma^{4}_{0}$.

In the case of the best linear moment function, we should consider the
conditional mean
$\mathbb{E}\left(\mathbb{M}\mathbf{Y}^{**}_{t},\mathbf{Z}^{**}_{t}|\mathcal{F}_{t-1}\right)$.
The conditional mean of $\mathbf{Y}^{**}_{t-1}$ can be determined from
$\mathbf{Y}^{**,-1}_{t-1}=c_{t}\left(\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{s=t}^{T-1}\mathbf{Y}^{*}_{s}\right)$,
where $c_{t}=\left(\frac{T-t}{T-t+1}\right)^{1/2}$. Using Lemma 4,
$\mathbb{E}\left(\mathbf{Y}^{**}_{t-1}|\mathcal{F}_{t-1}\right)$ can be
approximated by

$\displaystyle\mathbf{H}_{t}$
$\displaystyle=c_{t}\left(\left(\mathbf{I}_{n}-\frac{1}{T-t}\sum_{h=1}^{T-t}\mathbf{A}^{h}\right)\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\left(\mathbf{X}_{r}\boldsymbol{\beta}_{0}+\alpha_{r,0}\mathbf{1}_{n}\right)\right)$
$\displaystyle-
c_{t}\frac{1}{(T-t)(t-1)}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\sum_{s=1}^{t-1}\left(\mathbf{S}\mathbf{Y}^{*}_{s}-\mathbf{Z}^{*}_{s}\boldsymbol{\eta}_{0}-\alpha_{s0}\mathbf{1}_{n}\right),$

where
$\mathbf{Z}^{*}_{s}=\left(\mathbf{Y}^{*}_{s-1},\mathbb{M}\mathbf{Y}^{*}_{s-1},\mathbf{X}_{s}\right)$.
Thus, the best theoretical IV
$\mathbf{J}_{n}\mathbb{E}\left(\mathbf{Y}^{**}_{t-1}|\mathcal{F}_{t-1}\right)$
can be approximated by $\mathbf{J}_{n}\mathbf{H}_{t}$.666Note that when $t=1$,
we may simply use
$\mathbf{H}_{1}=c_{1}\left(\left(\mathbf{I}_{n}-\frac{1}{T-1}\sum_{h=1}^{T-1}\mathbf{A}^{h}\right)\mathbf{Y}^{*}_{0}-\frac{1}{T-1}\sum_{r=1}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\left(\mathbf{X}_{r}\boldsymbol{\beta}_{0}+\alpha_{r,0}\mathbf{1}_{n}\right)\right)$.
Similarly, the best IVs for $\mathbf{J}\mathbf{Z}^{**}_{t}$ can be taken as
$\mathbf{J}_{n}\mathbf{K}_{t}$, where
$\mathbf{K}_{t}=\left(\mathbf{H}_{t},\,\mathbb{M}\mathbf{H}_{t},\,\mathbf{X}^{*}_{t}\right)$.
Using (12), the best IV for $\mathbf{J}_{n}\mathbf{M}_{r}\mathbf{Y}^{**}_{t}$
is
$\mathbf{J}_{n}\mathbf{G}_{r}\left(\mathbf{K}_{t}\boldsymbol{\delta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right)$
for $r=1,2,\ldots,p$. Overall, we may use $\mathbf{J}_{n}\mathbf{Q}^{*}_{t}$
as the IV matrix for
$\mathbf{J}_{n}\left(\mathbb{M}\mathbf{Y}^{**}_{t},\,\mathbf{Z}^{**}_{t}\right)$,
where

$\displaystyle\mathbf{Q}^{*}_{t}=\left(\mathbf{G}_{1}\left(\mathbf{K}_{t}\boldsymbol{\delta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right),\ldots,\mathbf{G}_{p}\left(\mathbf{K}_{t}\boldsymbol{\delta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right),\,\mathbf{K}_{t}\right),\quad
t=1,2,\ldots,T-1.$ (27)

The feasible version of $\mathbf{Q}^{*}_{t}$ can be obtained by substituting
consistent estimators of the unknown parameters into (27). Note that if
$\mathbf{M}_{r}\mathbf{1}_{n}=\mathbf{1}_{n}$, i.e., when $\mathbf{M}_{r}$ is
row-normalized, we have
$\mathbf{J}_{n}\mathbf{M}_{r}=\mathbf{J}_{n}\mathbf{M}_{r}\mathbf{J}_{n}$.
This property suggests that the time fixed effects will dropped from
$\mathbf{J}_{n}\mathbf{Q}^{*}_{t}$ since
$\mathbf{J}_{n}\mathbf{1}_{n}=\mathbf{0}_{n}$. However, if $\mathbf{M}_{r}$’s
are not row normalized, then we also need an estimator of the time fixed
effects to get a feasible version of $\mathbf{Q}^{*}_{t}$. Let
$\hat{\boldsymbol{\vartheta}}_{t}=\mathbf{S}(\hat{\boldsymbol{\lambda}}_{n})\mathbf{Y}^{*}_{t}-\mathbf{Z}^{*}_{t}\hat{\boldsymbol{\eta}}_{n}$
be an estimator of
$\boldsymbol{\mu}_{0}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}+\alpha_{t0}\mathbf{1}_{n}$.
Under the normalization assumption of the form
$\mathbf{1}^{{}^{\prime}}_{n}\left(\boldsymbol{\mu}_{0}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}\right)=\mathbf{0}_{n}$,
we can estimate the time fixed effects by
$\hat{\alpha}_{t}=\frac{1}{n}\mathbf{1}^{{}^{\prime}}_{n}\hat{\boldsymbol{\vartheta}}_{t}$
for $t=1,2,\ldots,T$.777Note that when $T$ is large,
$\tilde{\boldsymbol{\mu}}_{0}=\left(\boldsymbol{\mu}_{0}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}\right)$
can be estimated by
$\hat{\tilde{\boldsymbol{\mu}}}_{n}=\frac{1}{T}\sum_{t=1}^{T}\left(\hat{\boldsymbol{\vartheta}}_{t}-\frac{1}{n}\mathbf{1}^{{}^{\prime}}_{n}\hat{\boldsymbol{\vartheta}}_{t}\mathbf{1}_{n}\right)$.
The following theorem provides our result on the best GMM estimator formulated
with the feasible versions of $\mathbf{Q}^{*}_{t}$ and $\mathbf{P}^{*}_{j}$
for $j=1,2,\ldots,p$.

######  Theorem 2.

Let $\hat{\mathbf{Q}}_{t}$ be the feasible version of $\mathbf{Q}^{*}_{t}$ for
$t=1,2,\ldots,T-1$, and $\hat{\mathbf{P}}^{*}_{j}$ be the feasible version of
$\mathbf{P}^{*}_{j}$ for $j=1,2,\ldots,p$. Consider the set of moment
functions $\mathbf{g}_{N}(\boldsymbol{\theta})$ formulated with
$\hat{\mathbf{Q}}_{t}$ and $\hat{\mathbf{P}}^{*}_{j}$. Then, the feasible best
GMM estimator defined by
$\hat{\boldsymbol{\theta}}^{*}_{N}=\operatorname{argmin}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbf{g}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\hat{\boldsymbol{\Omega}}^{-1}_{N}\mathbf{g}_{N}(\boldsymbol{\theta})$
has the following asymptotic distribution

$\displaystyle\sqrt{N}\left(\hat{\boldsymbol{\theta}}^{*}_{N}-\boldsymbol{\theta}_{0}\right)\xrightarrow{d}N\left(\mathbf{0}_{k_{z}+p},\,\boldsymbol{\Sigma}^{*-1}_{N}\right),$
(28)

where

$\displaystyle\boldsymbol{\Sigma}^{*}_{N}=\lim_{n,T\to\infty}\begin{pmatrix}\mathbf{C}^{*}_{N}/N&\mathbf{0}_{p\times
k_{z}}\\\ \mathbf{0}_{k_{z}\times p}&\mathbf{0}_{k_{z}\times
k_{z}}\end{pmatrix}+\operatorname{plim}_{n,T\to\infty}\frac{1}{N\sigma^{2}_{0}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right),$
(29)

with

$\displaystyle\mathbf{C}^{*}_{N}=\begin{pmatrix}\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{*s}_{1N}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{*s}_{1N}\mathbf{J}_{N}\right)\\\
\vdots&\ddots&\vdots\\\
\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{*s}_{pN}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{*s}_{pN}\mathbf{J}_{N}\right)\end{pmatrix}.$
(30)

###### Proof.

See Section D of Appendix. ∎

##  4 A Monte Carlo Study

In this section, we investigate the finite sample properties of the best GMM
estimator provided in Theorem 2. To that end, we consider
$y_{it}=h_{it}^{1/2}\operatorname{\varepsilon}_{it}$, and the following cases
for $h_{it}$:

$\displaystyle M_{1}:\,\log h_{it}=\sum_{j=1}^{n}\rho_{0}m_{ij}\log
y^{2}_{jt}+\gamma_{0}\log y^{2}_{j,t-1}+\sum_{j=1}^{n}\delta_{0}m_{ij}\log
y^{2}_{j,t-1}+\mathbf{x}^{{}^{\prime}}_{it}\boldsymbol{\beta}_{0}+\mu_{i0}+\alpha_{t0}$
$\displaystyle M_{2}:\,\log h_{it}=\sum_{j=1}^{n}\rho_{0}m_{ij}\log
y^{2}_{jt}+\gamma_{0}\log y^{2}_{j,t-1}+\sum_{j=1}^{n}\delta_{0}m_{ij}\log
y^{2}_{j,t-1}+\mathbf{x}^{{}^{\prime}}_{it}\boldsymbol{\beta}_{0}+\mu_{i0},$
$\displaystyle M_{3}:\,\log
h_{it}=\sum_{l=1}^{2}\sum_{j=1}^{n}\rho_{l0}m_{l,ij}\log
y^{2}_{jt}+\gamma_{0}\log
y^{2}_{j,t-1}+\sum_{l=1}^{2}\sum_{j=1}^{n}\delta_{l0}m_{l,ij}\log
y^{2}_{j,t-1}+\mathbf{x}^{{}^{\prime}}_{it}\boldsymbol{\beta}_{0}+\mu_{i0}+\alpha_{t0},$

where $\mu_{i0}$’s and $\alpha_{t0}$’s are i.i.d $N(0,1)$, and
$\mathbf{x}_{it}\sim$ i.i.d $N(\mathbf{0}_{2\times 1},\mathbf{I}_{2})$ with
$\boldsymbol{\beta}_{0}=(0.5,1)^{{}^{\prime}}$. For the first two models,
denoted $M_{1}$ and $M_{2}$, we consider two different temporal dependence
structures – a weakly temporal dependent model and a strongly persistent
model. Specifically, we set
$(\rho_{0},\gamma_{0},\delta_{0})^{{}^{\prime}}=\\{(0.2,0.2,-0.2)^{{}^{\prime}},\,(0.2,0.8,-0.2)^{{}^{\prime}}\\}$
in $M_{1}$ and $M_{2}$, respectively. Moreover, $M_{2}$ considers the case
without temporal fixed effects, i.e., $\alpha_{t0}=0$ for all $t$. In $M_{3}$,
including higher-order spatial lags, we set
$(\rho_{10},\rho_{20},\gamma_{0},\delta_{10},\delta_{20})=(0.6,0.2,0.1,0.01,0.01)^{{}^{\prime}}$.
That is, we have very weak temporal and spatiotemporal effects. In all cases,
we consider row-normalized queen contiguity spatial weights matrices, where
$\mathbf{M}_{1}$ has positive weights for the first-lag neighbors and
$\mathbf{M}_{2}$ for the second-lag neighbors. Furthermore, we consider two
distributions to generate the disturbance terms: (i)
$\operatorname{\varepsilon}_{it}\sim$ i.i.d $N(0,1)$ and (ii)
$\operatorname{\varepsilon}_{it}\sim$ i.i.d $t_{3}$, where $t_{3}$ is the
Student’s $t$ distribution with $3$ degrees of freedom. We set
$(n,T)=\\{(64,20),(100,40)\\}$, and the number of repetitions to $1000$ in all
cases. Thus, we considered 12 different model specifications in total.

The results of our Monte Carlo simulation study are reported in Tables 1 \- 3.
To evaluate the estimation performance, we report the average bias across all
replications and the mean absolute errors (MAE). For all simulation settings,
our theoretical findings are supported in the finite sample case. More
precisely, when $n$ and $T$ increase, our suggested GMM estimator reports
smaller bias and MAE in all cases. Comparing the performance with respect to
the error distribution, we see slightly lower MAEs in the heavy-tailed case.
These differences are insignificant in almost all cases ($\alpha=0.05$).
Overall, these results indicate that our suggested GMM estimator has good
finite sample properties in terms of bias and MAE.

Table 1: Average bias and mean absolute errors (MAE) of the estimated parameters of Model $M_{1}$. |  |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $N(0,1)$ |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $t_{3}$
---|---|---|---
|  | $n=64$ | $n=100$ | $n=64$ | $n=100$
|  | $T=20$ | $T=40$ | $T=20$ | $T=40$
| $\rho_{0}=0.2$ | 0.0034 | 0.0041 | -0.0001 | 0.0027
| $\gamma_{0}=0.2$ | 0.0001 | -0.0008 | -0.0004 | 0.0005
Bias | $\delta_{0}=-0.2$ | -0.0013 | -0.0004 | 0.0007 | -0.0009
| $\beta_{00}=0.5$ | -0.0033 | -0.0023 | -0.0020 | -0.0003
| $\beta_{10}=1$ | -0.0056 | -0.0007 | -0.0026 | -0.0035
| $\rho_{0}=0.2$ | 0.1142 | 0.0590 | 0.1183 | 0.0653
| $\gamma_{0}=0.2$ | 0.0266 | 0.0139 | 0.0253 | 0.0134
MAE | $\delta_{0}=-0.2$ | 0.0612 | 0.0321 | 0.0606 | 0.0348
| $\beta_{00}=0.5$ | 0.0527 | 0.0292 | 0.0581 | 0.0324
| $\beta_{10}=1$ | 0.0517 | 0.0287 | 0.0573 | 0.0318
Table 2: Average bias and mean absolute errors (MAE) of the estimated parameters of Model $M_{2}$. |  |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $N(0,1)$ |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $t_{3}$
---|---|---|---
|  | $n=64$ | $n=100$ | $n=64$ | $n=100$
|  | $T=20$ | $T=40$ | $T=20$ | $T=40$
| $\rho_{0}=0.2$ | 0.0214 | 0.0119 | 0.0204 | 0.0035
| $\gamma_{0}=0.8$ | -0.0014 | -0.0012 | -0.0025 | -0.0013
Bias | $\delta_{0}=-0.2$ | 0.0109 | -0.0019 | -0.0036 | 0.0033
| $\beta_{00}=0.5$ | -0.0031 | 0.0001 | -0.0039 | -0.0014
| $\beta_{10}=1$ | -0.0054 | -0.0018 | -0.0091 | 0.0002
| $\rho_{0}=0.2$ | 0.1096 | 0.0582 | 0.1161 | 0.0623
| $\gamma_{0}=0.8$ | 0.0372 | 0.0167 | 0.0383 | 0.0157
MAE | $\delta_{0}=-0.2$ | 0.1184 | 0.0627 | 0.1238 | 0.0643
| $\beta_{00}=0.5$ | 0.0526 | 0.0296 | 0.0553 | 0.0321
| $\beta_{10}=1$ | 0.0532 | 0.0276 | 0.0581 | 0.0308
Table 3: Average bias and mean absolute errors (MAE) of the estimated parameters of Model $M_{3}$. |  |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $N(0,1)$ |  $\operatorname{\varepsilon}_{it}\sim$ i.i.d $t_{3}$
---|---|---|---
|  | $n=49$ | $n=100$ | $n=49$ | $n=100$
|  | $T=20$ | $T=40$ | $T=20$ | $T=40$
| $\rho_{10}=0.6$ | 0.0112 | 0.0026 | 0.0150 | 0.0037
| $\rho_{20}=0.2$ | 0.0122 | 0.0041 | 0.0134 | 0.0039
| $\gamma_{0}=0.1$ | -0.0007 | -0.0004 | -0.0008 | 0.0003
Bias | $\delta_{10}=0.01$ | -0.0066 | -0.0014 | -0.0076 | -0.0012
| $\delta_{20}=0.01$ | -0.0059 | -0.0023 | -0.0070 | -0.0037
| $\beta_{00}=0.5$ | -0.0006 | 0.0022 | -0.0043 | -0.0013
| $\beta_{10}=1$ | -0.0064 | -0.0046 | -0.0059 | -0.0062
| $\rho_{10}=0.6$ | 0.0794 | 0.0458 | 0.0864 | 0.0479
| $\rho_{20}=0.2$ | 0.1183 | 0.0650 | 0.1271 | 0.0659
| $\gamma_{0}=0.1$ | 0.0254 | 0.0136 | 0.0261 | 0.0135
MAE | $\delta_{10}=0.01$ | 0.0514 | 0.0277 | 0.0525 | 0.0279
| $\delta_{10}=0.01$ | 0.0646 | 0.0338 | 0.0648 | 0.0351
| $\beta_{00}=0.5$ | 0.0559 | 0.0293 | 0.0556 | 0.0313
| $\beta_{10}=1$ | 0.0534 | 0.0306 | 0.0599 | 0.0337

##  5 Real-World Example: Intra-city housing market risk

The real-estate market is undoubtedly a financial market with the most
apparent spatial and temporal dependence. The location of a property, along
with size and condition, is an important price-determining influence. Hence,
there are pronounced spatial spillover effects in real-estate prices, in
addition to the natural temporal dependence. Furthermore, taxes may
significantly affect the market, such as property taxes or real-estate
transfer taxes. In general, however, taxes appear to play a subordinate role
in purchasing decisions – with one exception, namely, if the transfer taxes
change, some sales could be shifted for a certain period. If, for example, the
land transfer tax increases by one percentage point and one wants to buy a
property in January, it is profitable to conclude the purchase contract
already in December. This results in a shift of property sales from January to
December, and thus more sales in December and fewer sales in January than
expected. However, does this also impact the risk of the real-estate market?

Figure 2: Overview of the data set. Left: Average house prices for each of the
190 zip-code areas over the period from January 1995 to December 2015. Right:
Total number of real-estate transactions in all zip-code areas of Berlin. The
vertical bars indicate the time points of changes in the real-estate transfer
taxes in Berlin or the surrounding state of Brandenburg (red: increase from
3.5 to 4.5 % in Berlin; green: increase from 3.5 to 5 % in Brandenburg; blue:
increase from 4.5 to 5 % in Berlin; yellow: increase from 5 to 6 % in Berlin)

For the empirical analysis, we use monthly log-returns of the average sales
prices of all condominium sales in all postcode regions of Berlin from January
1995 to December 2015 (see Figure 2, left). The relative price per square
meter is determined for each zip-code area from an average of 6.31 sales per
month. The specific location of the German capital Berlin in the centre of
another federal state, Brandenburg, makes it a very intriguing example. The
surrounding area of Berlin is very well connected to the city center by public
transport and infrastructure, such that exogenous effects such as tax changes
in Brandenburg may have an impact on Berlin and vice versa. Every real estate
purchase in Germany is subject to the real estate transfer tax, which must be
paid once at the time of purchase. The amount of tax depends on the purchase
price. Until 1.9.2006, a unified tax rate of 3.5 per cent was applied in
Germany as a whole. Afterwards, each federal state could set its tax rate, and
there were gradual increases in all federal states. Specifically, Berlin
increased the tax rates from 3.5 to 4.5 per cent on 1.1.2007, from 4.5 to 5
per cent on 1.4.2012, and from 5 to 6 per cent on 1.1.2014. The shifting
effects described above can also be observed for Berlin, as it is shown in
Figure 2. In addition, an end-of-year effect is clearly visible due to other
accounting and tax reasons. This motivates why we would expect different
market risks at the end and beginning of a year. To estimate these temporal
effects, we consider a model without temporal fixed effects (i.e.,
$\alpha_{t0}=0$ for all $t$ in (2)) and model the temporal effects by
including yearly and monthly indicator variables as regressors.

Table 4: Estimated parameters of the dynamic ARCH process and diagnostic measures of the residuals. |  | Parameter | Estimate | Standard error | T-statistics
---|---|---|---|---|---
_Regressive effects_
|  Total number of transactions | $\beta_{1}$ | -0.1356 | 0.0184 | -7.3509
_End-of-year effects_
|  December | $\beta_{2}$ | -0.0516 | 0.0465 | -1.1098
|  January | $\beta_{3}$ | 0.2152 | 0.0729 | 2.9527
|  February | $\beta_{4}$ | 0.1465 | 0.0555 | 2.6379
_Yearly effects_
|  1997 | $\beta_{5}$ | -0.1447 | 0.0603 | -2.4002
|  1998 | $\beta_{6}$ | -0.2112 | 0.0735 | -2.8739
|  1999 | $\beta_{7}$ | -0.2003 | 0.0708 | -2.8292
|  2000 | $\beta_{8}$ | -0.1136 | 0.0585 | -1.9413
|  2001 | $\beta_{9}$ | -0.1194 | 0.0600 | -1.9897
|  2011 | $\beta_{10}$ | -0.1244 | 0.0674 | -1.8466
|  2012 | $\beta_{11}$ | -0.1161 | 0.0684 | -1.6974
|  2013 | $\beta_{12}$ | -0.1204 | 0.0684 | -1.7606
_Tax effects_
|  Month before tax increase | $\beta_{13}$ | -0.0559 | 0.1120 | -0.4988
_Spatiotemporal effects_
|  Spatial interaction (contiguity-based) | $\rho$ | 0.4032 | 0.1450 | 2.7801
|  Temporal interaction (first time lag) | $\gamma$ | 0.1913 | 0.0051 | 37.4971
|  Spatiotemporal interaction | $\delta$ | -0.0737 | 0.0313 | -2.3547
_Model diagnostics_
|  BIC |  | 82572.28 |  |
|  n |  | 190 |  |
|  T |  | 239 |  |
|  Percentage of locations with significantly (temp.) autocorrelated errors ($\alpha=5$ %) |  | 18 |  |
|  Percentage of time points with significantly (spatially) positive autocorrelated errors (Moran’s $I$, $\alpha=5$ %) |  | 0 |  |

We estimate a first-order version of our model in (2). We specify the spatial
weights matrix (row-standardized) based on the queen contiguity scheme, where
all adjacent neighbors are equally weighted. In Table 4 and Figure 3, we
present the estimated parameters of our dynamic spatiotemporal ARCH model and
the estimated volatility, respectively. The included covariates were selected
by stepwise excluding regressors, such that the Bayesian information criterion
is minimized. The overall market dynamic measured by the total numbers of
real-estate transactions has the most significant effect on the volatility
(see Table 4). The more transactions, the lower the log-volatility. In
addition, we observe lower volatilities at the end of each year, and
significantly higher volatilities in January and February. These effects
decrease from January to February ($0.2152$ to $0.1465$). The March effects
were already insignificant. Further, we observe two periods of significantly
lower risks compared to the remaining periods, namely 1997-2001 and 2011-2013.
The anticipated tax effects are not significant though. Thus, we could not
find evidence that the legal changes in the taxation framework affect the
volatility of log-returns.

As expected, we also see significant spatial and spatiotemporal spill-over
effects. The spatial ARCH parameter $\hat{\rho}=0.4032$ is on moderate level.
That is, an increase in the log-squared return in one location instantaneously
increases the log-volatility in the adjacent regions. Further, the temporal
dependence is composed of a purely temporal lag ($\hat{\gamma}=0.1913$) and a
spatiotemporal lag ($\hat{\delta}=-0.0737$), which in total indicate a
moderate temporal persistence.

When analyzing the estimated volatility, we clearly see temporal patterns with
reduced risks during the two above-mentioned periods. More interestingly,
there are several regions of higher volatility, mostly located at the outer
zip-code regions in the North and North-West, as shown in the top panels of
Figure 3, where the averaged volatility estimates are depicted. Moreover, the
average volatility over all zip-codes changes over time as shown in the first
figure of the top panels of Figure 3. Finally, to illustrate the estimated
$h_{it}$’s for one selected location, we show the results for Berlin-Tempelhof
(zip-code of the closed airport Berlin-Tempelhof) in April 2012, the month
after the increase of the real-estate transfer taxes from 4.5 to 5 %. We do
not see different patterns in the volatility estimates after the airport was
closed, marked by the dashed black line.

Figure 3: Estimated conditional volatility $\hat{h}_{it}$. Top left: The
spatially averaged volatilities levels are displayed over time
$\frac{1}{n}\sum_{i=1}^{n}\hat{h}_{it}$. The vertical grey line marks the Euro
introduction in January 2002. Top right: The temporally averaged volatilities
levels are displayed on the map $\frac{1}{T}\sum_{t=1}^{T}\hat{h}_{it}$.
Bottom left: $\hat{h}_{it}$ of one selected region over time (Berlin
Tempelhof, zip-code of the former airport, marked by red cross in the bottom
right plot). Bottom right: $\hat{h}_{it}$ at one selected time point (April
2012, first month after increase of real-estate transfer taxes from 4.5 to 5
%, shown by dashed red line in the bottom-left plot).

##  6 Conclusion

In this paper, we introduced a dynamic spatiotemporal ARCH model that allows
for unobserved heterogeneity over time and space. The model can be used to
describe the spatiotemporal clustering effect in the volatility of a random
process. As typically observed for spatial data, the model allows for
instantaneous spill-over effects across space, which is the main difference to
multivariate time-series GARCH models. In the latter case, spatial
interactions would only occur after one time lag. In addition to these
instantaneous spatial effects, the model includes temporal and spatiotemporal
autoregressive effects of the log-squared returns in the log-volatility
equation. While the temporal effect measures the dependence between the
current and past observation of the same spatial unit, the spatiotemporal
coefficients describe the dependence between the observation in one location
and its past observations at neighboring locations.

For our suggested dynamic spatiotemporal ARCH model, we obtain an estimation
equation by applying a log-square transformation together with an orthonormal
and a deviation from group-mean operator to eliminate the fixed effects. We
introduced a GMM estimation approach based on a set of linear and quadratic
moment functions of the transformed process. We establish the consistency and
asymptotic normality of our suggested GMM estimator under fairly general
assumptions for large and finite $T$ cases. Moreover, when the number of time
periods is large, we present an optimal set of moment functions that leads to
an efficient estimator.

We investigated the finite-sample performance of our suggested estimator in a
series of Monte-Carlo simulations under different model settings and error
distributions. Overall, the simulation results are in line with our
theoretical claims. In an empirical application, we illustrated the use of our
model for the log-returns of the intra-city real-estate prices in Berlin over
the period 1995 - 2015. Our estimation results show that the spatial, temporal
and spatiotemporal lags of the log-squared returns have statistically
significant effect on the log-volatility. This leads to temporal and spatial
spill-over effects. We showed that the average volatility of log-returns over
space and time varies significantly. Finally, our model allows us to estimate
the market risk in terms of the volatility in each location and time point.

In future studies, our model can be extended in a number of ways. First, we
considered additive time and space fixed effects in the log-volatility
equation. Instead of this additive structure, a log-volatility equation that
includes interactive fixed effects can be studied. Second, the spatial and
spatiotemporal lags in the log-volatility equation can be formulated with
time-varying spatial weights matrices. Finally, we can also allow for
potential endogeneity in the spatial weights instead of exogenous spatial
weights. All of these extensions can be explored in future studies.

Appendix

##  Appendix A Some Useful Lemmas

In this section, we provide four lemmas that are essential for our main
results.

######  Lemma 1.

Let $\mathcal{F}_{t-1}$ be the $\sigma$-algebra generated by
$\left(\mathbf{Y}_{0},\ldots,\mathbf{Y}_{t-1}\right)$ conditional on
$\left(\mathbf{X}_{1},\ldots,\mathbf{X}_{T},\boldsymbol{\mu}_{0},\boldsymbol{\alpha}_{0}\right)$,
and $\mathbf{B}$ be a non-stochastic $n\times n$ matrix that has row and
column sums uniformly bounded in absolute value. Then, under Assumption 1, we
have the following results.

  1. 1.

$\mathbb{E}\left(\mathbf{U}^{*^{\prime}}_{t}\mathbf{B}\mathbf{U}^{*}_{s}|\mathcal{F}_{t-1}\right)=0$
for $t\neq s$.

  2. 2.

$\frac{1}{N}\mathbf{U}^{{}^{\prime}}_{N}\mathbf{B}_{N}\mathbf{U}_{N}=\frac{\sigma^{2}_{0}}{n}\mathrm{tr}\left(\mathbf{B}\right)+O_{p}\left(N^{-1/2}\right)=O_{p}(1)$,
where $\mathbf{B}_{N}=\left(\mathbf{I}_{T-1}\otimes\mathbf{B}\right)$.

  3. 3.

$\mathbb{E}\left(\left(\mathbf{U}^{*^{\prime}}_{t}\mathbf{B}\mathbf{U}^{*}_{t}\right)^{2}\right)=\left(\mu_{4}-3\sigma^{4}_{0}\right)c^{4}_{t}\left(1+\frac{1}{(T-t)^{3}}\right)\text{vec}^{{}^{\prime}}_{D}(\mathbf{B})\text{vec}_{D}(\mathbf{B})+\sigma^{4}_{0}\left(\mathrm{tr}^{2}(\mathbf{B})+\mathrm{tr}\left(\mathbf{B}\mathbf{B}^{s}\right)\right)$,
where $c_{t}=\left(\frac{T-t}{T-t+1}\right)^{1/2}$.

  4. 4.

Under Assumption 6, we have
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\sum_{t=1}^{T-1}\mathbf{Q}^{{}^{\prime}}_{t}\mathbf{B}\mathbf{U}^{*}_{t}=\mathbf{0}$,
where $\mathbf{Q}_{t}$ is the IV matrix.

  5. 5.

Under Assumption 5, we have
$\frac{1}{N}\mathbf{Y}^{{}^{\prime}}_{N,-1}\mathbf{B}_{N}\mathbf{U}_{N}-\mathbb{E}\left(\frac{1}{N}\mathbf{Y}^{{}^{\prime}}_{N,-1}\mathbf{B}_{N}\mathbf{U}_{N}\right)=O_{p}\left(N^{-1/2}\right)$,
where
$\mathbb{E}\left(\frac{1}{N}\mathbf{Y}^{{}^{\prime}}_{N,-1}\mathbf{B}_{N}\mathbf{U}_{N}\right)=\frac{\sigma^{2}_{0}}{N}\mathrm{tr}\left(\mathbf{B}^{{}^{\prime}}\left(\sum_{h=1}^{T-1}\left(1-\frac{h}{T}\right)\mathbf{A}^{h-1}\right)\mathbf{S}^{-1}\right)$
is $O(T^{-1})$.

###### Proof.

See Lemma 2 in Lee and Yu, (2014). ∎

######  Lemma 2.

Under Assumption 1, we have the following results.

  1. 1.

For $l,j=1,2,\ldots,m$, we have

$\displaystyle\mathbb{E}\left(\mathbf{U}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}\cdot\mathbf{U}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}_{jN}\mathbf{J}_{N}\mathbf{U}_{N}\right)$
$\displaystyle=\sigma^{4}_{0}\mathrm{tr}\left(\mathbf{J}_{N}\mathbf{P}_{lN}\mathbf{J}_{N}\left(\mathbf{J}_{N}\mathbf{P}_{jN}\mathbf{J}_{N}\right)^{s}\right)+(\mu_{4}-3\sigma^{4}_{0})\text{vec}^{{}^{\prime}}_{D}\left(\mathbf{J}_{N}\mathbf{P}_{lN}\mathbf{J}_{N}\right)\text{vec}_{D}\left(\mathbf{J}_{N}\mathbf{P}_{jN}\mathbf{J}_{N}\right).$

  2. 2.

For $l=1,2,\ldots,m$, we have

$\displaystyle\mathbb{E}\left(\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}\cdot\mathbf{U}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}\right)=\mathbf{0}_{q\times
1}.$

###### Proof.

This lemma is a simple extension of Lemma 2 of Lee and Yu, (2014) to our
setting. Therefore, we omit its proof. ∎

######  Lemma 3.

Let $\mathbf{p}_{t}$ be $n\times 1$ column vector from the IV matrix
$\mathbf{Q}_{t}$ given in Assumption 6. Consider
$a_{t}=\mathbf{q}^{{}^{\prime}}_{t}\mathbf{U}^{*}_{t}+\mathbf{U}^{*^{\prime}}_{t}\mathbf{B}\mathbf{U}^{*}_{t}-\sigma^{2}_{0}\mathrm{tr}(\mathbf{B})$,
where $\mathbf{B}$ is a non-stochastic $n\times n$ matrix that has row and
column sums uniformly bounded in absolute value. Using Lemma 2, we have
$\mathrm{Var}\left(\sum_{t=1}^{T-1}a_{t}\right)=\sigma^{2}_{0}\sum_{t=1}^{T-1}\mathbf{q}^{{}^{\prime}}_{t}\mathbf{q}_{t}+T\left(\mu_{4}-3\sigma^{4}_{0}\right)\sum_{i=1}^{n}b^{2}_{ii}+T\sigma^{4}_{0}\mathrm{tr}\left(\mathbf{B}\mathbf{B}^{s}\right)$.
If
$\left\\{\frac{1}{N}\mathrm{Var}\left(\sum_{t=1}^{T-1}a_{t}\right)\right\\}$
is bounded away from zero, then

$\frac{\sum_{t=1}^{T-1}a_{t}}{\mathrm{Var}^{1/2}\left(\sum_{t=1}^{T-1}a_{t}\right)}\xrightarrow{d}N(0,\,1).$

###### Proof.

See the CLT results in Lee and Yu, (2014); Yu et al., (2008). ∎

######  Lemma 4.

Let
$\mathbf{Z}^{*}_{s}=\left(\mathbf{Y}^{*}_{s-1},\mathbb{M}\mathbf{Y}^{*}_{s-1},\mathbf{X}_{s}\right)$
and $c_{t}=\left(\frac{T-t}{T-t+1}\right)^{1/2}$. Then,

$\displaystyle\mathbf{Y}^{**,-1}_{t-1}$
$\displaystyle=c_{t}\left(\left(\mathbf{I}_{n}-\frac{1}{T-t}\sum_{h=1}^{T-t}\mathbf{A}^{h}\right)\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\left(\mathbf{X}_{r}\boldsymbol{\beta}_{0}+\alpha_{r,0}\mathbf{1}_{n}\right)\right)$
$\displaystyle-
c_{t}\frac{1}{(T-t)(t-1)}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\sum_{s=1}^{t-1}\left(\mathbf{S}\mathbf{Y}^{*}_{s}-\mathbf{Z}^{*}_{s}\boldsymbol{\eta}_{0}-\alpha_{s0}\mathbf{1}_{n}\right)$
(31) $\displaystyle-
c_{t}\frac{1}{(T-t)(t-1)}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\sum_{s=1}^{t-1}\mathbf{U}_{s}-c_{t}\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\mathbf{U}_{r}.$

###### Proof.

Using (7), we have

$\displaystyle\mathbf{Y}^{*}_{t}$
$\displaystyle=\mathbf{A}\mathbf{Y}^{*}_{t-1}+\mathbf{S}^{-1}\left(\mathbf{X}_{t}\boldsymbol{\beta}_{0}+\tilde{\boldsymbol{\mu}}_{0}+\alpha_{t0}\mathbf{1}_{n}+\mathbf{U}_{t}\right),$
(32)

where
$\tilde{\boldsymbol{\mu}}_{0}=\boldsymbol{\mu}+\mu_{\operatorname{\varepsilon}}\mathbf{1}_{n}$.
Using (32), we can expand $\mathbf{Y}^{*}_{t+h}$ for $h\geq 0$ as

$\displaystyle\mathbf{Y}^{*}_{t+h}$
$\displaystyle=\mathbf{A}^{h+1}\mathbf{Y}^{*}_{t-1}+\sum_{j=0}^{h}\mathbf{A}^{j}\mathbf{S}^{-1}\left(\mathbf{X}_{t+h-j}\boldsymbol{\beta}_{0}+\tilde{\boldsymbol{\mu}}_{0}+\alpha_{t+h-j,0}\mathbf{1}_{n}+\mathbf{U}_{t+h-j}\right),$
(33)

Therefore, we have

$\displaystyle\mathbf{Y}^{*}_{s}$
$\displaystyle=\mathbf{A}^{s-t+1}\mathbf{Y}^{*}_{t-1}+\sum_{j=0}^{s-t}\mathbf{A}^{j}\mathbf{S}^{-1}\left(\mathbf{X}_{s-j}\boldsymbol{\beta}_{0}+\tilde{\boldsymbol{\mu}}_{0}+\alpha_{s-j,0}\mathbf{1}_{n}+\mathbf{U}_{s-j}\right).$
(34)

Then, we can express $\sum_{s=t}^{T-1}\mathbf{Y}^{*}_{s}$ as

$\displaystyle\sum_{s=t}^{T-1}\mathbf{Y}^{*}_{s}$
$\displaystyle=\sum_{s=t}^{T-1}\mathbf{A}^{s-t+1}\mathbf{Y}^{*}_{t-1}+\sum_{s=t}^{T-1}\sum_{j=0}^{s-t}\mathbf{A}^{j}\mathbf{S}^{-1}\left(\mathbf{X}_{s-j}\boldsymbol{\beta}_{0}+\tilde{\boldsymbol{\mu}}_{0}+\alpha_{s-j,0}\mathbf{1}_{n}+\mathbf{U}_{s-j}\right)$
$\displaystyle=\sum_{h=1}^{T-t}\mathbf{A}^{h}\mathbf{Y}^{*}_{t-1}+\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\left(\mathbf{X}_{r}\boldsymbol{\beta}_{0}+\tilde{\boldsymbol{\mu}}_{0}+\alpha_{r,0}\mathbf{1}_{n}+\mathbf{U}_{r}\right).$
(35)

Thus, using
$\mathbf{Y}^{**,-1}_{t-1}=c_{t}\left(\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{s=t}^{T-1}\mathbf{Y}^{*}_{s}\right)$,
we obtain

$\displaystyle\mathbf{Y}^{**,-1}_{t-1}$
$\displaystyle=c_{t}\left(\left(\mathbf{I}_{n}-\frac{1}{T-t}\sum_{h=1}^{T-t}\mathbf{A}^{h}\right)\mathbf{Y}^{*}_{t-1}-\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\left(\mathbf{X}_{r}\boldsymbol{\beta}_{0}+\alpha_{r,0}\mathbf{1}_{n}\right)\right)$
$\displaystyle-
c_{t}\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\tilde{\boldsymbol{\mu}}_{0}-c_{t}\frac{1}{T-t}\sum_{r=t}^{T-1}\left(\sum_{h=0}^{T-r-1}\mathbf{A}^{h}\right)\mathbf{S}^{-1}\mathbf{U}_{r}.$
(36)

Note that we can express $\tilde{\boldsymbol{\mu}}_{0}$ as

$\displaystyle\tilde{\boldsymbol{\mu}}_{0}=\frac{1}{t-1}\sum_{s=1}^{t-1}\left(\mathbf{S}\mathbf{Y}^{*}_{s}-\mathbf{Z}^{*}_{s}\boldsymbol{\eta}_{0}-\alpha_{s0}\mathbf{1}_{n}-\mathbf{U}_{s}\right),$
(37)

where
$\mathbf{Z}^{*}_{s}=\left(\mathbf{Y}^{*}_{s-1},\mathbb{M}\mathbf{Y}^{*}_{s-1},\mathbf{X}_{s}\right)$.
Then, substituting (37) into (A) yields the result. ∎

##  Appendix B Details on Identification

In the GMM setting, the parameter vector is identified if
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{g}_{N}(\boldsymbol{\theta})=\mathbf{0}$
has the unique solution $\boldsymbol{\theta}_{0}$. Consider the linear moment
function
$\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})$.
From (9), we have
$\mathbf{J}_{n}\mathbf{U}^{*}_{t}(\boldsymbol{\theta})=\mathbf{J}_{n}\left(\mathbf{S}(\boldsymbol{\rho})\mathbf{Y}^{**}_{t}-\mathbf{Z}^{**}_{t}\boldsymbol{\eta}\right)$,
and
$\mathbf{Y}^{**}_{t}=\mathbf{S}^{-1}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}+\mathbf{U}^{*}_{t}\right)$.
Thus, we have

$\displaystyle\mathbf{J}_{n}\mathbf{U}^{*}_{t}(\boldsymbol{\theta})=\mathbf{J}_{n}\left(\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}+\mathbf{U}^{*}_{t}\right)-\mathbf{Z}^{**}_{t}\boldsymbol{\eta}\right)$

Note that we can express $\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}$ in the
following way:

$\displaystyle\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}$
$\displaystyle=\left(\mathbf{I}_{n}-\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\right)\mathbf{S}^{-1}=\left(\mathbf{I}_{n}+\mathbf{S}-\mathbf{S}-\sum_{l=1}^{p}\rho_{l}\mathbf{M}_{l}\right)\mathbf{S}^{-1}$
$\displaystyle=\left(\mathbf{S}-\sum_{l=1}^{p}(\rho_{l}-\rho_{l0})\mathbf{M}_{l}\right)\mathbf{S}^{-1}=\left(\mathbf{I}_{n}-\sum_{l=1}^{p}(\rho_{l}-\rho_{l0})\mathbf{G}_{l}\right).$

Then, we can express $\mathbf{J}_{n}\mathbf{U}^{*}_{t}(\boldsymbol{\theta})$
as

$\displaystyle\mathbf{J}_{n}\mathbf{U}^{*}_{t}(\boldsymbol{\theta})=\mathbf{J}_{n}\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}\right)+\mathbf{J}_{n}\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}\mathbf{U}^{*}_{t}-\mathbf{J}_{n}\mathbf{Z}^{**}_{t}\boldsymbol{\eta}$
$\displaystyle=\mathbf{J}_{n}\left(\mathbf{I}_{n}-\sum_{l=1}^{p}(\rho_{l}-\rho_{l0})\mathbf{G}_{l}\right)\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}\right)+\mathbf{J}_{n}\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}\mathbf{U}^{*}_{t}-\mathbf{J}_{n}\mathbf{Z}^{**}_{t}\boldsymbol{\eta}$
$\displaystyle=\mathbf{J}_{n}\left(\mathbf{Z}^{**}_{t}\left(\boldsymbol{\eta}_{0}-\boldsymbol{\eta}\right)-\sum_{l=1}^{p}(\rho_{l}-\rho_{l0})\mathbf{G}_{l}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}\right)\right)+\mathbf{J}_{n}\mathbf{S}(\boldsymbol{\rho})\mathbf{S}^{-1}\mathbf{U}^{*}_{t}.$

Define
$\mathbf{L}_{r,t}=\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t}\mathbf{1}_{n}\right)$,
$\mathbf{L}_{t}=\left(\mathbf{L}_{1,t},\ldots,\mathbf{L}_{p,t}\right)$, and
$\mathbf{L}_{N}=\left(\mathbf{L}^{{}^{\prime}}_{1},\ldots,\mathbf{L}^{{}^{\prime}}_{T-1}\right)^{{}^{\prime}}$.
Thus, we can write
$\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})$
as

$\displaystyle\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})$
$\displaystyle=\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{Z}_{N}\left(\boldsymbol{\eta}_{0}-\boldsymbol{\eta}\right)+\mathbf{L}_{N}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\rho}\right)+\mathbf{S}_{N}(\boldsymbol{\rho})\mathbf{S}^{-1}_{N}\mathbf{U}^{*}_{t}\right)$
$\displaystyle=\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{Z}_{N},\,\mathbf{L}_{N}\right)\left(\left(\boldsymbol{\eta}_{0}-\boldsymbol{\eta}\right)^{{}^{\prime}},\left(\boldsymbol{\rho}_{0}-\boldsymbol{\rho}\right)^{{}^{\prime}}\right)^{{}^{\prime}}+\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{S}_{N}(\boldsymbol{\rho})\mathbf{S}^{-1}_{N}\mathbf{U}_{N}.$

By Lemma 1, we have
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{S}_{N}(\boldsymbol{\rho})\mathbf{S}^{-1}_{N}\mathbf{U}_{N}=\mathbf{0}$.
Thus, we require that the equation
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{Z}_{N},\,\mathbf{L}_{N}\right)\left(\left(\boldsymbol{\eta}_{0}-\boldsymbol{\eta}\right)^{{}^{\prime}},\left(\boldsymbol{\rho}_{0}-\boldsymbol{\rho}\right)^{{}^{\prime}}\right)^{{}^{\prime}}=\mathbf{0}$
should have a unique solution at $\boldsymbol{\theta}_{0}$. If
$\operatorname{plim}_{n\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{Z}_{N},\,\mathbf{L}_{N}\right)$
has the full column rank, then we will have a unique solution at
$\boldsymbol{\theta}_{0}$.

##  Appendix C Proof of Theorem 1

In this section, we will first show that
$\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\mathbf{D}_{1N}+\mathbf{D}_{2N}+O_{p}(N^{-1/2})$.
Recall that

$\displaystyle\mathbf{g}_{N}(\boldsymbol{\theta})=\begin{pmatrix}\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{1N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\vdots\\\
\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}_{mN}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\end{pmatrix},$
(38)

where
$\boldsymbol{\theta}=\left(\boldsymbol{\rho}^{{}^{\prime}},\boldsymbol{\eta}^{{}^{\prime}}\right)^{{}^{\prime}}$,
$\mathbf{U}_{N}(\boldsymbol{\theta})=\left(\mathbf{U}^{*^{\prime}}_{1}(\boldsymbol{\theta}),\ldots,\mathbf{U}^{*^{\prime}}_{T-1}(\boldsymbol{\theta})\right)^{{}^{\prime}}$
and
$\mathbf{U}^{*}_{t}(\boldsymbol{\theta})=\left(\mathbf{S}(\boldsymbol{\rho})\mathbf{Y}^{**}_{t}-\mathbf{Z}^{**}_{t}\boldsymbol{\eta}-\alpha^{*}_{t}\mathbf{1}_{n}\right)$.
Note that
$\frac{\partial\mathbf{U}^{*}_{t}(\boldsymbol{\theta})}{\partial\rho_{j}}=-\mathbf{M}_{j}\mathbf{Y}^{**}_{t}$
for $j=1,\ldots,p$ and
$\frac{\partial\mathbf{U}^{*}_{t}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}=-\mathbf{Z}^{**}_{t}$.
Then, the components of
$\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}$
are

$\displaystyle\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}=-\begin{pmatrix}\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}^{s}_{1N}\mathbf{J}_{N}\mathbf{Z}_{N}\\\
\vdots\\\
\mathbf{U}^{{}^{\prime}}_{N}(\boldsymbol{\theta})\mathbf{J}_{N}\mathbf{P}^{s}_{mN}\mathbf{J}_{N}\mathbf{Z}_{N}\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Z}_{N}\end{pmatrix},\quad\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\rho}^{{}^{\prime}}}=\left(\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\rho_{1}},\ldots,\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\rho_{p}}\right),$
(39)

where

$\displaystyle\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\rho_{j}}=-\begin{pmatrix}\left(\mathbf{M}_{jN}\mathbf{Y}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{1N}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\vdots\\\
\left(\mathbf{M}_{jN}\mathbf{Y}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{mN}\mathbf{J}_{N}\mathbf{U}_{N}(\boldsymbol{\theta})\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{M}_{jN}\mathbf{Y}_{N}\end{pmatrix},$
(40)

with $\mathbf{M}_{jN}=(\mathbf{I}_{T-1}\otimes\mathbf{M}_{j})$ for
$j=1,2,\ldots,p$. Next, we will determine the probability limit of
$\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}$.
Consider
$\mathbf{U}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{Z}_{N}$
in
$\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}$.
Recall that
$\mathbf{Z}^{**}_{t}=\left(\mathbf{Y}^{**,-1}_{t-1},\mathbb{M}\mathbf{Y}^{**,-1}_{t-1},\mathbf{X}^{*}_{t}\right)$.
Then, we can express
$\mathbf{Z}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
in the following way

$\displaystyle\mathbf{Z}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
$\displaystyle=\sum_{t=1}^{T-1}\mathbf{Z}^{**^{\prime}}_{t}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}=\begin{pmatrix}\sum_{t=1}^{T-1}\mathbf{Y}^{**,-1^{\prime}}_{t-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\\\
\sum_{t=1}^{T-1}\left(\mathbb{M}\mathbf{Y}^{**,-1}_{t-1}\right)^{{}^{\prime}}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\\\
\sum_{t=1}^{T-1}\mathbf{X}^{*^{\prime}}_{t}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\end{pmatrix}.$
(41)

Using Lemma 1, it can be shown that

$\displaystyle\mathbb{E}\left(\frac{1}{N}\sum_{t=1}^{T-1}\mathbf{\mathbf{Y}}^{**,-1}_{t-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\right)=\frac{\sigma^{2}_{0}}{NT}\mathrm{tr}\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\right)$
$\displaystyle\mathbb{E}\left(\frac{1}{N}\sum_{t=1}^{T-1}\left(\mathbf{M}_{j}\mathbf{Y}^{**,-1}_{t-1}\right)^{{}^{\prime}}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\right)=\frac{\sigma^{2}_{0}}{NT}\mathrm{tr}\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{M}_{j}\right),$
$\displaystyle\mathbb{E}\left(\frac{1}{N}\sum_{t=1}^{T-1}\mathbf{X}^{*^{\prime}}_{t}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}^{*}_{t}\right)=\mathbf{0}_{k\times
1}.$

for $l=1,2,\ldots,m$ and $j=1,2,\ldots,p$. Define
$\mathbf{b}_{l\gamma}=\frac{\sigma^{2}_{0}}{N}\mathrm{tr}\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\right)$,
$\mathbf{b}_{l\delta_{j}}=\frac{\sigma^{2}_{0}}{N}\mathrm{tr}\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{M}_{s}\right)$
for $j=1,2,\ldots,p$, and

$\displaystyle\mathbf{b}_{l\boldsymbol{\eta}}=\left(\mathbf{b}_{l\gamma},\mathbf{b}_{l\delta_{1}},\mathbf{b}_{l\delta_{2}},\ldots,\mathbf{b}_{l\delta_{p}},\mathbf{0}_{k\times
1}\right).$ (42)

Then, Lemma 1 ensures that

$\displaystyle\frac{1}{N}\mathbf{Z}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}=\frac{1}{T}\mathbf{b}^{{}^{\prime}}_{l\boldsymbol{\eta}}+O_{P}(N^{-1/2}),$
(43)

for $l=1,2,\ldots,m$. Let
$\mathbf{b}_{\boldsymbol{\delta}}=\left(\mathbf{b}^{{}^{\prime}}_{1\boldsymbol{\eta}},\mathbf{b}^{{}^{\prime}}_{l\boldsymbol{\eta}},\ldots,\mathbf{b}^{{}^{\prime}}_{m\boldsymbol{\eta}}\right)^{{}^{\prime}}$.
Then, we have

$\displaystyle\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}=-\frac{1}{N}\begin{pmatrix}\mathbf{0}_{m\times
k_{z}}\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Z}_{N}\end{pmatrix}+\frac{1}{T}\begin{pmatrix}\mathbf{b}_{\boldsymbol{\delta}}\\\
\mathbf{0}_{k_{q}\times k_{z}}\end{pmatrix}+O_{p}(N^{-1/2}).$ (44)

Next, we consider
$\left(\mathbf{M}_{rN}\mathbf{Y}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
in
$\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\rho}^{{}^{\prime}}}$.
Recall that
$\mathbf{M}_{r}\mathbf{Y}^{**}_{t}=\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right)+\mathbf{G}_{r}\mathbf{U}^{*}_{t}$
for $r=1,\ldots,p$. This expression implies that
$\mathbf{M}_{rN}\mathbf{Y}_{N}=\mathbf{L}_{rN}+\mathbf{G}_{rN}\mathbf{U}_{N}$,
where $\mathbf{G}_{rN}=(\mathbf{I}_{T-1}\otimes\mathbf{G}_{r})$ and
$\mathbf{L}_{rN}=\left(\mathbf{L}^{{}^{\prime}}_{r1},\mathbf{L}^{{}^{\prime}}_{r2},\ldots,\mathbf{L}^{{}^{\prime}}_{r,T-1}\right)^{{}^{\prime}}$
with
$\mathbf{L}_{rt}=\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right)$.
Thus, we can express
$\left(\mathbf{M}_{rN}\mathbf{Y}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
as

$\displaystyle\left(\mathbf{M}_{rN}\mathbf{Y}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}=\mathbf{L}^{{}^{\prime}}_{rN}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}+\left(\mathbf{G}_{rN}\mathbf{U}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$

By Lemma 1, we have
$\frac{1}{N}\left(\mathbf{G}_{rN}\mathbf{U}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}=\frac{\sigma^{2}_{0}}{N}\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{rN}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\right)+O_{p}(N^{-1/2})$
for $r=1,2,\ldots,p$ and $l=1,2,\ldots,m$. Note that we can express
$\mathbf{L}^{{}^{\prime}}_{rN}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
as

$\displaystyle\mathbf{L}^{{}^{\prime}}_{rN}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}=\sum_{t=1}^{T-1}\mathbf{L}^{{}^{\prime}}_{rt}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}_{t}=\sum_{t=1}^{T-1}\left(\mathbf{G}_{r}\left(\mathbf{Z}^{**}_{t}\boldsymbol{\eta}_{0}+\alpha^{*}_{t0}\mathbf{1}_{n}\right)\right)^{{}^{\prime}}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}_{t}$
$\displaystyle=\sum_{t=1}^{T-1}\left(\mathbf{G}_{r}\left(\gamma_{0}+\sum_{s=1}^{p}\delta_{s0}\mathbf{M}_{s}\right)\mathbf{Y}^{**,-1}_{t-1}+\mathbf{G}_{r}\mathbf{X}^{*}_{t}\boldsymbol{\beta}_{0}+\alpha^{*}_{t0}\mathbf{G}_{r}\mathbf{1}_{n}\right)^{{}^{\prime}}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\mathbf{U}_{t}$

Then, using Lemma 1, we obtain

$\displaystyle\frac{1}{N}\mathbf{L}^{{}^{\prime}}_{rN}\mathbf{J}_{N}\mathbf{P}^{s}_{lN}\mathbf{J}_{N}\mathbf{U}_{N}$
$\displaystyle=\frac{\sigma^{2}_{0}}{NT}\mathrm{tr}\left(\left(\mathbf{G}_{r}\left(\gamma_{0}+\sum_{j=1}^{p}\delta_{j0}\mathbf{M}_{j}\right)\right)\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\right)\right).$

Let
$\mathbf{b}_{l\rho_{r}}=\frac{\sigma^{2}_{0}}{NT}\mathrm{tr}\left(\left(\mathbf{G}_{r}\left(\gamma_{0}+\sum_{j=1}^{q}\delta_{j0}\mathbf{M}_{j}\right)\right)\left(\sum_{h=1}^{T-1}(T-h)\mathbf{A}^{h-1}\mathbf{S}^{-1}\mathbf{J}_{n}\mathbf{P}^{s}_{l}\mathbf{J}_{n}\right)\right)$
for $r=1,2,\ldots,p$ and $l=1,2,\ldots,m$. Define the following matrices:

$\displaystyle\mathbf{C}_{N}=\begin{pmatrix}\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{s}_{1N}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{s}_{1N}\mathbf{J}_{N}\right)\\\
\vdots&\ddots&\vdots\\\
\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{s}_{mN}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{s}_{mN}\mathbf{J}_{N}\right)\end{pmatrix},$
(45)
$\displaystyle\mathbf{b}_{\boldsymbol{\rho}}=\begin{pmatrix}\mathbf{b}_{1\rho_{1}}&\mathbf{b}_{1\rho_{2}}&\ldots&\mathbf{b}_{1\rho_{p}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{b}_{m\rho_{1}}&\mathbf{b}_{m\rho_{2}}&\ldots&\mathbf{b}_{m\rho_{p}}\\\
\mathbf{0}_{k_{q}\times 1}&\mathbf{0}_{k_{q}\times
1}&\ldots&\mathbf{0}_{k_{q}\times 1}\end{pmatrix}.$ (46)

Let $\mathbf{L}_{N}=\left(\mathbf{L}_{1N},\ldots,\mathbf{L}_{pN}\right)$. The
preceding analysis indicates that

$\displaystyle\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}=-\frac{1}{N}\begin{pmatrix}\sigma^{2}_{0}\mathbf{C}_{N}\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{L}_{N}\end{pmatrix}+\frac{1}{T}\begin{pmatrix}\mathbf{b}_{\boldsymbol{\rho}}\\\
\mathbf{0}_{k_{q}\times p}\end{pmatrix}.$ (47)

Substituting (44) and (47) into
$\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\left(\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\rho}^{{}^{\prime}}},\,\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\eta}^{{}^{\prime}}}\right)$,
we obtain

$\displaystyle\frac{1}{N}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}=\mathbf{D}_{1N}+\mathbf{D}_{2N}+o_{p}(1),$
(48)

where

$\displaystyle\mathbf{D}_{1N}=-\frac{1}{N}\begin{pmatrix}\sigma^{2}_{0}\mathbf{C}_{N}&\mathbf{0}_{m\times
k_{z}}\\\
\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{L}_{N}&\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Z}_{N}\end{pmatrix}=O(1),$
(49)
$\displaystyle\mathbf{D}_{2N}=\frac{1}{T}\begin{pmatrix}\mathbf{b}_{\boldsymbol{\rho}}&\mathbf{b}_{\boldsymbol{\delta}}\\\
\mathbf{0}_{k_{q}\times p}&\mathbf{0}_{k_{q}\times
k_{z}}\end{pmatrix}=O(T^{-1}).$ (50)

Under our set of assumptions, Lemma 3 and Cramer-Wold device suggest that
$\frac{1}{\sqrt{N}}\frac{\partial\mathbf{g}_{N}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}^{{}^{\prime}}}\xrightarrow{d}N\left(\mathbf{0},\,\operatorname{plim}_{n\to\infty}\boldsymbol{\Omega}_{N}\right)$.
Then, the asymptotic distribution of the GMM estimator in our setting follows
from the asymptotic argument in Lee, (2007, Proposition 2).

##  Appendix D Proof of Theorem 2

When $T$ is large, the precision matrix of
$\sqrt{N}\left(\hat{\boldsymbol{\theta}}_{N}-\boldsymbol{\theta}_{0}\right)$
reduces to

$\displaystyle\mathbf{D}^{{}^{\prime}}_{1N}\boldsymbol{\Omega}^{-1}_{N}\mathbf{D}_{1N}$
$\displaystyle=\frac{1}{N}\begin{pmatrix}\mathbf{C}^{{}^{\prime}}_{N}\left(\boldsymbol{\Delta}_{mN}+\frac{\mu_{4}-3\sigma^{4}_{0}}{\sigma^{4}_{0}}\boldsymbol{\omega}^{{}^{\prime}}_{mN}\boldsymbol{\omega}_{mN}\right)^{-1}\mathbf{C}_{N}&\mathbf{0}_{p\times
k_{z}}\\\ \mathbf{0}_{k_{z}\times p}&\mathbf{0}_{k_{z}\times
k_{z}}\end{pmatrix}$
$\displaystyle+\frac{1}{N\sigma^{2}_{0}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right).$
(51)

As shown by Lee and Yu, (2014), the maximum of
$\mathbf{C}^{{}^{\prime}}_{N}\left(\boldsymbol{\Delta}_{mN}+\frac{\mu_{4}-3\sigma^{4}_{0}}{\sigma^{4}_{0}}\boldsymbol{\omega}^{{}^{\prime}}_{mN}\boldsymbol{\omega}_{mN}\right)^{-1}\mathbf{C}_{N}$
is obtained by choosing $\mathbf{P}^{*}_{j}$ for $j=1,2,\ldots,p$. Thus, the
maximum is

$\displaystyle\mathbf{C}^{*}_{N}=\begin{pmatrix}\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{*s}_{1N}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{*s}_{1N}\mathbf{J}_{N}\right)\\\
\vdots&\ddots&\vdots\\\
\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{1N}\mathbf{J}_{N}\mathbf{P}^{*s}_{pN}\mathbf{J}_{N}\right)&\ldots&\mathrm{tr}\left(\mathbf{G}^{{}^{\prime}}_{pN}\mathbf{J}_{N}\mathbf{P}^{*s}_{pN}\mathbf{J}_{N}\right)\end{pmatrix}.$
(52)

Next , we consider the second term on the right hand side of (D). It follows
from Lemma 1 (4) that
$\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\left(\mathbf{L}_{N},\mathbf{Z}^{**}_{N}\right)=\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\mathbf{Q}^{{}^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}^{*}_{N}$,
where
$\mathbf{Q}^{*}_{N}=(\mathbf{Q}^{*^{\prime}}_{1},\ldots,\mathbf{Q}^{*^{\prime}}_{T-1})^{{}^{\prime}}$.
Thus, we have

$\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{M}_{\mathbf{Q}}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)\leq\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\mathbf{Q}^{*^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}^{*}_{N},$

suggesting that $\mathbf{Q}^{*}_{N}$ is the best IV matrix. From Lemma 5 of
Lee and Yu, (2014), it also follows that

$\displaystyle\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)^{{}^{\prime}}\mathbf{J}_{N}\left(\mathbf{L}_{N},\,\mathbf{Z}_{N}\right)=\operatorname{plim}_{n,T\to\infty}\frac{1}{N}\mathbf{Q}^{*^{\prime}}_{N}\mathbf{J}_{N}\mathbf{Q}^{*}_{N}.$
(53)

Thus, the results in (52) and (53) suggest that that the precision matrix of
$\sqrt{N}\left(\hat{\boldsymbol{\theta}}^{*}_{N}-\boldsymbol{\theta}_{0}\right)$
is given by $\boldsymbol{\Sigma}^{*}_{N}$. Finally, the asymptotic
distribution of
$\sqrt{N}\left(\hat{\boldsymbol{\theta}}^{*}_{N}-\boldsymbol{\theta}_{0}\right)$
can be shown by following the argument given in Theorem 2 of Lee and Yu,
(2014). Therefore, the details are omitted.

## References

  * Abadir and Magnus, (2005) Abadir, K. M. and Magnus, J. R. (2005).  Matrix Algebra.  Cambridge University Press, New York.
  * Bashar, (2021) Bashar, O. H. (2021).  An intra-city analysis of house price convergence and spatial dependence.  The Journal of Real Estate Finance and Economics, 63(4):525–546.
  * Billé et al., (2017) Billé, A. G., Benedetti, R., and Postiglione, P. (2017).  A two-step approach to account for unobserved spatial heterogeneity.  Spatial Economic Analysis, 12(4):452–471.
  * Bollerslev et al., (1992) Bollerslev, T., Chou, R. Y., and Kroner, K. F. (1992).  Arch modeling in finance: A review of the theory and empirical evidence.  Journal of Econometrics, 52(1):5–59.
  * Chang and Diao, (2021) Chang, Z. and Diao, M. (2021).  Inter-city transport infrastructure and intra-city housing markets: Estimating the redistribution effect of high-speed rail in Shenzhen, China.  Urban Studies.
  * Engle, (1982) Engle, R. F. (1982).  Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation.  Econometrica, 50(4):987–1007.
  * Engle and Bollerslev, (1986) Engle, R. F. and Bollerslev, T. (1986).  Modelling the persistence of conditional variances.  Econometric Reviews, 5(1):1–50.
  * Gupta and Robinson, (2015) Gupta, A. and Robinson, P. M. (2015).  Inference on higher-order spatial autoregressive models with increasingly many parameters.  Journal of Econometrics, 186(1):19–31.
  * Hølleland and Karlsen, (2020) Hølleland, S. and Karlsen, H. A. (2020).  A stationary spatio-temporal GARCH model.  Journal of Time Series Analysis, 41(2):177–209.
  * Holmes et al., (2017) Holmes, M. J., Otero, J., and Panagiotidis, T. (2017).  A pair-wise analysis of intra-city price convergence within the Paris housing market.  The Journal of Real Estate Finance and Economics, 54(1):1–16.
  * Jacquier et al., (1994) Jacquier, E., Polson, N. G., and Rossi, P. E. (1994).  Bayesian analysis of stochastic volatility models.  Journal of Business & Economic Statistics, 12(4):371–389.
  * Kelejian and Prucha, (2010) Kelejian, H. H. and Prucha, I. R. (2010).  Specification and estimation of spatial autoregressive models with autoregressive and heteroskedastic disturbances.  Journal of Econometrics, 157:53–67.
  * Kim et al., (1998) Kim, S., Shephard, N., and Chib, S. (1998).  Stochastic volatility: Likelihood inference and comparison with ARCH models.  The Review of Economic Studies, 65(3):361–393.
  * Lee, (2004) Lee, L.-f. (2004).  Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models.  Econometrica, 72(6):1899–1925.
  * Lee, (2007) Lee, L.-f. (2007).  GMM and 2SLS estimation of mixed regressive, spatial autoregressive models.  Journal of Econometrics, 137(2):489–514.
  * Lee and Liu, (2010) Lee, L.-f. and Liu, X. (2010).  Efficient GMM estimation of high order spatial autoregressive models with autoregressive disturbances.  Econometric Theory, 26(1):187–230.
  * Lee and Yu, (2010) Lee, L.-f. and Yu, J. (2010).  A spatial dynamic panel data model with both time and individual fixed effects.  Econometric Theory, 26:564–597.
  * Lee and Yu, (2014) Lee, L.-f. L. and Yu, J. (2014).  Efficient GMM estimation of spatial dynamic panel data models with fixed effects.  Journal of Econometrics, 180(2):174–197.
  * LeSage and Pace, (2009) LeSage, J. P. and Pace, R. K. (2009).  Introduction to Spatial Econometrics (Statistics: A Series of Textbooks and Monographs.  Chapman and Hall/CRC, London.
  * McMillen, (2014) McMillen, D. (2014).  Local quantile house price indices.  Journal of Urban Economics.
  * Meen, (1999) Meen, G. (1999).  Regional house prices and the ripple effect: A new interpretation.  Housing Studies, 14(6):733–753.
  * Merk and Otto, (2021) Merk, M. S. and Otto, P. (2021).  Directional spatial autoregressive dependence in the conditional first-and second-order moments.  Spatial Statistics, 41:100490.
  * Otto, (2019) Otto, P. (2019).  spGARCH: An R-Package for Spatial and Spatiotemporal ARCH and GARCH models.  The R Journal, 11(2):401–420.
  * Otto and Schmid, (2020) Otto, P. and Schmid, W. (2020).  Spatial and spatiotemporal GARCH models – a unified approach.
  * Otto et al., (2018) Otto, P., Schmid, W., and Garthoff, R. (2018).  Generalised spatial and spatiotemporal autoregressive conditional heteroscedasticity.  Spatial Statistics, 26:125–145.
  * Robinson, (2009) Robinson, P. M. (2009).  Large-sample inference on spatial dependence.  Econometrics Journal, 12.
  * Sandmann and Koopman, (1998) Sandmann, G. and Koopman, S. J. (1998).  Estimation of stochastic volatility models via Monte Carlo maximum likelihood.  Journal of Econometrics, 87(2):271 – 301.
  * Sato and Matsuda, (2017) Sato, T. and Matsuda, Y. (2017).  Spatial autoregressive conditional heteroskedasticity models.  Journal of the Japan Statistical Society, 47(2):221–236.
  * Sato and Matsuda, (2021) Sato, T. and Matsuda, Y. (2021).  Spatial extension of generalized autoregressive conditional heteroskedasticity models.  Spatial Economic Analysis, 16(2):148–160.
  * Shephard, (1994) Shephard, N. (1994).  Partial non-Gaussian state space.  Biometrika, 81(1):115–131.
  * Taşpınar et al., (2021) Taşpınar, S., Doğan, O., Chae, J., and Bera, A. K. (2021).  Bayesian inference in spatial stochastic volatility models: An application to house price returns in Chicago.  Oxford Bulletin of Economics and Statistics, 83:1243–1272.
  * Tobler, (1970) Tobler, W. R. (1970).  A computer movie simulating urban growth in the Detroit region.  Economic Geography, 46:234–240.
  * Yu et al., (2008) Yu, J., de Jong, R., and fei Lee, L. (2008).  Quasi-maximum likelihood estimators for spatial dynamic panel data with fixed effects when both n and t are large.  Journal of Econometrics, 146(1):118 – 134.
  * Zhang and Yi, (2017) Zhang, L. and Yi, Y. (2017).  Quantile house price indices in Beijing.  Regional Science and Urban Economics, 63:85–96.