Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
54,098 characters · 11 sections · 34 citation commands
Generalized Spatial and Spatiotemporal ARCH Models
\def\spacingset#1{ {#1}} \spacingset{1}
{\it Keywords:} Spatial GARCH, spatiotemporal statistics, unified approach, variance clusters, real estate prices.
\spacingset{1.45}
Recent literature have dealt with the extension of generalized autoregressive conditional heteroscedasticity (GARCH) models to spatial and spatiotemporal processes (e.g. Otto16_arxiv,Otto18_spARCH,Otto19_statpapers,Sato17,Sato18b,Sato18a). Whereas the classical ARCH model is defined as a process over time, these random processes have a multidimensional support. Thus, they allow for spatially dependent second-order moments, while the local means are uncorrelated and constant in space (see Otto19_statpapers). Sato17,sato2020spatial introduced a random process incorporating elements of GARCH and exponential GARCH (E-GARCH) processes, which is, however, neither a GARCH nor an E-GARCH process. Moreover, Otto16_arxiv only focussed on spatial ARCH processes without considering the influences from the realized, conditional variance at neighbouring locations. Direct extensions of GARCH and E-GARCH processes to spatial settings do not exist among current research.
{In this paper, we introduce a completely novel generalized spatial ARCH model (spGARCH). Because a general definition of this model is used, time-series GARCH models (Bollerslev86), the previously introduced spatial ARCH (Otto16_arxiv) and the hybrid GARCH (Sato17,Sato18a) are included. This definition also allows us to define a spatial logarithmic and exponential spatial GARCH model, which will be the subject of a future paper. Moreover,} other GARCH-type models, like threshold or multivariate GARCH models, can easily be constructed. This unified spatial GARCH process is a completely new class of models in spatial econometrics, {for which we derive consistent estimators based on a non-linear least-squares approach.} In addition, all models are computationally implemented in one library, the R-package spGARCH {(version $> 2.0$)}.
From a practical perspective, this unified spatial GARCH model can be used to model spill-over effects in the conditional variances {across the spatial units}. That means that an increasing variance in a certain region of the considered space would lead to an increase or decrease in the adjacent regions, depending on the direction (sign) of the spatial dependence. {Compared to previous spatiotemporal GARCH models, these spill-overs are instantaneous.} Local climate risks, such as fluctuations in the temperature and precipitation, or financial risks in spatially constrained markets, such as real estate or labour, could be modelled using this approach. Furthermore, spatial GARCH-type models can be used as error models for any linear or non-linear spatial regression model to account for local model uncertainties (i.e., areas in which the considered models perform worse than in others). Such model uncertainties can be considered to be a kind of local risk.
The remainder of this paper is structured as follows. In the next section, we introduce the {generalized framework of spatial and spatiotemporal autoregressive conditional heteroscedasticity models} and discuss two examples nested within this approach, more precisely, the novel spatial GARCH (as an equivalent to the time-series GARCH models) and the hybrid spatial GARCH processes by Sato17,sato2020spatial. Following from there, a non-linear least-squares procedure is introduced for this model class. These theoretical sections are followed up with a discussion of the insights gained from simulation studies. The paper then supplies a real-world example, namely the real estate prices in the German capital city of Berlin. In Section (ref), we stress some important extensions for future research before concluding the paper.
Let $\left\{Y(\boldsymbol{s}) \in \mathds{R}: \boldsymbol{s} \in D_{\boldsymbol{s}} \right\}$ be a univariate stochastic process, where $D_{\boldsymbol{s}}$ represents a set of possible locations in a $q$-dimensional space. Thus, spatial and spatiotemporal models are both covered by this approach. With regards to spatiotemporal processes, the temporal dimension can be easily considered as one of the $q$ dimensions. In addition, time-series GARCH models are included for $q = 1$.
Let $\boldsymbol{s}_1, \ldots, \boldsymbol{s}_n$ denote all locations, and let $\boldsymbol{Y}$ stand for the vector of observations $\left(Y\left(\boldsymbol{s}_i\right)\right)_{i = 1, \ldots, n}$. The commonly applied spatial autoregressive (SAR) model implies that the conditional variance $Var(Y(\boldsymbol{s}_i) | Y(\boldsymbol{s}_j), j \neq i)$ is constant (cf. Cressie93,Cressie11) and does not depend on the observations of neighbouring locations. This approach is extended by assuming the changes in the volatility can spill over to neighbouring regions and that conditional variances can vary over space, resulting in clusters of high and low variance. As in time-series ARCH models developed by Engle82, the vector of observations is given by the non-linear relationship
where $\boldsymbol{h} = (h(\boldsymbol{s}_1), \ldots, h(\boldsymbol{s}_n))'$ and $\boldsymbol{\varepsilon} = (\varepsilon(\boldsymbol{s}_1), \ldots, \varepsilon(\boldsymbol{s}_n))'$ is a noise component, which is later specified in more detail.
{Moreover}, we assume that a known function $f$ exists, which relates $\boldsymbol{h}$ to a vector $\boldsymbol{F} = (f(h(\boldsymbol{s}_1)), \ldots, f(h(\boldsymbol{s}_n)))^\prime$. {This general approach is beneficial because} different spatial GARCH-type processes can be defined by choosing $f$ and a suitable model of $\boldsymbol{F}$. For instance, they could have additive or multiplicative dynamics, or the spill-over effects in the conditional variances could be global or locally constrained to direct neighbouring observations. {In this paper, we initially focus on generalized spatial GARCH models (spGARCH), which are analogously defined to the time-series GARCH models and have additive dynamics (cf. Bollerslev86). Besides, previously introduced spatial ARCH models are nested within this general approach} (e.g., Otto18_spARCH,sato2020spatial; see also Examples (ref) and (ref) in Section (ref)).
{Below,} we introduce a general approach covering some important spatial and spatiotemporal GARCH-type models, namely the spatial ARCH model of Otto16_arxiv,Otto18_spARCH and the hybrid model of sato2020spatial. For these models, vector $\boldsymbol{F}$ is chosen as
with a measurable function $\boldsymbol{\gamma}(\boldsymbol{x}) = (\gamma_1(\boldsymbol{x}), \ldots, \gamma_n(\boldsymbol{x}))^\prime$ and $\boldsymbol{Y}^{(2)} = (Y(\boldsymbol{s}_1)^2, \ldots , Y(\boldsymbol{s}_n)^2)^\prime$. The weighting matrices $\mathbf{W}_1 = (w_{1,ij})_{i,j = 1, \ldots, n}$ and $\mathbf{W}_2 = (w_{2,ij})_{i,j = 1, \ldots, n}$ are assumed to be non-negative with zeros on the diagonal (i.e., $w_{v,ij} \ge 0$ and $w_{v,ii}=0$ for all $i,j = 1, \ldots, n$ and $v = 1, 2$). Moreover, let $\boldsymbol{\alpha} = (\alpha_i)_{i = 1, \ldots, n}$ be a positive vector.
First, we discuss under what conditions the process is well-defined. {To do this, we make use of the Banach fixed point theorem for random processes. The field of random fixed point theorems has been {studied} by several authors (e.g., hanvs1957reduzierende,bharucha1976fixed,tan1997random).
In the following, we make use of the notation $\mathbf{E} = \mbox{diag}(\varepsilon(\boldsymbol{s}_1)^2, \ldots, \varepsilon(\boldsymbol{s}_n)^2)$. {Considering} the operator
defined on $I\!\!R^n$ with a norm $||.||${, the following conditions can be derived such that the process is well-defined.}
{The proof of this theorem is given in the Appendix.} Note that the condition ((ref)) is fulfilled if, for example, $\boldsymbol{\gamma}$ satisfies a Lipschitz condition with constant $L_1$, $(f(z_i))_{i=1,...,n}$ satisfies a Lipschitz condition with a constant $L_2$ and $L_n := L_1 ||\mathbf{W}_1 \mathbf{E}|| + L_2 ||\mathbf{I} - \mathbf{W}_2|| < 1$ where we make use of the matrix norm induced by the vector norm. However, in order to guarantee that $L_n$ does not depend on $n$, we need stronger conditions. If we take the 1-norm and if the matrices $\mathbf{W}_1$ and $\mathbf{W}_2$ are row-standardized then $||\mathbf{I} - \mathbf{W}_2|| < 2$. To ensure that $||\mathbf{W}_1 \mathbf{E}||$ is bounded we have to assume that the $\varepsilon(\boldsymbol{s}_i)$ are uniformly bounded. We refer to Otto18_spARCH where this problem is discussed for a spatial ARCH process in more detail.
{Moreover, it is important to note that} the operator $\boldsymbol{T}(\omega)$ is continuous if $f$ and $\gamma$ are continuous.
Further, the solution of ((ref)) which reflects $\boldsymbol{h}$ should be non-negative such that the process $\boldsymbol{Y}$ is a well-defined real-valued process. In many applications, the functions $f$ and $\boldsymbol{\gamma}$ are defined to be zero for negative values{, such that $\boldsymbol{h}$ is always positive if $\boldsymbol{\alpha} > 0$.} We will come back to this point later. } In addition, the fixed-point theorem of Banach implies that the sequence $\boldsymbol{z}_m = \boldsymbol{T}(\boldsymbol{z}_{m-1})$, $m \ge 1$ converges to $\boldsymbol{h}$ for given $\boldsymbol{\alpha}$, $\mathbf{E}$, $\mathbf{W}_1$, and $\mathbf{W}_2$. Consequently, this result {represents one way} to simulate such a process.
Below, we discuss some important properties of this process including the following condition for stationarity.
Moreover, the observations $Y(\boldsymbol{s})$ are uncorrelated with a mean of zero, as we will show in the following theorem. Thus, spatial GARCH models are suitable error models for use with other linear or non-linear spatial regression models, such as spatial autoregressive or spatial error models (see also Elhorst10), without affecting the mean equation. In this way, locally varying model uncertainties can be captured.
In the spatial setting, however, the conditional variance $Var(Y(\boldsymbol{s}_i) | Y(\boldsymbol{s}_j), j \neq i)$ is not exactly equal to $h(\boldsymbol{s}_i)$ (see Otto19_statpapers). Nevertheless, the interpretation of $\boldsymbol{h}$ is similar to the conditional variance. In locations $\boldsymbol{s}$, where $h(\boldsymbol{s})$ is large, the conditional variance is also large and vice versa (see Otto19_statpapers, Fig. 1). That means that the local risk or level of uncertainty of this particular region is high compared to its neighbours. Such regions could be identified via $\boldsymbol{h}$; this could be of interest in terms of the valuation of real estate or other immovable assets since it provides insights into an individual location's risk.
In addition to this, the spatial GARCH coefficients measure potential risk spill-overs from neighbouring locations. It is worth noting that in the case of directional spatial processes, $\boldsymbol{h}$ is equal to the conditional variances. Thus, it can be interpreted in the same way as with time-series GARCH models (see Otto19_statpapers).
This {general framework} allows for a large range of GARCH-type models. Depending on the definition of $f$ and $\boldsymbol{\gamma}$, the resulting spatial GARCH-type models have different stochastic properties. We discuss some important special cases below, starting with the spatial ARCH model (Otto16_arxiv,Otto18_spARCH,Otto19_statpapers), which is a direct extension of the ARCH process of Engle82 to spatial and spatiotemporal processes. It was originally introduced by Otto16_arxiv. For more details on its stochastic properties, we refer to Otto19_statpapers.
{The process is well-defined if $||\mathbf{W}_1 \mathbf{E}|| < 1$. This is an immediate consequence of Theorem (ref).} Indeed, the spatial ARCH process can be easily extended to a spatial GARCH process by considering the realized values of $h(\cdot)$ in adjacent locations. {This novel spatial GARCH process is defined in the following example.}
Since $\boldsymbol{Y}^{(2)} = \mathbf{E} \boldsymbol{h}$, the quantity $\boldsymbol{h}$ can be specified as
if the inverse exists. For this simple example, there is a unique solution if {$||\mathbf{W}_1 \mathbf{E} + \mathbf{W}_2|| < 1$}, as it is already expressed in Theorem (ref). {Alternatively, the condition is fulfilled if the process is directional. In this case, $\mathbf{W}_1$ and $\mathbf{W}_2$ are lower or upper triangular matrices (cf. Basak18,otto2019stochastic,merk2021directional).}
Contrary to this approach, Sato17,Sato18b,sato2020spatial have considered a slightly different choice of $\boldsymbol{h}$, and have used the log-transformation to avoid any non-negativity problems of $\boldsymbol{h}$. Thus, their model combines the GARCH and the E-GARCH attempts. For that reason, we will {refer to} it as the hybrid model {(H-spGARCH)}. Let $\boldsymbol{h}_L = ( \log(h(\boldsymbol{s}_i)) )_{i=1,...,n}$ and $\boldsymbol{Y}^{(2)}_L = ( \log(Y(\boldsymbol{s}_j)^2 ) )_{i=1,...,n}$.
The process has a unique solution if {$||\mathbf{W}_1 + \mathbf{W}_2|| < 1$.} This is {also} an immediate consequence of Theorem (ref). It is obtained by setting $\gamma(\boldsymbol{x}) = (\log(x_i) )_{i=1,\ldots,n}$, $f(x)=\log(x)$, and considering the right side of (ref) to be a function of $\log(z)$. We see that the condition on the existence of a solution is much simpler than for the spGARCH process, since it only depends on the weight matrices and not on the random matrix $\mathbf{E}$. This simplification is due to the fact that we have an additive decomposition of the function $\boldsymbol{\gamma}$, i.e., $\boldsymbol{\gamma}(\mathbf{E} \boldsymbol{z}) = \boldsymbol{\phi}_1(\mathbf{E}) + \boldsymbol{\phi}_2(\boldsymbol{z})$ with certain functions $\boldsymbol{\phi}_1$ and $\boldsymbol{\phi}_2$. This functional equation is solved by the logarithm function. However, the behavior of the H-spGARCH is different to that of the spGARCH. Thus, the one or the other could be preferable for empirical applications.
In the following section, we firstly discuss the choice of the weight matrices in more detail. In a general setting, $\mathbf{W}_1$ and $\mathbf{W}_2$ have $n(n-1)$ free parameters, while only $n$ values are observed. In spatial econometrics, these matrices are therefore usually replaced with a parametric model to control the influence of adjacent regions. Alternatively, they might instead be estimated using statistical learning approaches, e.g., lasso-type estimators under the assumption of a certain degree of sparsity. {In this section of the paper, however, we will refer to a classical parametric model. For this, we develop an estimation method based on non-linear least squares estimators and show the consistency of these estimators.}
There is great flexibility in the choice of the weight matrices (see Getis09 for an overview). In practice, these are usually dependent upon additional parameters and spatial locations. Frequently, it is assumed that $\mathbf{W}_1 = \rho \mathbf{W}^{*}_1$ and $\mathbf{W}_2 = \lambda \mathbf{W}^{*}_2$ with the predefined, known matrices $\mathbf{W}^{*}_1$ and $\mathbf{W}^{*}_2$. That is, $\mathbf{W}^{*}_1$ and $\mathbf{W}^{*}_2$ describe the structure of the spatial dependence, with the weights as a multiple of these specific matrices. In settings such as these, it is easy to test whether a random process exhibits such a spatial dependence, by testing the parameters $\rho$ and $\lambda$. As with time-series GARCH models, $\rho$ measures the extent to which a volatility shock in one region spills over to neighbouring regions, while $\rho + \lambda$ gives an impression how fast this effect will fade out in space (see, e.g., Campbell97). A more general approach can be obtained by choosing $\mathbf{W}_{k} = \text{diag}(\rho_1, \ldots, \rho_1, \ldots, \rho_k, \ldots, \rho_k) {\mathbf{W}_{\cdot}^{*}}$ as the weights for $k \in \{1,2\}$. Here, different areas are weighted in different ways. For instance, all counties of state $i$ are weighted by $\rho_i$, while counties of another state, $j$, get a different weighting factor, $\rho_j$. Alternatively, $\mathbf{W}_{k}^{*}$ could be chosen as $(K_\theta(\boldsymbol{s}_i - \boldsymbol{s}_j))_{i,j = 1, \ldots, n}$ for $k \in \{1,2\}$ with a known function $K$. In this case, the spatial correlation depends on the distance between two locations. For instance, inverse distance weighting schemes $K(\boldsymbol{x}) = ||\boldsymbol{x}||^{-k}$ with $k$ being estimated, or anisotropic weighting schemes dependent upon the bearing between two locations.
Below, we assume that the weight matrices have the structure
Thus, the model has three parameters to be estimated, $\boldsymbol{\vartheta} = (\rho, \lambda, \alpha)^\prime$. Let $\boldsymbol{\vartheta}_0$ denote the true parameters. In the following we use the symbol $||\boldsymbol{x}||_2$ for the Euclidean norm of a vector $\boldsymbol{x}$ and $||\mathbf{A}||_\infty = max_{1 \le i \le n} \sum_{j=1}^{n} |a_{ij}|$ for the matrix norm of an $n \times n$ matrix $\mathbf{A}$, which is induced by a maximum norm.
One possible method for estimating the parameters is the non-linear least-squares approach (NLSE). Squaring the components of (ref) and taking the logarithms, we get that for $i=1, \ldots, n$
with $\eta(\boldsymbol{s}_i) = \mbox{log}(\varepsilon(\boldsymbol{s}_i)^2) - E( \mbox{log}(\varepsilon(\boldsymbol{s}_i)^2) )$. Now $\eta(\boldsymbol{s}_i), i=1,...,n$ is a white noise process. Moreover, it follows with $\tau(x) = f(\mbox{exp}(x))$ that
where $ \tilde{\boldsymbol{\gamma}}(\boldsymbol{x}) = ( \gamma_i(\mbox{exp}(x_1),..., \mbox{exp}(x_n)) )_{i=1,...,n}$. Note that $(\mathbf{I} - \lambda \mathbf{W}_2^*)^{-1}$ exists if $||\lambda \mathbf{W}_2^*|| < 1$.
Now, let $( c_i(\lambda) )_{i=1,...,n} = (\mathbf{I} - \lambda \mathbf{W}_2^*)^{-1} \boldsymbol{1}$ and $( \boldsymbol{d}_i(\lambda)^\prime )_{i=1,...,n} = (\mathbf{I} - \lambda \mathbf{W}_2^*)^{-1} \mathbf{W}_1^*$. In order to denote the dependence on $\boldsymbol{\vartheta}$ we write $h_{\boldsymbol{\vartheta}}(\boldsymbol{s}_i)$, $i=1,...,n$. Then, \[ \mbox{log}(h_{\boldsymbol{\vartheta}}(\boldsymbol{s}_i)) = \tau^{-1}(\alpha \, c_i(\lambda) + \rho \, \boldsymbol{d}_i(\lambda)^\prime \, \tilde{\boldsymbol{\gamma}}(\mbox{log}(\boldsymbol{Y}^{(2)} ) ) ) . \] Here, we assume that $c = E( \mbox{log}(\varepsilon(\boldsymbol{s}_i)^2) )$ is a known quantity. Using $H_i = \mbox{log}(Y(\boldsymbol{s}_i)^2) - E( \mbox{log}(\varepsilon(\boldsymbol{s}_i)^2) )$ and $\boldsymbol{H} = ( H_i )_{i=1,...,n}$ the estimators of the parameters $\alpha$, $\lambda$, and $\rho$ are obtained by minimizing the non-linear sum of squares \[ \sum_{i=1}^n \left( H_i - \mbox{log}(h_{\boldsymbol{\vartheta}}(\boldsymbol{s}_i) \right)^2 = \sum_{i=1}^n \left( H_i - \tau^{-1}(\alpha \, c_i(\lambda) + \rho \, \boldsymbol{d}_i(\lambda)^\prime \, \tilde{\boldsymbol{\gamma}}(\boldsymbol{H} + c \boldsymbol{1} ) ) \right)^2 \] with respect to $\boldsymbol{\vartheta}$.
Although $\tau^{-1}$ is a known function, this minimization problem is complex. Thus, we will impose further assumptions which are fulfilled for all relevant special cases. We will suppose that $\boldsymbol{\gamma}(\boldsymbol{x}) = ( \gamma( x_i ) )$ with a known function $\gamma$. Consequently, $\tilde{\boldsymbol{\gamma}}(\boldsymbol{x}) = ( \tilde{\gamma}(x_i) )_{i=1,...,n}$ with $\tilde{\gamma}(x) = \gamma( \mbox{exp}(x) )$, which leads to the easier minimization of
Note that since the $(i,i)$-th element of $\mathbf{W}_1^*$ is zero it follows that $\boldsymbol{d}_i(\lambda)^\prime \, \left( \tilde{\gamma}( H_v + c ) \right)_{v=1,...,n}$ is no function of $H_i$. Minimization problems of that type have been studied in detail in, e.g., Amemiya85,Potscher97,Newey94. There, sufficient conditions are given for the {consistency} and asymptotic normality of the resulting estimators under various conditions. Note that, in the present case, $\{ H_i \}$ is a strictly stationary process. Moreover, in most papers on this topic, the regression function is assumed to be a deterministic function depending on certain parameters. In the present case, however, it is a function depending on the observations $\{ H_i \}$ which makes the analysis of the asymptotic behaviour of the estimators much harder. Further, it must be noted that a spatial problem is present. The positions $\boldsymbol{s}_i$ are points in a space and we need a certain distance measure between these points to assess the dependence of the observations.
Note that the solution of ((ref)) does not have to be unique. For more details, we refer to Section 4 of Amemiya85.
Moreover, for a spGARCH process it holds that $f^\prime(x) x = x$ and $h_{\boldsymbol{\vartheta}}(\boldsymbol{s}_i) \ge \alpha_l > 0$ and thus the above condition is fulfilled. For a H-spGARCH process we have that $f^\prime(x) x = 1$ and thus it is fulfilled as well. Further, for a H-spGARCH process the condition ((ref)) can be easily seen to be fulfilled since for a strictly stationary process $\{ Y(\boldsymbol{s}_i ) \}$ the quantity $E(\log(h_{\boldsymbol{\vartheta}}(\boldsymbol{s}_i)))$ does not depend on $i$ at all.
{Moreover, to prove the consistency of the local minimum, the roots of the first derivative of the sum of non-linear squares with respect to the parameters must be zero, i.e.,} \[ \frac{\partial Q_n(\boldsymbol{\vartheta})}{\partial \boldsymbol{\vartheta}} = 0 . \]
{Like in the previous section, we will now consider a special case of the general framework, namely an spGARCH model. That is,} we choose $f(x) = x I_{ [0, \infty ]}(x)$ while $\gamma$ is an arbitrary function satisfying certain conditions.
Note that the assumption ((ref)) is a statement about the topological structure of the underlying space. Moreover, ((ref)) shows that it can be interpreted as an assumption on the underlying autocorrelation structure of the process. It is fulfilled if $\mathbf{W}_1^*$ and $\mathbf{W}_2^*$ are sparse to limit the spatial dependence to a manageable degree. Here the choice of the weight matrices is restricted. It is also satisfied if the autocorrelation is weak.
Here, we assume the simple parametric setting given by (ref). {We simulated a spatial GARCH process as specified in Example (ref) and the} weighting matrices $\mathbf{W}^*_1$ and $\mathbf{W}^*_2$ were set as row-standardized Rook contiguity matrices, for which the upper diagonal elements were set to zero to avoid negative values of $h(\boldsymbol{s}_i)$. Thus, {Theorem (ref) is fulfilled. In practice, such processes are relevant to model directional processes, for instance.}
The simulation study is performed on a $d \times d$ spatial unit grid (i.e., $D_{\boldsymbol{s}} = \{\boldsymbol{s} = (s_1, s_2)' \in \mathds{Z}^2 : 1 \leq s_1, s_2 \leq d \}$), resulting in $n = d^2$ observations, with $m = 10000$ replications. {The size of this spatial field has been successively increased with $d \in \{5, 10, 15\}$. Moreover, we have considered different settings depending on the data-generating parameters $\boldsymbol{\theta_0} = (\rho_0, \lambda_0, \alpha_0)'$. Whereas the unconditional variance parameter $\alpha_0$ equals 1 for all settings (i.e., variance level which is independent of the spatial location), $\rho_0$ and $\lambda_0$ varied across the settings. To be precise, $\rho_0 \in \{0.2, 0.5, 0.8\}$ and $\lambda_0 \in \{0.2, 0.5, 0.8\}$ to have settings with a weak, moderate and large dependence in the conditional spatial heteroscedasticity.}
{For all three parameters $\rho$, $\lambda$, and $\alpha$, the average RMSE is shown Table (ref). In general, our theoretical results are confirmed. That is, the parameters can correctly be estimated and the RMSE is decreasing with an increasing number of observations.} {Moreover, the non-linear least squares approach can efficiently be implemented and runs very fast on standard computers, even if the number of observations is large. On average, the computing time for estimating the parameters ranges from 3.4 to 3.8 seconds for 25 observations, 3.2 - 4.8 seconds for 100 observations, and 5.5 - 8.5 seconds for 225 observations on a standard notebook. This method is implemented in the R package spGARCH (version $> 2.0$).}
In markets that are constrained in space, one can typically expect to find locally varying risks. Typical examples of such markets are real estate and labour. For the former, the property prices are highly dependent upon the location of the real estate and prices in the surrounding areas. Similarly, for the latter, this market is also often constrained in space due to the limited mobility of labourers.
On the one hand, we observe conditional mean levels that vary in space, so-called spatial clusters. That is, both clustered areas of higher prices and lower prices can be observed. On the other hand, we may also expect to find locally varying risks relating to price, which can be considered as local volatility clusters. The proposed spatial GARCH-type models are capable of capturing such spatial dependencies in the conditional variance. This motivates why we consider condominium prices at a fine spatial scale. In particular, we will analyze the relative changes in Box-Cox transformed prices from 1995--2014 across all Berlin ZIP-Code regions (i.e., $n = 190$). The data are depicted in Figure (ref). The sample mean for these price changes is 0.8103 with a median of 0.6965. In total, the price changes range from -2.5650 to 7.3131. We can observe a spatial cluster of positive values in the north-western ZIP-Code regions.
{However, this fine spatial scale of all ZIP-Code regions causes another problem, namely exogenous regressors are often not available or cannot be assigned in the given small-scale resolution. For instance, the average household income could play an important role in the increase of the condominium prices, but the place of living (in terms of ZIP-Code areas) does not usually coincide with the place of work. There is no reliable way to associate quantities like personal or household income with ZIP-Code areas. Moreover, local infrastructure like schools, leisure facilities, parks is not limited to the residents of the respective ZIP-code areas. Thus, modelling an empirical process on such a small spatial scale is typically prone to heteroscedasticity induced by latent variables.}
{To illustrate these effects, we applied the developed spatial GARCH model to the residuals of a spatial autoregressive model, briefly SAR (see, e.g. Halleck15,Lee04), with and without exogenous regressors. More precisely, we select the regressors from a set of potential covariates available for different spatial scales, including the number of crimes (in so-called life-world oriented spaces, LOR, 161 units), number of schools, kinder gardens (ZIP-code level), percentage of migrants (LOR level), number of inhabitants (LOR level), size of areas used for infrastructure, living, water, vegetation (district level, 12 units), and average net income per household (district level). The included regressors were chosen such that the Akaike information criterion is minimal. The results of these two models are shown in Table (ref) along with the spGARCH coefficients of the error process. Regarding the residuals of these two mean models with and without regressors (before fitting a spGARCH model to the residuals), we observe that both of them are not autocorrelated in space (Moran's $I = -0.0375$ with a p-value of $0.7811$ for the intercept-only model, and $I = -0.0124$ with $p = 0.5676$ for the model including covariate effects). That is, the spatial correlation of original data could fully be modelled ($I = 0.5104$). However, looking at the absolute values of the residuals, we observe that there is significant autocorrelation for both models ($I = 0.1027$ ($p = 0.0042$) and $I = 0.0808$ ($p = 0.0181$) for the intercept-only model and regressive model, respectively).}
{In the final regressive model, only four regressors were selected, namely a proxy for the available free space (i.e., the proportion of settlement area to total area), the net household income in each district, and linear trends in the east-west direction, and north-south direction have been included (i.e., the coordinates of the centroids of each ZIP-code unit). While a significant negative trend can be observed from west to east, the increase in the north-south direction seems to be of minor importance. Furthermore, the average household income clearly influences the price development of condominiums in Berlin. Condominiums in high-class districts in terms of average income have increased relatively less than condominiums in lower-income areas.}
{It is worth noting that the dependence in the conditional heteroscedasticity could partly be covered by these covariates and the spatial autocorrelation in the absolute residuals was reduced. Nevertheless, latent effects are present which were not modelled by including these regressors. Thus, an spGARCH model has been fitted in a second step to the residuals of both models. The obtained spGARCH parameters can be interpreted as local model uncertainties and simultaneously cover latent variables, which could not be included due to the fine spatial scale of ZIP-Code levels.}
{In both cases, we observe significant positive dependence in the conditional second moments. More precisely, the GARCH effects amount to $\hat{\rho} = 0.2136$ and $\hat{\rho} = 0.2022$ with $\hat{\lambda}$ roughly equal to 0.70 for the intercept-only and regression model, respectively. These parameters can similarly be interpreted as in the time-series case, although $\boldsymbol{h}$ does not necessarily coincide with the conditional second moments (see Otto19_statpapers). Analyzing the residuals of this combined model shows that the remaining dependence in the heteroscedasticity could be explained. The squared residuals are no longer significantly autocorrelated in space.}
{The resulting conditional variance is visualized for the model with regressors in Figure (ref). The highest values of $h(\boldsymbol{s}_i)$ can be observed for the outer regions in north-west and another cluster is located in the southern city centre. This indicates that the highest uncertainty in the price changes is observed for these regions, while there is a band around the city centre where the price changes could more accurately be predicted by the regression model (i.e., the estimated values of $h(\boldsymbol{s}_i)$ are lower). These results seem to be very reasonable because the real-estate market was changing the most in these areas. First, due to increasing prices in the centre, new land for building has been created outside the city -- mostly along the regional transport tracks which are mainly going in an east-west direction. Second, the major airport of West Berlin, namely the airport Berlin-Tempelhof, was closed in 2008 changing the atmosphere in this region, some ZIP-code areas changed from regions in the flight paths to calm regions very close to the city centre, while others were not affected. This explains the second cluster in the centre.}
Recently, a few papers have introduced spatial ARCH and GARCH-type models that allow the modelling of an instantaneous spatial autoregressive dependence of heteroscedasticity. In this paper, we propose a {generalized spatial ARCH model} that {additionally} covers all previous approaches. Due to the flexible definition of the model as a set of functions, we can derive a common estimation strategy for all these spatial GARCH-type models. {It} is based on non-linear least squares.
{In the second part of the paper, we confirmed our theoretical findings on the consistency of the estimators by means of Monte Carlo simulation studies. The estimation method is computationally implemented in the R package spGARCH.} Eventually, the use of the model was demonstrated through an empirical example. More precisely, this paper has shown how the model uncertainties of local price changes in the real estate market in Berlin can be described using an spGARCH model {as residuals' process}. Though all proposed models are uncorrelated and have a zero mean, potential interactions between the error process and the mean equation should be analyzed in greater detail in future research.
In addition, we want to stress that the dependence structure does not necessarily have to be interpreted in a spatial sense. Thus, we briefly discuss a further example below, on which the “spatial” proximity could also be defined as the edges of networks. In such cases, $\mathbf{W}_1$ and $\mathbf{W}_2$ would be interpreted as adjacency matrices. For instance, one might consider the financial returns of several stocks as a network, where the only assets that are connected are those that are correlated above a certain threshold. This could create a financial network, as shown in Figure (ref). Thus, spGARCH models can be used to analyze various forms of information, whether that might be volatility, risk, or spill-overs from one stock to another, if these assets are close to one another within a certain network. In future research, attempts for modelling volatility clusters within networks, using spatial GARCH models, should be analyzed in greater detail.
{Up to now}, we have assumed that suitable functions of the {spGARCH} model framework are known. Hence, it is possible to maximize certain goodness-of-fit criteria in order to obtain the best-fitting model. However, these functions can also be estimated using a non-parametric approach; for instance by penalized or classical B-splines. {Besides, further choices of $f$ have not been discussed in this paper yet, including choices of $f$ to obtain E-spGARCH or logarithmic spGARCH models. Also, multivariate models remain open for future research. This} will be the subject of {some} forthcoming papers.