spGARCH: An R-Package for Spatial and Spatiotemporal ARCH models

TL;DR

This paper introduces the spGARCH R-package for modeling spatial and spatiotemporal ARCH processes, including new models and estimation methods, demonstrated through simulations and a real-data application.

Contribution

It presents new spatial ARCH models, including an exponential variant and a complex-valued form, with maximum-likelihood estimators and practical implementation in R.

Findings

Successful simulation of spatial ARCH models

Estimation of models using the spGARCH package

Application to real spatial data example

Abstract

In this paper, a general overview on spatial and spatiotemporal ARCH models is provided. In particular, we distinguish between three different spatial ARCH-type models. In addition to the original definition of Otto et al. (2016), we introduce an exponential spatial ARCH model in this paper. For this new model, maximum-likelihood estimators for the parameters are proposed. In addition, we consider a new complex-valued definition of the spatial ARCH process. From a practical point of view, the use of the R-package spGARCH is demonstrated. To be precise, we show how the proposed spatial ARCH models can be simulated and summarize the variety of spatial models, which can be estimated by the estimation functions provided in the package. Eventually, we apply all procedures to a real-data example.

Tables

Table 1. Table 1: Overview of all types of spARCH models implemented in the spGARCH package.

Process type	Definition of $𝒉$	Comments
spARCH	${𝒉}_O=\alpha \mathrm{𝟏}+\rho 𝐖{\left(𝐈-𝐀\right)}^{-1}\left(\alpha {𝜺}^{\left(2\right)}\right)$	$𝜺$ is simulated from multivariate normal distribution (MN) truncated on the interval $[-1/\sqrt[4]{{\Vert {\rho }^2{𝐖}^2\Vert }_1},1/\sqrt[4]{{\Vert {\rho }^2{𝐖}^2\Vert }_1}]$
spARCH (oriented)	${𝒉}_O=\alpha \mathrm{𝟏}+\rho 𝐖{\left(𝐈-𝐀\right)}^{-1}\left(\alpha {𝜺}^{\left(2\right)}\right)$	$𝜺\sim \text{MN}(\mathrm{𝟎},𝐈)$, $𝐖$ must be a strictly triangular weighting matrix
spatial E-ARCH	$\ln {𝒉}_E=𝐒\left(\alpha \mathrm{𝟏}+\rho b𝐖\ln \left|𝒀\right|\right)$	$𝜺\sim \text{MN}(\mathrm{𝟎},𝐈)$, but moments of $𝒀$ differ from the moments of classical spARCH process (cf. Otto et al., 2016)
spARCH (complex)	${𝒉}_C=\alpha \mathrm{𝟏}+\rho 𝐖{\left(𝐈-𝐀\right)}^{-1}\left(\alpha {𝜺}^{\left(2\right)}\right)$	$𝜺\sim \text{MN}(\mathrm{𝟎},𝐈)$, but complex-valued $𝒀$

Table 2. Table 2: Summary of the main functions of the spGARCH package.

Function	Description
Main functions
sim.spARCH()	Simulation of spARCH and E-spARCH processes
qml.spARCH()	Quasi-maximum-likelihood estimation for spARCH models
qml.SARspARCH()	Quasi-maximum-likelihood estimation for SAR models with spARCH residuals
Generic methods
summary()	Summary of an object of ‘spARCH’ class generated by qml.spARCH() or qml.SARspARCH()
print()	Printing method for ‘spARCH’ class or summary.spARCH class
fitted()	Extracts the fitted values of an object of ‘spARCH’ class
residuals()	Extracts the residuals of an object of ‘spARCH’ class
logLik()	Extracts the log-likelihood of an object of ‘spARCH’ class
extractAIC()	Extracts the AIC of an object of ‘spARCH’ class
plot()	Provides several descriptive plots of the residuals of an object of ‘spARCH’ class

Table 3. Table 3: Overview of spatial models, which can be fitted by qml.spARCH() and qml.SARspARCH() .

Function	formula	type	Resulting model
qml.spARCH()	y $\sim$ 0	"gaussian"	spARCH model (see (1) and (2))
qml.spARCH()	y $\sim$ 1	"gaussian"	spARCH model with an additional intercept for the mean equation
qml.spARCH()	y $\sim$ a + b	"gaussian"	Linear Regression with regressors a and b and spARCH residuals
qml.spARCH()	y $\sim$ a + b:c	"gaussian"	Linear Regression with more complex expressions and spARCH residuals
qml.spARCH()	y $\sim$ 0	"exp"	E-spARCH model (see (1) and (5))
qml.spARCH()	y $\sim$ 1	"exp"	E-spARCH model with an additional intercept for the mean equation
qml.spARCH()	y $\sim$ a + b	"exp"	Linear Regression with regressors a and b and E-spARCH residuals
qml.spARCH()	y $\sim$ a + b:c	"exp"	Linear Regression with more complex expressions and E-spARCH residuals
qml.SARspARCH()	y $\sim$ 0	"gaussian"	SAR model without an intercept, but with spARCH residuals (see (8) and (9))
qml.SARspARCH()	y $\sim$ 1	"gaussian"	SAR model with an intercept and spARCH residuals
qml.SARspARCH()	y $\sim$ a + b	"gaussian"	SAR model with an intercept and the regressors a and b and spARCH residuals
qml.SARspARCH()	y $\sim$ a + b:c	"gaussian"	SAR model with more complex expressions and spARCH residuals
qml.SARspARCH()	y $\sim$ 0	"exp"	SAR model without an intercept, but with E-spARCH residuals (see (8) and (9))
qml.SARspARCH()	y $\sim$ 1	"exp"	SAR model with an intercept and E-spARCH residuals
qml.SARspARCH()	y $\sim$ a + b	"exp"	SAR model with an intercept and the regressors a and b plus E-spARCH residuals
qml.SARspARCH()	y $\sim$ a + b:c	"exp"	SAR model with more complex expressions and E-spARCH residuals

Table 4. Table 4: Overview of all included regressors and factor loading for the 10 common factors. The regressors were divided into 5 subgroups to allow for distinctions between the factors.

	F. 1	F. 2	F. 3	F. 4	F. 5	F. 6	F. 7	F. 8	F. 9	F. 10
$P{M}_{2.5}$ concentration	0.69	0.72
$S{O}_2$ concentration	0.33	-0.03
$N{O}_2$ concentration	0.13	-0.12
$CO$ concentration	0.31	0.05
$P{M}_{10}$ concentration	0.07	0.44
$O_3$ concentration	1.00	-0.02
Solar radiation			0.60	0.44
Precipitation			-0.08	-0.26
Outdoor temperature			1.00	-0.05
Temperature differences			0.32	0.94
Ambient maximal temperature			0.08	-0.39
$C{H}_{2}O$	-0.23	0.32
Percentage of current smokers					0.47	-0.85
Percentage of former smokers					0.92	0.37
Smoke some days					-0.07	-0.62
Never smoked					-0.96	0.25
Aerobic activity							-0.05	0.58
Exercises							0.41	0.33
Physical activity index							-0.09	0.99
Alcohol consumption					0.04	0.62
Binge drinking					0.07	0.44
Heavy drinking					0.43	0.02
High cholesterol									0.00	1.00
Cholesterol checked									0.55	0.00
Overweight (BMI 25.0-29.9)									0.99	0.09
Obese (BMI 30.0 - 99.8)									-0.75	0.01
Blood stool test									0.56	-0.23
Sigmoidoscopy									0.14	-0.16
High blood pressure									0.03	0.79
Flu shot							0.81	-0.13
Pneumonia vaccination							0.51	-0.26
Health care coverage							0.58	0.18
Seatbelt use							-0.58	0.10

Equations

$$\mathbf{Y}=\mathrm{diag}(\mathbf{h}{)}^{1/2}\mathbf{\epsilon },$$

$$Var(Y({\mathbf{s}}_i)|Y({\mathbf{s}}_1),\dots ,Y({\mathbf{s}}_{i-1})),$$

$${\mathbf{h}}_O=\alpha 1+\rho \mathbf{W}\mathrm{diag}(\mathbf{Y})\mathbf{Y},$$

$${\mathbf{Y}}^{(2)}=\alpha {\left(\mathbf{I}-\mathbf{A}\right)}^{-1}{\mathbf{\epsilon }}^{(2)},$$

$$\mathbf{A}=\rho \mathrm{diag}\left(\epsilon ({\mathbf{s}}_1{)}^2,\dots ,\epsilon ({\mathbf{s}}_n{)}^2\right)\mathbf{W}.$$

$${\mathbf{h}}_O=\alpha 1+\rho \alpha \mathbf{W}{\left(\mathbf{I}-\mathbf{A}\right)}^{-1}{\mathbf{\epsilon }}^{(2)}.$$

$$(\mathbf{I}-{\mathbf{A}}^2{)}^{-1}$$

$$a=\left\{\begin{array}[]{cc}\infty&\exists k>0:\rho\mathbf{W}^{k}=\mathbf{0}\\ {1}/{\sqrt[4]{\rho^{2}||\mathbf{W}^{2}||_{1}}}&\rho^{2}||\mathbf{W}^{2}||_{1}>0\end{array}\right.\,,$$

$$\ln {\mathbf{h}}_E=\alpha 1+\rho \mathbf{W}{g}_b(\mathbf{\epsilon }),$$

$$g_b(\mathbf{\epsilon })=(\ln |\epsilon ({\mathbf{s}}_1){|}^b,\dots ,\ln |\epsilon ({\mathbf{s}}_n){|}^b{)}^{\prime }$$

$$\ln h_E({\mathbf{s}}_i)=\alpha +\sum _{v=1}^n\rho b{w}_{iv}\ln |\epsilon ({\mathbf{s}}_v)|\text{for}i=1,\dots ,n.$$

$$\ln {\mathbf{h}}_E=\mathbf{S}\left(\alpha 1+\rho b\mathbf{W}\ln |\mathbf{Y}|\right)$$

$$\mathbf{S}=(s_{ij}{)}_{i,j=1,\dots ,n}={\left(\mathbf{I}+\frac{1}{2}\rho b\mathbf{W}\right)}^{-1}.$$

$$h_C({\mathbf{s}}_i)=\alpha +\sum _{v=1}^n\rho w_{iv}Y({\mathbf{s}}_v{)}^2.$$

$$\mathbf{W}=\left(\begin{array}[]{cccc}\mathbf{W}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\ \mathbf{I}&\mathbf{W}_{2}&\cdots&\mathbf{0}\\ \vdots&\vdots&\ddots&\vdots\\ \mathbf{I}&\mathbf{I}&\cdots&\mathbf{W}_{T}\\ \end{array}\right)\,.$$

$$\rho\left(\begin{array}[]{cccc}\mathbf{W}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\ \mathbf{0}&\mathbf{W}_{2}&\cdots&\mathbf{0}\\ \vdots&\vdots&\ddots&\vdots\\ \mathbf{0}&\mathbf{0}&\cdots&\mathbf{W}_{T}\\ \end{array}\right)+\phi_{1}\left(\begin{array}[]{cccc}\mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\ \mathbf{I}&\mathbf{0}&\cdots&\mathbf{0}\\ \vdots&\vdots&\ddots&\vdots\\ \mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\ \end{array}\right)+\ldots\,$$

$$\mathbf{Y}=\lambda \mathbf{B}\mathbf{Y}+\mathbf{X}\mathbf{\beta }+\mathbf{u}$$

$$\mathbf{u}=\mathrm{diag}(\mathbf{h}{)}^{1/2}\mathbf{\epsilon }.$$

$$f_{\mathbf{Y}}(\mathbf{y})$$

$$(\hat{\alpha },\hat{\rho })=\underset{\alpha >0,\rho \ge 0}{\mathrm{argmax}}\ln |\det \left({\left(\frac{\partial y_j/\sqrt{h_j}}{\partial y_i}\right)}_{i,j=1,\dots ,n}\right)|+\sum _{i=1}^n\ln f_{\epsilon }(y_i).$$

$$\frac{\partial y_{j}/\sqrt{h_{j}}}{\partial y_{i}}=\left\{\begin{array}[]{ccc}1\,/\,\sqrt{h}_{j}&\mbox{for}&i=j\\ -\frac{y_{i}y_{j}}{h_{j}^{3/2}}\rho w_{ji}&\mbox{for}&i\neq j\end{array}\right.\,.$$

$$\frac{\partial y_{j}/\sqrt{h_{j}}}{\partial y_{i}}=\left\{\begin{array}[]{ccc}1\,/\,\sqrt{h_{j}}&\mbox{for}&i=j\\ -\frac{by_{j}}{2y_{i}h_{j}^{3/2}}\rho s_{ji}w_{ji}&\mbox{for}&i\neq j\end{array}\right.\,$$

$$h_j=\exp \left(\sum _{v=1}^n{s}_{jv}\left(\alpha +\rho w_{jv}\ln |y_v|\right)\right).$$

$$\ln |\det \left(\text{diag}\left(\frac{h_1}{y_1^2},\dots ,\frac{h_n}{y_n^2}\right)-\rho {\mathbf{W}}^{\prime }\right)|+\sum _{i=1}^n\ln \frac{y_i^2}{h_i^{3/2}}$$

$$\ln |\det \left(\mathrm{diag}\left(\frac{2h_1}{b},\dots ,\frac{2h_n}{b}\right)-\rho {\mathbf{S}}^{\prime }\circ {\mathbf{W}}^{\prime }\right)|+\sum _{i=1}^n\ln \frac{b}{2h_i^{3/2}}.$$

$$\ln {\mathbf{h}}_E=\mathbf{S}\left(\alpha 1+\rho b\mathbf{W}\ln |\mathbf{Y}|\right)$$

$$\mathbf{S}=(s_{ij}{)}_{i,j=1,\dots ,n}={\left(\mathbf{I}+\frac{1}{2}\rho b\mathbf{W}\right)}^{-1}.$$

$$\det \left(\mathbf{I}+\frac{1}{2}\rho b\mathbf{W}\right)\ge 1+\frac{1}{2}\rho b\det (\mathbf{W})>0.$$

$$\mathbf{S}=(s_{ij}{)}_{i,j=1,\dots ,n}={\left(\mathbf{I}+\frac{1}{2}\rho b\mathbf{W}\right)}^{-1}.$$

$$\ln {\mathbf{h}}_E=\mathbf{S}\left(\alpha 1+\rho b\mathbf{W}\ln |\mathbf{Y}|\right)$$

Full text

Keywords: Spatial ARCH model, exponential ARCH model, R, SARspARCH model, spatiotemporal statistics, variance clusters

1 Introduction

Whereas autoregressive conditional heteroscedasticity (ARCH) models are applied widely in time series analysis, especially in financial econometrics, spatial conditional heteroscedasticity has not been seen as critical issue in spatial econometrics up to now. Although it is well-known that classical least squares estimators are biased for spatially correlated data as well as for spatial data with an inhomogeneous variance across space, there are just a few papers proposing statistical models accounting for spatial conditional heteroscedasticity in terms of the ARCH and GARCH models of Engle, (1982) and Bollerslev, (1986). The first extensions to spatial models attempted were time series models incorporating spatial effects in temporal lags (see Borovkova and Lopuhaa, 2012 and Caporin and Paruolo, 2006, for instance). Instantaneous spatial autoregressive dependence in the conditional second moments, i.e., the conditional variance in each spatial location is influenced by the variance nearby, has been introduced by Otto et al., (2016); Otto et al., 2018a ; Otto et al., 2018b . Their models allow for these instantaneous effects but require certain regularity conditions. In this paper, we propose an alternative specification of spatial autoregressive conditional heteroscedasticity based on an exponential definition of the conditional variance. This new model can be seen as the spatial equivalent of the exponential GARCH model by Nelson, (1991). Other recent papers propose a mixture of these two approaches (see Sato and Matsuda, 2017; Sato and Matsuda, 2018a ). Moreover, all these models can be used in spatiotemporal settings (see Otto et al., 2016; Sato and Matsuda, 2018b ).

In addition to the novel spatial exponential ARCH model, this paper demonstrates the use of the R-package spGARCH. From this practical point of view, the simulation of several spatial ARCH-type models as well as the estimation of a variety of spatial models with conditional heteroscedasticity are shown. There are several packages implementing geostatistical models, kriging approaches, and other spatial models (cf. Cressie, 1993; Cressie and Wikle, 2011). One of the most powerful packages used to deal with models of spatial dependence is spdep, written by Bivand and Piras, (2015). It implements most spatial models in a user-friendly way, such as spatial autoregressive models, spatial lag models, and so forth (see, also, Elhorst, 2010 for an overview). These models are typically called spatial econometrics models, although they are not tied to applications in economics. In contrast, the package gstat provides functions for geostatistical models, variogram estimation, and various kriging approaches (see Pebesma, 2004 for details). For dealing with big geospatial data, the Stem package uses an expectation-maximization (EM) algorithm for fitting hierarchical spatiotemporal models (see Cameletti, 2015 for details). For a distributed computing environment, the MATLAB software D-STEM from Finazzi and Fasso, (2014) also provides powerful tools for dealing with heterogeneous spatial supports, large multivariate data sets, and heterogeneous spatial sampling networks. Additionally, these fitted models are suitable for spatial imputation. Contrary to these EM approaches, Bayesian methods for modeling spatial data are implemented in the R-INLA package (see Rue et al., 2009 for technical details of the integrated nested Laplace approximations and Martins et al., 2013 for recently implemented features). Along with this package, the R-INLA project provides several functions for diverse spatial models incorporating integrated nested Laplace approximations.

In contrast to the above mentioned software for spatial models, the prevalent R-package for time series GARCH-type models is rugarch from Ghalanos, (2018). Since spGARCH has been developed mainly to deal with spatial data, we aim to provide a package which is user-friendly for researchers and data scientists working in applied spatial science. Thus, the package is coordinated with the objects and ideas of R packages for spatial data rather than packages for dealing with time series.

We structure the paper as follows. In the next Section 2, we discuss all covered spatial and spatiotemporal ARCH-type models. In addition, we introduce a novel exponential spatial ARCH model, which has weaker regularity conditions than the other spatial ARCH models. In the subsequent section, parameter estimation based on the maximum-likelihood principle is discussed for both the previously proposed spatial ARCH models as well as the new exponential spatial ARCH model. After these theoretical sections, we demonstrate the use of the R-package spGARCH in Section 4. Further, we fit a spatial autoregressive model with exogenous regressors and spatial ARCH residuals for a real-world data set. In particular, we analyze prostate cancer incidence rates in southeastern U.S. states. Section 6 concludes the paper.

2 Spatial ARCH-type models

Let $\left\{Y(\boldsymbol{s})\in\mathds{R}:\boldsymbol{s}\in D\right\}$ be a univariate stochastic process having a spatial autoregressive structure in the conditional variance. The process is defined in a multidimensional space $D$, which is typically a subset of the $q$-dimensional real numbers $\mathds{R}^{q}$, as space is usually finite. For dealing with spatial lattice data, $D$ is subset of the $q$-dimensional integers $\mathds{Z}^{q}$. For both cases, it is important that the subset contains a $q$-dimensional rectangle of positive volume (cf. Cressie and Wikle, 2011). Moreover, this definition is suitable for modeling spatiotemporal data, as one might assume that $D$ is the product set $\mathds{R}^{k}\times\mathds{Z}^{l}$ with $k+l=d$.

To define spatial models, in particular areal spatial models such as the simultaneous autoregressive (SAR) models, it is convenient to consider a vector of observations $\boldsymbol{Y}=(Y(\boldsymbol{s}_{1}),\ldots,Y(\boldsymbol{s}_{n}))^{\prime}$ at all locations $\boldsymbol{s}_{1},\ldots,\boldsymbol{s}_{n}$. For spatial ARCH models, we specify this vector as

[TABLE]

an analogue to the well-known time series ARCH models (cf. Engle, 1982; Bollerslev, 1986). However, note that the vector $\boldsymbol{h}$ does not necessarily coincide with the conditional variance

[TABLE]

as the variance in any location $\boldsymbol{s}_{j}$ also depends on $Y(\boldsymbol{s}_{i})$ for $j\neq i$ (see Otto et al., 2016 for details). We now distinguish between several spatial ARCH-type models via the definition of $\boldsymbol{h}$.

2.1 Spatial ARCH model

First, we define this vector $\boldsymbol{h}$ in such a way as to be analogous to the definition in Otto et al., (2016); Otto et al., 2018a . For this model, the vector $\boldsymbol{h}_{O}$ is given by

[TABLE]

where $\mathrm{diag}(\boldsymbol{a})$ is a diagonal matrix with the entries of $\boldsymbol{a}$ on the diagonal. In order to be consistent with the implementation in the R-package spGARCH, we focus on the special case with two parameters $\alpha$ and $\rho$, whereas Otto et al., (2016) proposed a more general model with a vector $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{n})^{\prime}$ and the first-order spatial lag $\mathbf{W}\mathrm{diag}(\boldsymbol{Y})\boldsymbol{Y}$.

For this definition, there is a one-to-one relation between $\boldsymbol{Y}$ and $\boldsymbol{\varepsilon}$ via the squared observations $\boldsymbol{Y}^{(2)}=(Y(\boldsymbol{s}_{1})^{2},\ldots,Y(\boldsymbol{s}_{n})^{2})^{\prime}$ and squared errors $\boldsymbol{\varepsilon}^{(2)}=(\varepsilon(\boldsymbol{s}_{1})^{2},\ldots,\varepsilon(\boldsymbol{s}_{n})^{2})^{\prime}$ with

[TABLE]

where $\mathbf{W}$ is a predefined spatial weighting matrix and

[TABLE]

Thus,

[TABLE]

It is important to assume that the spatial weighting matrix is a non-stochastic, positive matrix with zeros on the main diagonal to ensure that a location is not influenced by itself (cf. Elhorst, 2010; Cressie and Wikle, 2011). The vector of random errors is denoted by $\boldsymbol{\varepsilon}$. Due to the complex dependence implied by the weighting matrix $\mathbf{W}$, $\boldsymbol{h}_{O}$ is not necessarily positive; thus, $\mathrm{diag}(\boldsymbol{h})^{1/2}$ does not necessarily have a solution in the real numbers such that the process in (1) is well-defined. This is only the case if the condition of the following lemma is fulfilled.

Lemma 1** (Otto et al., 2016).**

Suppose that $\alpha\geq 0$, $\rho\geq 0$ and that $\det(\mathbf{I}-\mathbf{A}^{2})\neq 0$. If all elements of the matrix

[TABLE]

are nonnegative, then all components of $\boldsymbol{Y}^{(2)}$ are nonnegative, i.e., $Y(\boldsymbol{s}_{i})^{2}\geq 0$ for $i=1,\ldots,n$. Moreover, $h_{O}(\boldsymbol{s}_{i})\geq 0$ for $i=1,\ldots,n$.

It is important to note that $\mathbf{A}$ depends on both the weighting matrix and the realizations of the errors. In order to ensure that this condition is fulfilled, Otto et al., (2016) propose to truncate the support of the error distribution on the interval $(-a,a)$ with

[TABLE]

where $||\cdot||_{1}$ denotes the matrix norm based on the $l_{1}$ vector norm.

There are two cases in which the support of the errors does not need to be constrained. If $\rho=0$, the process coincides with a spatial white noise process such that $a$ equals $\infty$. Moreover, all entries of $\boldsymbol{h}$ are non-negative if $\mathbf{W}$ is similar to a strictly triangular matrix. Then, $\mathbf{W}$ is nilpotent. This case covers the classical time-series ARCH($p$) models introduced by Engle, (1982) as well as the so-called oriented spARCH processes. For these processes, the spatial dependence has a certain direction, e.g., observations are only influenced by observations in a southward direction or by observations which are closer to an arbitrarily chosen center. The setting also covers recent time-series GARCH models incorporating spatial information (e.g., Borovkova and Lopuhaa, 2012; Caporin and Paruolo, 2006).

Of course, the truncated support of the errors has an impact on the extent of the spatial dependence on the conditional variances. Obviously, the support need not be constrained regarding $\rho=0$. However, this support decreases with increasing values of $\rho$. For instance, if $\rho=1$, then the parameter $a$ is equal to $0.968$ for Rook’s contiguity matrices on a two-dimensional lattice. As a measure of the spatial dependence of the variance, one might consider Moran’s $I$ for the squared observations (see Moran, 1950). Moreover, we observe that the growth rate of $I$ decreases with increasing spatial weights. This trend can be explained by the compact support of the errors. Since there cannot be large variations $\varepsilon(\boldsymbol{s}_{i})$ in absolute terms, there also cannot be large spatial clusters of high or low variance. To illustrate this behavior, Figure 1 depicts Moran’s $I$ for simulated observations $\boldsymbol{Y}$ and their squares for $\rho\in\{0,0.05,\ldots,2\}$. For the Monte Carlo simulation study, we simulate $n=400$ observation on a two-dimensional lattice $D=\{\boldsymbol{s}=(s_{1},s_{2})^{\prime}\in\mathds{Z}^{2}:0\leq s_{1},s_{2}\leq 20\}$. The weighting matrix is a common Rook’s contiguity matrix, and the simulation is done for $10^{5}$ replications. Although the exact distribution of Moran’s statistic is bounded, the standardized statistic is asymptotically normally distributed for the “majority of spatial structures” (Tiefelsdorf and Boots, 1995, see also Cliff and Ord, 1981). Thus, the asymptotic 95% confidence intervals are plotted in Figure 1, as well.

2.2 Exponential Spatial ARCH model

Next, we consider an exponential spatial ARCH process (E-spARCH). In this setting, we define the natural logarithm of $\boldsymbol{h}_{E}=(h_{E}(\boldsymbol{s}_{1}),\ldots,h_{E}(\boldsymbol{s}_{n}))^{\prime}$ as

[TABLE]

with a function $g_{b}:\mathds{R}^{n}\rightarrow\mathds{R}^{n}$. Like Nelson, (1991), we assume that

[TABLE]

for positive values of $b$. For this definition, there is a one-to-one relation between $\boldsymbol{Y}$ and $\boldsymbol{\varepsilon}$, as we show in the following theorem.

Theorem 1**.**

Suppose that $\alpha>0$, $\rho\geq 0$, and $w_{ij}\geq 0$ for all $i,j=1,\ldots,n$ and $g_{b}(\boldsymbol{\varepsilon})=(\ln|\varepsilon(\boldsymbol{s}_{1})|^{b},\ldots,\ln|\varepsilon(\boldsymbol{s}_{n})|^{b})^{\prime}$. Then there exists one and only one $Y(\boldsymbol{s}_{1}),\ldots,Y(\boldsymbol{s}_{n})$ that corresponds to each $\varepsilon(\boldsymbol{s}_{1}),\ldots,\varepsilon(\boldsymbol{s}_{n})$ for $b>0$.

At location $\boldsymbol{s}_{i}$, the value of $h_{E}(\boldsymbol{s}_{i})$ is then given by

[TABLE]

For this definition of $g_{b}$, one could rewrite $\ln\boldsymbol{h}$ as

[TABLE]

with

[TABLE]

In contrast to the spARCH process described in Section 2.1, Corollary 1 shows that the entries of $\boldsymbol{h}_{E}$ are positive for all $\rho\geq 0$ and $\alpha>0$. Hence, the process is well-defined and there are no further restrictions needed, as in the case for the spARCH model.

Corollary 1**.**

Assume that the assumptions of Theorem 1 are fulfilled, then $\boldsymbol{h}_{E}(\boldsymbol{s}_{i})\geq 0$ for all $i=1,\ldots,n$.

For all proofs, we refer to the Appendix.

2.3 Complex Spatial ARCH model

Now, we propose a complex-valued spARCH process. In order to obtain a solution of $\mathrm{diag}(\boldsymbol{h})^{1/2}$ in the $n$-dimensional space of real numbers for the model defined in (2), all elements of the matrix $(\mathbf{I}-\mathbf{A}^{2})^{-1}$ must be nonnegative (see Otto et al., 2016). For the complex spARCH process, we relax the assumption that there should be a solution to $\mathrm{diag}(\boldsymbol{h})^{1/2}$ in the real numbers and also consider complex solutions. Thus, the definition of $\boldsymbol{h}$ coincides with $\boldsymbol{h}_{O}$ of the original model, i.e.,

[TABLE]

2.4 Spatiotemporal ARCH model

Finally, we show that spatiotemporal processes are covered directly by these approaches. For spatiotemporal data, the vector $\boldsymbol{s}$ simply includes both the spatial location $\boldsymbol{s}_{s}$ and the point in time $t$, i.e., $\boldsymbol{s}=(\boldsymbol{s}_{s},t)^{\prime}$. In addition, it is important to assume that future observations do not influence past observations, i.e., the weights $w_{ij}$ must be zero if $t_{j}\geq t_{i}$. However, the dimension of the weighting matrix $\mathbf{W}$ might become very large for this representation. More precisely, the matrix has dimension $NT\times NT$, where $N$ is the total number of spatial locations and $T$ stands for the total number of time points. From a computational perspective, this is not necessarily a drawback since $\mathbf{W}$ is usually sparse and could also have a block diagonal structure. Moreover, it is often reasonable to assume that $h(\boldsymbol{s}_{i})$ is only influenced by the neighbors of $\boldsymbol{s}_{s,i}$ at the same point of time and by past observations at the same location. Then the weighting matrix would have the following structure

[TABLE]

Indeed, it is plausible to weight the spatial and temporal lags differently by replacing $\rho\mathbf{W}$ by a sum

[TABLE]

with positive weights $\phi_{k}$ for all temporal lags $1\leq k\leq p$.

2.5 Spatial ARCH Disturbances

Since all conditional and unconditional odd moments of spatial ARCH processes are equal to zero, these ARCH-type models can easily be added to any kind of (spatial) regression model without influencing the mean equation as well as the spatial dependence in the first conditional and unconditional moments. This makes the spatial ARCH models flexible tools for dealing with conditional spatial heteroscedasticity in the residuals of spatial models. For instance, one can consider spatial autoregressive models for $\boldsymbol{Y}$, i.e.,

[TABLE]

with $\boldsymbol{u}$ following either a spatial ARCH model with the original definition $\boldsymbol{h}_{O}$ or the exponential model with $\boldsymbol{h}_{E}$. Thus,

[TABLE]

Further, we call this model the SARspARCH model. For $\lambda=0$, the model collapses to a simple linear regression model; if, additionally, $\boldsymbol{\beta}=\boldsymbol{0}$, the model coincides with the previously discussed ARCH models. Thus, these coefficients can be used for testing against nested models.

In contrast to other models for heteroscedastic errors, such as the SARAR or SARMA models, which assume spatial autoregressive or spatial moving average error terms (cf. Kelejian and Prucha, 2010; Fingleton, 2008; Haining, 1978), the SARspARCH model does not affect the spatial autocorrelation of the process, just the spatial heteroscedasticity, because all conditional and unconditional odd moments are equal to zero. Thus, $\lambda\mathbf{B}$ can be interpreted directly as the spatial dependence of the process, while $\rho\mathbf{W}$ describes the spatial dependence in the second conditional moments. Moreover, these two parts can be interpreted separately, as we will demonstrate in the last section via an empirical example.

3 Parameter Estimation

The parameters of a spatial ARCH process can be estimated by the maximum-likelihood approach. To obtain the joint density for $\boldsymbol{Y}=k(\boldsymbol{\varepsilon})$, the Jacobian matrix of $k^{-1}$ at the observed values $\boldsymbol{y}$ must be computed (e.g., Bickel and Doksum, 2015). If $f_{\boldsymbol{\varepsilon}}$ is the distribution of the error process, then the joint density $f_{\boldsymbol{Y}}$ of $\boldsymbol{Y}$ is given by

[TABLE]

If the residuals are additionally independent and identically distributed, the parameter estimates can be obtained from the maximization of the log-likelihood as follows

[TABLE]

The Jacobian matrix, of course, depends on the definition of $\boldsymbol{h}$. For the spARCH process, this Jacobian matrix can be specified as

[TABLE]

In contrast, the Jacobian matrix for the E-spARCH process is slightly different, namely

[TABLE]

with

[TABLE]

From a computational perspective, the computation of the log determinant of this matrix is feasible for large data sets. To be precise, the log-determinant is equal to

[TABLE]

for the spARCH process. Similarly, it is given by

[TABLE]

for the E-spARCH process, where $\circ$ stands for the Hadamard product.

In the spGARCH package, we implemented the iterative maximization algorithm with inequality constraints proposed by Ye, (1988), which is implemented in the R-package Rsolnp (see Ghalanos and Theussl, 2012). It is important to note that the log determinant of the Jacobian also depends on the parameters in such a way that it needs to be computed in each iteration (see, also, Theorem 13.7.3 of Harville, (2008) for the computation of a determinant for the sum of a diagonal matrix and an arbitrary matrix), but $\mathbf{W}$, and therefore $\mathbf{S}\circ\mathbf{W}$, are usually sparse. Thus, the required time for the estimation of the parameters depends mainly on the dimension and sparsity of $\mathbf{W}$.

4 Overview of the R-Package spGARCH

The R-package spGARCH provides several basic functions for the analysis of spatial data showing spatial conditional heteroscedasticity. In particular, the process can be simulated for arbitrarily chosen weighting matrices according to the definitions in Section 2. Moreover, we implement a function for the computation of the maximum-likelihood estimators. To generate a user-friendly output, the object generated by the estimation function can easily be summarized by the generic summary() function. We also provide all common generic methods, such as plot(), print(), logLik(), and so forth. To maximize the computational efficiency, the actual version of the package contains compiled C++ code (using the packages Rcpp and RcppEigen, cf. Eddelbuettel and François, 2011; Bates and Eddelbuettel, 2013). A brief overview of the package and its main functions is given in Table 2. Further, we focus on the two main aspects of the package, i.e., the simulation (described in detail in Section 4.1) and estimation (Section 4.2) aspects of the spARCH, E-spARCH, and SARspARCH processes.

4.1 Simulation of ARCH-type stochastic processes

The simulations of all spatial ARCH-type models are implemented in one function, namely, the sim.spARCH() function. The different definitions of the model are specified via the argument type. The use of sim.spARCH() is very similar to how a basic random number generator is used, meaning that the first argument n is the number of generated values and all further arguments specify the parameters of the spARCH process. For instance, one might simulate an oriented spARCH process (meaning $\mathbf{W}$ is triangular) on a $d\times d$ spatial lattice with $\rho=0.7$ and $\alpha=1$ using the following lines.

⬇

1R> require("spdep")

2R> rho <- 0.7

3R> alpha <- 1

4R> d <- 50

5R> n <- d^2

6R> nblist <- cell2nb(d, d, type = "queen")

7R> W <- nb2mat(nblist)

8R> W[upper.tri(W)] <- 0

9R> Y <- sim.spARCH(n = n, rho = rho, alpha = alpha, W = W, type = "gaussian", control = list(seed = 5515))

To build the spatial weighting matrix, we used cell2nb() from the spdep package, returning an nb object of a $d\times d$ lattice (see Cressie, 1993; Bivand and Piras, 2015). Further, we converted the nb object into a contiguity matrix, as sim.spARCH() requires either a matrix (class matrix) or a sparse matrix (class dgCMatrix) as an argument. Usually, spatial weighting matrices are sparse by construction. Thus, $\mathbf{W}$ is always converted internally to a dgCMatrix matrix or rather to a SparseMatrix object defined in the eigen library in C++. Via the control parameter, a random seed might be passed to the simulation function. If not, a random seed is assigned randomly from a uniform distribution and printed in console in order that one might reproduce the result even without having a random seed specified in advance. We prefer to print a single number in the console rather than returning to the random number generator (RNG) state as an attribute of the returned vector. Thus, a random seed might either be passed as an optional argument to sim.spARCH() or set before calling sim.spARCH() by set.seed().

There are several types of spatial ARCH processes which can be simulated by sim.spARCH(). They are all specified by the argument type. If

• •

type = "gaussian", then the original spARCH process according to the definition in Otto et al., (2016) is simulated.

• –

If there exists a permutation such that $\mathbf{W}$ is a strictly triangular matrix, then the function simulates automatically an oriented spARCH process with independent and identically gaussian distributed errors.

• –

If there is no such permutation, then the errors are simulated from a truncated normal distribution with $a=1/\sqrt[4]{\rho^{2}||\mathbf{W}||_{1}}$.

• •

type = "exp", an E-spARCH process is simulated with an user-specified value of $b$ (default 2) and standard normal random errors.

• •

type = "complex", complex solutions of $\mathrm{diag}(\boldsymbol{h})^{1/2}$ are considered in order to simulate the spARCH process.

Figure 2 illustrates the behavior of different types of spatial ARCH processes. All of them are simulated with the same parameters and random seeds in such a manner that the vector $\boldsymbol{\varepsilon}$ is identical for all types of processes, except for the spARCH process with the truncated normal errors. In the first row, the spatial weighting is achieved via a strictly triangular Queen’s contiguity matrix, which means that the spatial dependence has its origin in the upper left corner. To the contrary, $\mathbf{W}$ presents a classical Queen’s contiguity matrix in the second row. We additionally plot a spatial white noise process for comparison, as we used a rather unconventional two-color scheme. Using this kind of color scheme, one might distinguish between positive and negative observations, such that it is easier to see the spatial volatility clusters. Areas of smaller volatility are characterized by rather evenly gray pixels, whereas clusters of high volatility have rather intense colors. Moreover, the colors fluctuate irregularly between blue and red.

4.2 Maximum-likelihood estimation

Other important functions of the package are the qml.spARCH() and qml.SARspARCH() functions, which implement a quasi-maximum-likelihood estimation algorithm (QML). As for the sim.spARCH() function, many spARCH models are covered in the qml.spARCH() and qml.SARspARCH() function. Thus, the user needs to specify which particular spARCH model is to be fitted via the argument type. Moreover, the model for the mean equation is a user-specified formula, making the use of the estimation functions similar to the use of the common lm() or glm() functions.

In general, the estimators exhibited good performances for a variety of error distributions in simulation studies, although the likelihood function was derived under the normality assumption. This is not surprising, as the maximum-likelihood estimators have good properties under mild assumptions for the error processes of a variety of similar spatial econometrics models (cf. Lee, 2004; Lee and Yu, 2012; Lee and Yu, 2010b ; Lee and Yu, 2010a ). Thus, we refer to the approach as the QML approach, and the name of the estimation functions start with qml instead of ml. In the following paragraphs, we start the simulation of one specific sample, which is then used further to illustrate the log-likelihood functions as well as to demonstrate parameter estimation.

Compared to the E-spARCH processes, the likelihood functions of spARCH models are rather flat around the global maximum. This behavior is illustrated for simulated processes in Figure 3. The observations for the E-spARCH process have been simulated as follows.

⬇

1R> nblist <- cell2nb(20, 20, type = "queen")

2R> W <- nb2mat(nblist)

3R> y <- sim.spARCH(n = 20^2, rho = 0.5, alpha = 1, W = W, type = "exp", control = list(seed = 5515))

To simulate an oriented process, the entries of $\mathbf{W}$ above the diagonal must be set to zero and the argument type must be changed to "gaussian", i.e.,

⬇

1R> W[upper.tri(W)] <- 0

2R> y2 <- sim.spARCH(n = 20^2, rho = 0.5, alpha = 1, W = W, type = "gaussian", control = list(seed = 5515))

To estimate the parameters of an intercept-free E-spARCH model without any regressors, the formula passed to the function qml.spARCH() should be specified as y $\sim$ 0. In addition, a data.frame can be passed via the data argument to the qml functions. Although the likelihood function of a spARCH process is flat, good estimates can be obtained through iterative maximization. Otto et al., (2016) analyze the performance of the estimators in detail. The algorithm implemented in the packages is based on the Rsolnp package, allowing for both equality and inequality parameter constraints (cf. Ghalanos and Theussl, 2012).

The results of the estimation procedure are returned via an object of the class ‘spARCH’, for which we provide additionally several generic functions. First, there is a summary() function for the ‘spARCH’ object. The summary shows all important estimation results, i.e., the parameter estimates, standard errors, test statistics, and asymptotic p-values, including significance stars. The estimation of the above simulated E-spARCH process would return the following results.

⬇

1R> spARCH_object <- qml.spARCH(y ~ 0, W = W, type = "exp")

2R> summary(spARCH_object)

3 Call:

4qml.spARCH(formula = y ~ 0, W = W, type = "exp")

5

6 Residuals:

7 Min. 1st Qu. Median Mean 3rd Qu. Max.

8-2.6867629 -0.6197315 -0.0053580 -0.0002615 0.5708346 2.8576621

9

10 Coefficients:

11 Estimate Std. Error t value Pr(>|t|)

12alpha 0.919324 0.128544 7.1518 8.564e-13 ***

13rho 0.402998 0.056519 7.1304 1.001e-12 ***

14---

15Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

16

17 AIC: 543.01, BIC: 539.01 (Log-Likelihood: -269.51)

18

19 Moran’s I (residuals): -0.028568, p-value: 0.31795

20

21 Moran’s I (squared residuals): 0.035239, p-value: 0.14479

The standard errors are estimated as Cramer-Rao bounds from the Hessian matrix of the log-likelihood function. For triangular weighting matrices, the estimators are asymptotically normally distributed (Otto et al., 2016). In addition to the Akaike and Bayesian Schwarz information criteria, the results of Moran’s test on the residuals and squared residuals are reported for the spatial autocorrelation of the residuals. However, it is possible to use functions like AIC() or BIC(), since there is a logLik() method for the objects from class ‘spARCH’. Additionally, the fitted values and residuals can be extracted by fitted() and residuals(), respectively.

To analyze the residuals, we provide additionally several descriptive plots via the generic plot() function. The first two plots are produced by moran.plot() imported from the package spdep. They inspect the spatial autocorrelation of the residuals and the squared residuals. In addition, the error distribution is depicted in the third graphic by a normal Q-Q-plot. The output obtained for the above numerical example is given below and in Figure 4.

⬇

1R> AIC(spARCH_object)

2[1] 543.0126

3R> BIC(spARCH_object)

4[1] 550.9956

5R> par(mfcol = c(1,3))

6R> plot(spARCH_object)

7Reproduce the results as follows:

8 eps <- residuals(x)

9 W <- as.matrix(x$W)

10 moran.plot(eps, mat2listw(W), zero.policy = TRUE,

11 xlab = "Residuals", ylab = "Spatially Lagged Residuals")

12Reproduce the results as follows:

13 eps <- residuals(x)

14 W <- as.matrix(x$W)

15 moran.plot(eps, mat2listw(W), zero.policy = TRUE,

16 xlab = "Residuals", ylab = "Spatially Lagged Residuals")

17Reproduce the results as follows:

18 eps <- residuals(x)

19 std_eps <- (eps - mean(eps))/sd(eps)

20 qqnorm(eps, ylab = "Standardized Residuals")

21 qqline(eps)

The mean equation can be specified as formula for all models, i.e., the spARCH, E-spARCH, and SARspARCH models. Thus, there is a huge variety of possible spatial ARCH models as well as regression models with spARCH residuals which can be fitted by the estimation functions. In addition to linear models of the form y $\sim$ a + b, more sophisticated models can also be fitted, e.g., models with interactions y $\sim$ a + b:c, factor models y $\sim$ factor, polynomial models y $\sim$ poly(a, 3), seasonally or regularly varying models of the form y $\sim$ sin(t) + cos(t) or y $\sim$ sin(long) + cos(long) + sin(lat) + cos(lat), and so forth. We also implement an extractAIC() method for ‘spARCH’ objects, such that one might also use step() for stepwise model selection. Table 3 provides an overview of possible combinations of the arguments formula and type and shows the resulting models, which can be fitted by the functions qml.spARCH() and qml.SARspARCH(), respectively.

5 Real-data example: prostate cancer incidence rates

Below, the focus is on the incidence rates (2008–2012) for prostate cancer provided by the Centers for Disease Control and Prevention (U.S. Department of Health and Human Services, Centers for Disease Control and Prevention and National Cancer Institute, 2015). In particular, we analyze the incidence rates in all counties of several southeastern U.S. states, namely Arkansas, Louisiana, Mississippi, Tennessee, North and South Carolina, Georgia, Alabama, and Florida. This area also covers the counties along the Mississippi River collectively known as “cancer alley” (see Nitzkin, 1992; Brent, 2010; Berry, 2003). All rates are age-adjusted to the 2000 U.S. standard population (cf. U.S. Department of Health and Human Services, Centers for Disease Control and Prevention and National Cancer Institute, 2015).

As explanatory variables, we included a large set of environmental, climate, behavioral, and health covariates, which might have an influence on incidence rates for prostate cancer. For instance, we consider air pollution, such as $PM_{2.5}$, $PM_{10}$, $SO_{2}$, $NO_{2}$, $CO$, $O_{3}$, and $CH2O$, as potential environmental hazard factors. Moreover, we account for smoking, drinking, sport activities, and further healthcare-related variables as potential influences on the cancer incidence rates. In total, we account for 34 explanatory variables, which were obtained by inverse-distance-kriging from spatial points processes. Most of the variables are correlated, so we performed a factor analysis on 5 subgroups to identify 10 common factors. The factor loadings are summarized in Table 4. Eventually, the final explanatory factors were chosen by minimizing the Bayesian information criterion using the generic function step() as follows.

⬇

1R> out <- step(qml.SARspARCH(formula, B = B, W = W, type = "gaussian"), k = log(length(Y)))

The formula object simply defines a linear model between the logarithmic incidence rates and all factors. Further, matrix $\mathbf{B}$ describes the predefined spatial dependence structure in the mean equation. For this analysis, $\mathbf{B}$ has been chosen as a row-standardized contiguity matrix of the direct neighbors. For the spatial dependence in the spatial ARCH term of the residuals, we also included all neighbors up to order 4. Hence, $\mathbf{W}$ is the row-standardized matrix of the sum of the first-, second-, third-, and fourth-lag neighbors.

By minimizing the BIC criterion, the $2^{\text{nd}}$ and $10^{\text{th}}$ factor has been selected. Whereas the $2^{\text{nd}}$ factor has positive loadings mainly for fine particulate matters, $PM_{2.5}$ and $PM_{10}$, the $10^{\text{th}}$ describes the tendency for high blood pressure and cholesterol in the county’s population. However, note that this analysis is based on aggregated data rather than individual patients; hence, the selected factors cannot be interpreted as carcinogenic factors.

Using the generic summary() for the ‘spARCH’ class, the estimated model can be summarized as follows.

⬇

1 Call:

2qml.SARspARCH(formula = formula, B = B, W = W, type = "gaussian",

3 data = NULL)

4

5 Residuals:

6 Min. 1st Qu. Median Mean 3rd Qu. Max.

7-0.7492270 -0.1079639 -0.0001509 -0.0005261 0.1121190 0.6404564

8

9 Coefficients:

10 Estimate Std. Error t value Pr(>|t|)

11alpha (spARCH) 0.0203839 0.0042674 4.7766 1.783e-06 ***

12rho (spARCH) 0.3782104 0.1309656 2.8879 0.003879 **

13lambda (SAR) 0.6768133 0.0356765 18.9708 < 2.2e-16 ***

14(Intercept) 1.5388985 0.1702222 9.0405 < 2.2e-16 ***

15X_factor_scores[, 2] 0.0192857 0.0069917 2.7584 0.005809 **

16X_factor_scores[, 10] -0.0205693 0.0064450 -3.1915 0.001415 **

17---

18Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

19

20 AIC: -1734.2, BIC: -1746.2 (Log-Likelihood: 873.11)

21

22 Moran’s I (residuals): -0.022899, p-value: 0.32023

23

24 Moran’s I (squared residuals): 0.021409, p-value: 0.00050052

First, we see that the model has a significant spatial autocorrelation in the mean equation since $\hat{\lambda}$ (lambda (SAR)) differs significantly from zero. This implies that there are clusters of higher prostate cancer incidence rates and, vice versa, lower incidence rates. Second, the error process shows conditional, autoregressive heteroscedasticity in space, which is captured by the spARCH component of the model, i.e., $\hat{\rho}=0.378$ and $\hat{\alpha}=0.020$. This can be interpreted as differences in the local uncertainty of the model. Hence, there are regions where the model predicts the true incidence rates more accurately, and there are regions with a worse fit. This can also be interpreted as local risks coming from unobserved, hidden factors. Note additionally that it is important to account for spatial conditional heteroscedasticity, as the estimates of spatial autoregressive models are biased if the error variance is not homogeneous across space. Inspecting the residuals, one can see that the spatial autocorrelation has been fully captured by the model, as Moran’s $I$ of the residuals is close to zero. In contrast, there is a weak spatial dependence in the squared residuals. To inspect the reason for this dependence graphically, the function plot() can be used to produce the plots shown in Figure 5.

After fitting the model, one also may include further regressors or estimate an intercept-only model via update(). For illustration, we added the percentage of positive results for a prostate-specific antigen (PSA) test in each county as an additional explanatory variable by

⬇

1R> out2 <- update(out, . ~ . + PSA_test)

The PSA test is used for prostate cancer screening, meaning that there should definitely be a positive dependence between the PSA test and the incidence rates. In fact, the estimated parameter is positive, and the AIC is lower compared to the previous model. To be precise, the updated parameters are

⬇

1 Estimate Std. Error t value Pr(>|t|)

2alpha (spARCH) 0.0199281 0.0043105 4.6231 3.78e-06 ***

3rho (spARCH) 0.3902185 0.1280266 3.0479 0.0023041 **

4lambda (SAR) 0.6643605 0.0366748 18.1149 < 2.2e-16 ***

5(Intercept) 1.1349551 0.2301554 4.9313 8.17e-07 ***

6X_factor_scores[, 2] 0.0198504 0.0069903 2.8397 0.0045159 **

7X_factor_scores[, 10] -0.0224035 0.0065828 -3.4034 0.0006656 ***

8PSA_test 0.0095962 0.0042728 2.2459 0.0247125 *

6 Summary and discussion

This paper examines spatial models for autoregressive conditional heteroscedasticity. In contrast to previously proposed spatial GARCH models, these models allow for instantaneous autoregressive dependence in the second conditional moments. Previous approaches only allowed for spatial dependence in the first temporal lag. However, these models are also captured by the spatial ARCH approach, since temporal dependence can be included by appropriate choices of the weighting matrix. In addition to discussing previously proposed models, we introduced a novel spatial exponential ARCH model, for which the probability density has been derived and maximum-likelihood estimators discussed.

In addition to this theoretical model, we focus on the computational implementation of all considered spatial ARCH models in the R-package spGARCH. In particular, the simulation and estimation has been demonstrated. Regarding maximum-likelihood estimation, a broad range of spatial models are implemented in the package. Furthermore, the spatial weights matrices, as well as the mean model, can easily be specified by the user, providing a flexible and easy-to-use tool for spatial ARCH models. All estimation functions return an object for class ‘spARCH’, for which several generic functions are provided, such as summary(), plot(), and AIC(). This setup also allows the use of the R-base functions, such as step() for stepwise model selection or update() for updating the results of different mean models. Eventually, the use of these functions are demonstrated by an empirical example, namely county-level incidence rates of prostate cancer.

In the future, the package should be extended for further spatial ARCH-type models. Along this vein, a class for model specifications should be added alongside the actual implementations via arguments for the fitting functions. In that way, the package can be aligned to common time series ARCH packages, such as the rugarch package. Furthermore, the package could benefit from robust estimation methods, another focus for future research.

7 Appendix

Proof of Theorem 1.

For this definition of $g_{b}$, one could rewrite $\ln\boldsymbol{h}$ as

[TABLE]

with

[TABLE]

Since $w_{ij}\geq 0$ for all $i,j=1,\ldots,n$, $\mathbf{W}$ is positive definite and it holds that

[TABLE]

Thus, the relation between $Y(\boldsymbol{s}_{1}),\ldots,Y(\boldsymbol{s}_{n})$ and $\varepsilon(\boldsymbol{s}_{1}),\ldots,\varepsilon(\boldsymbol{s}_{n})$ is given by (1) and (11). ∎

Proof of Corollary 1.

For $\rho\geq 0$, $b\geq 0$, and $w_{ij}\geq 0$ for all $i,j$, the inverse

[TABLE]

is a non-negative matrix. Thus,

[TABLE]

is positive for $\alpha>0$. ∎