Dynamic Spatiotemporal ARCH Models

Dynamic Spatiotemporal ARCH Models

\def\spacingset#1{ {#1}} \spacingset{1}

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

\spacingset{1.45}

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:1992, Engle: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 Tobler70. From a statistical perspective, this spatial dependence may occur in the means and the volatility of a random process Otto: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. gupta2015inference,merk2021directional). 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,lee2010efficient,Lee: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.\footnote{In the literature, the log-square transformation approach is also used for the estimation of (i) the cross-sectional spatial stochastic volatility models Taspinar:2021, and (ii) the cross-sectional spatial ARCH/GARCH models Sato:2017,Otto:2019,Otto:2020,Takaki:2021.} In a similar manner, we may alternatively consider the QML estimation approach rather than the GMM approach for our transformed model Yu:2008,Yu:2010,Hol: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 large\footnote{Lesage: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:2008,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 JPR:1994, Shephard:1994, Kim:1998, 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 mcmillen2014local,bille2017two,zhang2017quantile). In Section (ref), 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 (ref) 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 chang2021inter. Recently, holmes2017pair and bashar2021intra 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.

The rest of the paper proceeds in the following way. In Section (ref), we state our model specification and discuss its properties. In Section (ref), we provide the details of the GMM estimation approach for our model, and formally establish its large sample properties. In Section (ref), we investigate the finite sample properties of our suggested algorithm through an extensive Monte Carlo study. In Section (ref), we provide the details of our empirical application on Berlin's house price returns. In Section (ref), we offer our concluding comments with some directions for future studies. Some technical results are relegated to \hyperref[app]{Appendix}.

Model Specification

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

for $i=1,2,\hdots,n$ and $t=1,\hdots 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 (ref), where $\{m_{l,ij}\}_{l=1}^p$, for $i,j=1,\hdots,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,\hdots,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 (ref), $\mathbf{x}_i$ is a $k\times1$ 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},\hdots,\mu_{n0})^{'}$ and $\boldsymbol{\alpha}_0=(\alpha_{10},\hdots,\alpha_{T0})^{'}$. 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},\hdots,y_{n0})^{'}$ is observable.

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

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 (ref) and (ref) as

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

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

The Estimation Approach

The elements of $\boldsymbol{\operatorname{\varepsilon}}^{*}_t$ in (ref) 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.

where $\mathbf{U}_t=(u_{1t},\hdots,u_{nt})^{'}=\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}^{'}_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,\hdots,\mathbf{c}_T\right)$ be an $n\times T$ matrix, where $\mathbf{c}_t$ is an $n\times1$ vector for $t=1,\hdots,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,\hdots,\mathbf{c}^{*}_{T-1}\right)=\left(\mathbf{c}_1,\hdots,\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,\hdots,T-1$. If we apply $\mathbf{F}_{T,T-1}$ to our model in (ref) in a similar manner, we obtain

for $t=1,\hdots,T-1$, where $\left(\mathbf{Y}^{**}_1,\hdots,\mathbf{Y}^{**}_{T-1}\right)=\left(\mathbf{Y}^{*}_1,\hdots,\mathbf{Y}^{*}_{T}\right)\mathbf{F}_{T,T-1}$, $\left(\mathbf{Y}^{**,-1}_0,\hdots,\mathbf{Y}^{**,-1}_{T-2}\right)=\left(\mathbf{Y}^{*}_0,\hdots,\mathbf{Y}^{*}_{T-1}\right)\mathbf{F}_{T,T-1}$, $\left(\mathbf{X}^{*}_{l1},\hdots,\mathbf{X}^{*}_{k,T-1}\right)=\left(\mathbf{X}_{l1},\hdots,\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,\hdots,k$, $\left(\alpha^{*}_1,\hdots,\alpha^{*}_{T-1}\right)=\left(\alpha_1,\hdots,\alpha_T\right)\mathbf{F}_{T,T-1}$ and $\left(\mathbf{U}^{*}_{1},\hdots,\mathbf{U}^{*}_{T-1}\right)=\left(\mathbf{U}_{1},\hdots,\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

Let $N=n(T-1)$ and $\mathbf{U}_N=(\mathbf{U}^{*'}_1,\hdots,\mathbf{U}^{*'}_{T-1})^{'}$. Note that we can express $\left(\mathbf{U}^{*}_{1},\hdots,\mathbf{U}^{*}_{T-1}\right)=\left(\mathbf{U}_{1},\hdots,\mathbf{U}_{T}\right)\mathbf{F}_{T,T-1}$ as $\left(\mathbf{U}^{*'}_{1},\hdots,\mathbf{U}^{*'}_{T-1}\right)^{'}=\left(\mathbf{F}^{'}_{T,T-1}\otimes\mathbf{I}_n\right)\left(\mathbf{U}^{'}_{1},\hdots,\mathbf{U}^{'}_{T}\right)^{'}$.\footnote{Note 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}^{'}\otimes \mathbf{A})\text{vec}(\mathbf{B})$, where $\text{vec}(\mathbf{B})$ denotes the vectorization of the matrix $\mathbf{B}$ Karim:2005. This property can be applied to $\left(\mathbf{U}^*_1,\mathbf{U}^{*}_2,\hdots,\mathbf{U}^*_{T-1}\right)=\left(\mathbf{U}_1,\mathbf{U}_2,\hdots,\mathbf{U}_T\right)\mathbf{F}_{T,T-1}$ by setting $\mathbf{D}=\left(\mathbf{U}^*_1,\mathbf{U}^{*}_2,\hdots,\mathbf{U}^*_{T-1}\right)$, $\mathbf{C}=\mathbf{F}_{T,T-1}$, $\mathbf{B}=\left(\mathbf{U}_1,\mathbf{U}_2,\hdots,\mathbf{U}_T\right)$ and $\mathbf{A}=\mathbf{I}_n$. } Then, it follows that

indicating that the elements of $\mathbf{U}_N$ are uncorrelated. Among the orthonormal transformations, Lee: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 (ref) 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}^{'}_n\right)$ to get

where we used the fact that $\mathbf{J}_n\mathbf{1}_n=\mathbf{0}_{n}$. Our GMM estimation approach is based on (ref). 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,\hdots,T-1$. That is, we need IVs for the following variables:

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

We use $\boldsymbol{\rho}_0=\left(\rho_{10},\hdots,\rho_{p0}\right)^{'}$ and $\boldsymbol{\delta}_0=\left(\delta_{10},\hdots,\delta_{p0}\right)^{'}$ to denote the true parameter values, and $\boldsymbol{\rho}=\left(\rho_{1},\hdots,\rho_{p}\right)^{'}$ and $\boldsymbol{\delta}=\left(\delta_{1},\hdots,\delta_{p}\right)^{'}$ 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 (ref) can be expressed as

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}^{**'}_1,\hdots,\mathbf{Z}^{**'}_{T-1}\right)^{'}$. Then, using (ref), we have

where $\boldsymbol{\eta}_0=(\gamma_0,\boldsymbol{\delta}^{'}_0,\boldsymbol{\beta}^{'}_0)^{'}$, and $\mathbf{G}_{r}=\mathbf{M}_r\mathbf{S}^{-1}$. We can use (ref) 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,\hdots,T-1$, and the time lag variables $\mathbf{Y}^{*}_{0},\hdots,\mathbf{Y}^{*}_{t-1}$ as IVs. Similarly, we can use $\mathbf{M}_j\mathbf{X}^{*}_s$ for $s=1,\hdots,T-1$, and $\mathbf{M}_j\mathbf{Y}^{*}_s$ for $s=0,1,\hdots,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,\hdots,T-1$, where $k_q\geq k+2p+1$. For example, we may choose $\mathbf{Q}_t$ as

where $\mathbb{M}^2\mathbf{Y}^{*}_{t-1}=\left(\mathbf{M}^2_1\mathbf{Y}^{*}_{t-1},\hdots,\mathbf{M}_1\mathbf{M}_p\mathbf{Y}^{*}_{t-1},\mathbf{M}_2\mathbf{M}_1\mathbf{Y}^{*}_{t-1},\hdots,\mathbf{M}_1\mathbf{M}_p\mathbf{Y}^{*}_{t-1},\hdots,\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}^{'}_1,\hdots,\mathbf{Q}^{'}_{T-1})^{'}$, $\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

where $\boldsymbol{\theta}=\left(\boldsymbol{\rho}^{'},\boldsymbol{\eta}^{'}\right)^{'}$, $\mathbf{U}_N(\boldsymbol{\theta})=\left(\mathbf{U}^{*'}_1(\boldsymbol{\theta}),\hdots,\mathbf{U}^{*'}_{T-1}(\boldsymbol{\theta})\right)^{'}$ 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},\hdots,\alpha^{*}_{T-1})^{'}$ 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: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 (ref). 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

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

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}^{'}$. 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 (ref), it can be shown that\footnote{In applying Lemma (ref), we use the fact that $\mathrm{tr}\left(\mathbf{A}^{'}\mathbf{B}\right)=\text{vec}^{'}\left(\mathbf{A}\right)\text{vec}\left(\mathbf{B}\right)=\text{vec}^{'}\left(\mathbf{B}\right)\text{vec}\left(\mathbf{A}\right)$, where $\mathbf{A}$ and $\mathbf{B}$ are any two $N\times N$ matrices.}

where $\mu_4$ is the fourth moment of $u_{it}$, $\boldsymbol{\omega}_{mN}=\left(\text{vec}_D\left(\mathbf{J}_N\mathbf{P}^{'}_{1N}\mathbf{J}_N\right),\hdots,\text{vec}_D\left(\mathbf{J}_N\mathbf{P}^{'}_{mN}\mathbf{J}_N\right)\right)$ and

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

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

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

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

where $\tau_1=\left(\sum_{l=1}^{p}|\rho_{l}|\right)\times\max_{1\leq l\leq p}\left\Vert \mathbf{M}_{l}\right\Vert<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}\Vert\mathbf{M}_{l}\Vert\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 (ref) 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 (ref)). The full column rank condition in Assumption (ref) gives the identification condition based on the linear moment function in our setting. See Appendix (ref) for the details on the identification condition in our setting.

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

Our estimator defined in (ref) 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}^{**'}_1,\hdots,\mathbf{Y}^{**'}_{T-1}\right)^{'}$ and $\mathbf{Z}_N=\left(\mathbf{Z}^{**'}_1,\hdots,\mathbf{Z}^{**'}_{T-1}\right)^{'}$. Then, the 2SLS estimator is

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

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,\hdots,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,\hdots,p$. Then, the initial GMM estimator is given by $\tilde{\boldsymbol{\theta}}_N=\operatorname{argmin}_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\mathbf{g}^{'}_N(\boldsymbol{\theta})\mathbf{g}_N(\boldsymbol{\theta})$ where

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}}^{'}_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: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 (ref) depends on the IV matrix $\mathbf{Q}_N$ and the quadratic moment matrices $\mathbf{P}_{lN}$ for $l=1,2,\hdots,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 (ref) indicates that $\frac{1}{N}\frac{\partial \mathbf{g}_N(\boldsymbol{\theta}_0)}{\partial\boldsymbol{\theta}^{'}}=\mathbf{D}_{1N}+O_p(N^{-1/2})$, where

Then, (ref) shows that the precision matrix of $\sqrt{N}\left(\hat{\boldsymbol{\theta}}_N-\boldsymbol{\theta}_0\right)$ is

Since the above precision matrix has the same form as the one given in Lee: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}^{'}_N\left(\boldsymbol{\Delta}_{mN}+\frac{\mu_4-3\sigma^4_0}{\sigma^4_0}\boldsymbol{\omega}^{'}_{mN}\boldsymbol{\omega}_{mN}\right)^{-1}\mathbf{C}_N$. As shown in Lee:2014, these matrices are

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 (ref), $\mathbb{E}\left(\mathbf{Y}^{**}_{t-1}|\mathcal{F}_{t-1}\right)$ can be approximated by

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$.\footnote{Note 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 (ref), 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,\hdots,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

The feasible version of $\mathbf{Q}^{*}_t$ can be obtained by substituting consistent estimators of the unknown parameters into (ref). 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}^{'}_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}^{'}_n\hat{\boldsymbol{\vartheta}}_t$ for $t=1,2,\hdots,T$.\footnote{Note 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}^{'}_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,\hdots,p$.

A Monte Carlo Study

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

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\times1},\mathbf{I}_2)$ with $\boldsymbol{\beta}_0=(0.5,1)^{'}$. 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)^{'}=\{(0.2,0.2, -0.2)^{'},\,(0.2,0.8,-0.2)^{'}\}$ 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)^{'}$. 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 (ref) - (ref). 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.

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?

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 (ref), 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 (ref). 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 (ref)) and model the temporal effects by including yearly and monthly indicator variables as regressors.

We estimate a first-order version of our model in (ref). We specify the spatial weights matrix (row-standardized) based on the queen contiguity scheme, where all adjacent neighbors are equally weighted. In Table (ref) and Figure (ref), 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 (ref)). 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 (ref), 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 (ref). 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.

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.