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

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 logarithmic 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. Moreover, spatial GARCH models are briefly discussed. 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.


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 andLopuhaa 2012 andCaporin andParuolo 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). Further details and derivations can also be found in . 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 log-GARCH model by Pantula (1986); Geweke (1986); Milhøj (1987). Other recent papers propose a mixture of these two approaches (see Matsuda 2017, 2018b). Moreover, all these models can be used in spatiotemporal settings (see Sato and Matsuda 2018a).
In addition to the novel spatial logarithmic 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.2, we discuss all covered spatial and spatiotemporal ARCH-type models. In addition, we introduce a novel logarithmic 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 logarithmic spatial ARCH model. Furthermore, spatial GARCH models are briefly discussed. However, the focus of this paper should be on ARCH-type models. After these theoretical sections, we demonstrate the use of the R-package spGARCH in Section 2.5. 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 2.7 concludes the paper.

Spatial ARCH-type models
Let {Y (s) ∈ R : s ∈ D} 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 R q , as space is usually finite. For dealing with spatial lattice data, D is subset of the q-dimensional integers 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 R k × 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 Y = (Y (s 1 ), . . . , Y (sn)) at all locations s 1 , . . . , sn. For spatial ARCH models, we specify this vector as an analogue to the well-known time series ARCH models (cf. Engle 1982;Bollerslev 1986). However, note that the vector h does not necessarily coincide with the conditional variance as the variance in any location s j also depends on Y (s i ) for j = i (see Otto et al. 2018 for details). We now distinguish between several spatial ARCH-type models via the definition of h.

Spatial ARCH model
First, we define this vector h in such a way as to be analogous to the definition in . For this model, the vector h O is given by where diag(a) is a diagonal matrix with the entries of 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 α and ρ, whereas  proposed a more general model with a vector α = (α 1 , . . . , αn) and the first-order spatial lag Wdiag(Y )Y . For this definition, there is a one-to-one relation between Y and ε via the squared observations where W is a predefined spatial weighting matrix and Thus, 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. 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 . Suppose that α ≥ 0, ρ ≥ 0 and that det( It is important to note that A depends on both the weighting matrix and the realizations of the errors. In order to ensure that this condition is fulfilled,  propose to truncate the support of the error distribution on the interval (−a, a) with where || · || 1 denotes the matrix norm based on the Manhattan norm.
There are two cases in which the support of the errors does not need to be constrained. If ρ = 0, the process coincides with a spatial white noise process such that a equals ∞. Moreover, all entries of h are non-negative if W is similar to a strictly triangular matrix. Then, 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 ρ = 0. However, this support decreases with increasing values of ρ. For instance, if ρ = 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 ε(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 Y and their squares for ρ ∈ {0, 0.05, . . . , 2}. For the Monte Carlo simulation study, we simulate n = 400 observation on a two-dimensional lattice D = {s = (s 1 , s 2 ) ∈ Z 2 : 0 ≤ s 1 , s 2 ≤ 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 andBoots 1995, see also Cliff andOrd 1981). Thus, the asymptotic 95% confidence intervals are plotted in Figure 1, as well.

Spatial Log-ARCH model
Next, we consider an logarithmic spatial ARCH process (log-spARCH). In this setting, we define the natural logarithm of h E = (h E (s 1 ), . . . , h E (sn)) as with a function g b : R n → R n . Like Nelson (1991), we assume that for positive values of b. For this definition, there is a one-to-one relation between Y and ε, as we show in the following theorem.
At location s i , the value of h E (s i ) is then given by For this definition of g b , one could rewrite ln h as In contrast to the spARCH process described in Section 2.2.1, Corollary 1 shows that the entries of h E are positive for all ρ ≥ 0 and α > 0. Hence, the process is well-defined and there are no further restrictions needed, as in the case for the spARCH model. For all proofs, we refer to the Appendix.

Complex Spatial ARCH model
Now, we propose a complex-valued spARCH process. In order to obtain a solution of diag(h) 1/2 in the n-dimensional space of real numbers for the model defined in (2), all elements of the matrix (I − A 2 ) −1 must be nonnegative (see . For the complex spARCH process, we relax the assumption that there should be a solution to diag(h) 1/2 in the real numbers and also consider complex solutions. Thus, the definition of h coincides with h O of the original model, i.e.,

Spatiotemporal ARCH model
Finally, we show that spatiotemporal processes are covered directly by these approaches. For spatiotemporal data, the vector s simply includes both the spatial location ss and the point in time t, i.e., s = (ss, t) . 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 ≥ t i . However, the dimension of the weighting matrix W might become very large for this representation. More precisely, the matrix has dimension N T × N T , 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 W is usually sparse and could also have a block diagonal structure. Moreover, it is often reasonable to assume that h(s i ) is only influenced by the neighbors of s s,i at the same point of time and by past observations at the same location. Then the weighting matrix would have the

Process type
Definition of h Comments (2) ) ε ∼ MN(0, I), W must be a strictly triangular weighting matrix spatial log-ARCH ln h E = S (α1 + ρbW ln |Y |) ε ∼ MN(0, I), but moments of Y differ from the moments of classical spARCH process (cf. ) spARCH (complex) following structure Indeed, it is plausible to weight the spatial and temporal lags differently by replacing ρW by a sum with positive weights φ k for all temporal lags 1 ≤ k ≤ p.

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 Y , i.e., with u following either a spatial ARCH model with the original definition h O or the logarithmic model with h E . Thus, Further, we call this model the SARspARCH model. For λ = 0, the model collapses to a simple linear regression model; if, additionally, β = 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, λB can be interpreted directly as the spatial dependence of the process, while ρ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.

Generalized Spatial ARCH Models
Additionally, one may include spatially lagged observations of h to construct spatial GARCH-type models. For instance, a spatial GARCH model is given by whereW is a second spatial weighting matrix and λ is the corresponding spatial GARCH parameter. Obviously, the spatial GARCH-type models requires that (I − λW) is non-singular. In a similar manner, h LG can be specified as to define a spatial log-GARCH model. For theoretical details of spatial GARCH-type models, we refer to Otto and Schmid (2019) introducing a unified spatial GARCH model covering various spatial ARCH and GARCH models. Moreover, Otto and Schmid (2019) introduce an exponential spatial GARCH model allowing for asymmetry via an alternative definition of g in (5). To be precise, g is given by for the exponential spatial GARCH model.

Parameter Estimation
The parameters of a spatial ARCH process can be estimated by the maximum-likelihood approach.
To obtain the joint density for Y = k(ε), the Jacobian matrix of k −1 at the observed values y must be computed (e.g., Bickel and Doksum 2015). If fε is the distribution of the error process, then the joint density f Y of Y is given by If the residuals are additionally independent and identically distributed, the parameter estimates can be obtained from the maximization of the log-likelihood as follows The Jacobian matrix, of course, depends on the definition of h. For the spARCH process, this Jacobian matrix can be specified as In contrast, the Jacobian matrix for the log-spARCH process is slightly different, namely 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 for the spARCH process. Similarly, it is given by for the log-spARCH process, where • 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 W, and therefore S • W, are usually sparse. Thus, the required time for the estimation of the parameters depends mainly on the dimension and sparsity of W.
Certainly, the choice of the weighting matrices are an important design choice of the models, which has to be prespecified. However, the true structure of W is rarely known in practice. Moreover, all estimated parameters depend on the selection of this matrix. Hence, inference on these parameters and the coefficients itself must be interpreted based on the choice of W. For empirical applications, one might gain insights on the structure of W by looking at spatial autocorrelation functions or variograms. It is worth noting that the observations are uncorrelated for spatial ARCH models, so one should also look at squared observations. Then, the weighting scheme is typically chosen from a set of candidate schemes by maximizing certain goodness-of-fit criteria, like information criteria or out-of-sample prediction errors. For instance, W could be chosen as contiguity matrix, i.e., two locations are connected having positive weights, if they share a common border or if their distance is less than a certain threshold. For instance, in studies in spatial econometrics or epidemiology, the spatial domain is often a set of municipalities or counties (e.g., Amin et al. 2014;Buettner 2003). In this case, contiguity matrices are straightforward and if these binary matrices are additionally row-standardized, Wdiag(Y )Y can be interpreted as average of the squared neighboring observations. Alternatively, W can be specified as k-nearest-neighbor matrix, i.e., only the k nearest locations get positive weights, or as inverse-distance matrix, i.e., the weight between two locations is based on the distance between these locations. Further choices of W are discussed by . Finally, it is worthy to mention that the weights could also depend on exogenous variables or other factors. For instance, they could incorporate economic disparities, e.g., differences in the gross domestic products, poverty rates, household incomes etc., or other covariates, like the wind direction and speed when modeling spatial dependence of air pollutants (cf. Merk and Otto 2019). For spatiotemporal autoregressive processes, there are also some approaches to estimate the entire spatial dependence structure using machine learning methods (e.g., Lam and Souza 2016;Otto and Steinert 2018).

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.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 2.5.1) and estimation (Section 2.5.2) aspects of the spARCH, log-spARCH, and SARspARCH processes.

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 W is triangular) on a d × d spatial lattice with ρ = 0.7 and α = 1 using the following lines.

Function
Description Main functions sim.spARCH() Simulation of spARCH and log-spARCH processes sim.spGARCH() Simulation of spGARCH, E-spGARCH, and log-spGARCH 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 <-sim.spARCH(n = n, rho = rho, alpha = alpha, W = W, + type = "spARCH", control = list(seed = 5515)) To build the spatial weighting matrix, we used cell2nb() from the spdep package, returning an nb object of a d × 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, 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 = "spARCH", then the original spARCH process according to the definition in  is simulated.
-If there exists a permutation such that 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/ 4 ρ 2 ||W|| 1 .
• type = "log-spARCH", an log-spARCH process is simulated with an user-specified value of b (default 2) and standard normal random errors.
• type = "complex-spARCH", complex solutions of diag(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 ε is −20 −10 0 10 20 Above left: spatial white noise for comparison; center: oriented spARCH (type = gaussian ); right: Below left: spARCH with truncated normal errors (type = gaussian ); center: spatial E-ARCH (type = exp ), right: complex spARCH (type = complex ).

Figure 2:
Simulations on a two-dimensional lattice for triangular matrices (above) and nontriangular matrices (below). For all simulations, we set ρ = 0.7 and α = 1, and W is chosen to be the Queen's contiguity matrix.
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, 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.
As pointed out in Section 2.2, spatiotemporal ARCH models are directly covered if time is considered as one dimension of the q-dimensional space D. Thus, a two-dimensional spatiotemporal process {Y t (s) : t = 1, . . . , T ; s ∈ Ds} with Ds being a d × d unit grid and T = 20 points of time could be simulated as follows. Note that the true spatial orientation is preserved in this representation (4 × 4 grid), i.e., plots appearing close to each other are also located close to each other in space and they are, therefore, more related than more distant locations. An alternative representation as consecutive spatial random fields is shown in Figure  4.

Figure 3:
Simulated spatiotemporal log-ARCH process depicted as individual time series. The process has been simulated on a 4 × 4 spatial unit grid and for 20 points in time. The spatial ARCH parameter ρ equals 0.9 and α = 1.
The spatial weighting scheme has been chosen as block diagonal matrix like proposed above, i.e., constant matrices W 1 = W 2 = . . . = W T = Ws along the diagonal define the instantaneous spatial interactions, while constant weights of 0.2 below the diagonal define the extend of the temporal dependence. For instance, a central location would be weighted equally by all eight spatial neighbors with a weight of 0.125 and by the observation of the same location at the previous time point by 0.2. Thus, W describes the structure and the extend of the dependence in space and time. Eventually, ρ has been chosen as 0.8 and we simulated the process as logarithmic spatial ARCH process.
The resulting simulation is depicted in Figures 3 and 4. Whereas the simulated values are shown as time series plots placed at their correct spatial locations in Figure 3, 4 depicts the observations as consecutive spatial random fields. Note that the same color coding has been chosen for both representations. On the one hand side, one can observe spatial volatility clusters (e.g. in the pre-last plot in Figure 4, t = 19, in which the conditional variance is low in the upper left corner, whereas the conditional variance of the remaining locations is high). On the other hand, temporal volatility cluster can be observed as well. For example, at location 16, the variance is high at the first and last five time points, while it is lower between t = 6 and t = 14.
For sake of completeness, we briefly demonstrate the simulation of spatial GARCH-type models using the sim.spGARCH() function. Like for sim.spARCH(), the type of the spatial GARCH model t = 1 t = 2 t = 3 t = 4 t = 5 ...

t = 19 t = 20
The first five and the last two time points are plotted as spatial random fields, i.e., the simulations are shown in their natural temporal ordering. An alternative representation as time series in their true spatial ordering is shown in Figure 3.

Figure 4:
Simulated spatiotemporal log-ARCH process depicted as consecutive spatial random fields. The process has been simulated on a 4 × 4 spatial unit grid and for 20 points in time. The spatial ARCH parameter ρ equals 0.9 and α = 1.
can be chosen by the argument type. More precisely, there are the following options • type = "spGARCH" for simulation of spatial GARCH models according to the definition in (10), • type = "e-spGARCH" for simulation of exponential spatial GARCH models according to the definition in (12) with g, • type = "log-spGARCH" for simulation of logarithmic spatial GARCH models according to the definition in (12) with g b , and • type = "complex-spGARCH" for simulation of a complex-valued spatial GARCH model.
To simulate a spatial GARCH process, two spatial weights matrices need to be specified via the arguments W1 and W2. Moreover, two parameters ρ and λ are passed to the simulation function by the arguments rho and lambda. For instance, a spatial GARCH model can be simulated on a d × d spatial unit grid as follows R> require("spdep") R> rho <-0.5 R> lambda <-0.3 R> alpha <-1 R> d <-20 R> nblist <-cell2nb(d, d, type = "rook") # Rook s contiguity matrix R> W_1 <-nb2mat(nblist) R> W_2 <-W_1 R> Y <-sim.spGARCH(rho = rho, lambda = lambda, alpha = alpha, + W1 = W_1, W2 = W_2, type = "spGARCH") Similarly, spatial log-GARCH processes and exponential spatial GARCH processes can be simulated by changing the argument type. In this case, the parameters b must be provided for the log-GARCH or Θ and ζ for the e-spGARCH, respectively. These parameters can easily passed to sim.spGARCH() by the arguments b, theta, and zeta.

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 andYu 2012, 2010b,a). 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 log-spARCH processes, the likelihood functions of spARCH models are rather flat around the global maximum. This behavior is illustrated for simulated processes in Figure 5. The observations for the log-spARCH process have been simulated as follows.

R> nblist
<-cell2nb(20, 20, type = "queen") R> W <-nb2mat(nblist) R> y <-sim.spARCH(n = 20^2, rho = 0.5, alpha = 1, W = W, + type = "log-spARCH", control = list(seed = 5515)) To simulate an oriented process, the entries of W above the diagonal must be set to zero and the argument type must be changed to "spARCH", i.e., R> W[upper.tri(W)] <-0 R> y2 <-sim.spARCH(n = 20^2, rho = 0.5, alpha = 1, W = W, + type = "spARCH", + control = list(seed = 5515)) To estimate the parameters of an intercept-free log-spARCH model without any regressors, the formula passed to the function qml.spARCH() should be specified as y 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.  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 log-spARCH process would return the following results.

Residuals:
Min. The standard errors are estimated as Cramer-Rao bounds from the Hessian matrix of the loglikelihood function. For triangular weighting matrices, the estimators are asymptotically normally distributed . 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 6.  (1,3)) R> plot(spARCH_object) Reproduce the results as follows: eps <-residuals(x) W <-as.matrix(x$W) moran.plot(eps, mat2listw(W), zero.policy = TRUE, xlab = "Residuals", ylab = "Spatially Lagged Residuals") Reproduce the results as follows: eps <-residuals(x) W <-as.matrix(x$W) moran.plot(eps, mat2listw(W), zero.policy = TRUE, xlab = "Residuals", ylab = "Spatially Lagged Residuals") Reproduce the results as follows: eps <-residuals(x) std_eps <-(eps -mean(eps))/sd(eps) qqnorm(eps, ylab = "Standardized Residuals") qqline(eps) The mean equation can be specified as formula for all models, i.e., the spARCH, log-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 a + b, more sophisticated models can also be fitted, e.g., models with interactions y a + b:c, factor models y factor, polynomial models y poly(a,3), seasonally or regularly varying models of the form y sin(t) + cos(t) or y 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.

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). To reproduce the example, the logarithmic incidence rates as well as several covariates are included in the package.
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 P M 2.5 , P M 10 , SO 2 , N O 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. Note that the factor scores are directly included in the dataset prostate_cancer. Eventually, the final explanatory factors were chosen by minimizing the Bayesian information criterion using the generic function step() as follows.
R> data(prostate_cancer) R> out <-step(qml.SARspARCH(formula, B = B, W = W, type = "spARCH", + data = prostate_cancer), k = log(length(Y))) The formula object simply defines a linear model between the logarithmic incidence rates and all factors. Further, matrix B describes the predefined spatial dependence structure in the mean equation. For this analysis, 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, 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 nd and 10 th factor has been selected. Whereas the 2 nd factor has positive loadings mainly for fine particulate matters, P M 2.5 and P M 10 , the 10 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.
Call: qml.SARspARCH(formula = formula, B = B, W = W, type = "spARCH", Function formula type Resulting model qml.spARCH() y 0 "spARCH" spARCH model (see (1) and (2)) qml.spARCH() y 1 "spARCH" spARCH model with an additional intercept for the mean equation qml.spARCH() y a + b "spARCH" Linear Regression with regressors a and b and spARCH residuals qml.spARCH() y a + b:c "spARCH" Linear Regression with more complex expressions and spARCH residuals qml.spARCH() y 0 "log-spARCH" log-spARCH model (see (1) and (5)) qml.spARCH() y 1 "log-spARCH" log-spARCH model with an additional intercept for the mean equation qml.spARCH() y a + b "log-spARCH" Linear Regression with regressors a and b and log-spARCH residuals qml.spARCH() y a + b:c "log-spARCH" Linear Regression with more complex expressions and log-spARCH residuals qml.SARspARCH() y 0 "spARCH" SAR model without an intercept, but with spARCH residuals (see (8) and (9)) qml.SARspARCH() y 1 "spARCH" SAR model with an intercept and spARCH residuals qml.SARspARCH() y a + b "spARCH" SAR model with an intercept and the regressors a and b and spARCH residuals qml.SARspARCH() y a + b:c "spARCH" SAR model with more complex expressions and spARCH residuals qml.SARspARCH() y 0 "log-spARCH" SAR model without an intercept, but with log-spARCH residuals (see (8) and (9)) qml.SARspARCH() y 1 "log-spARCH" SAR model with an intercept and log-spARCH residuals qml.SARspARCH() y a + b "log-spARCH" SAR model with an intercept and the regressors a and b plus log-spARCH residuals qml.SARspARCH() y a + b:c "log-spARCH" SAR model with more complex expressions and log-spARCH residuals  First, we see that the model has a significant spatial autocorrelation in the mean equation sinceλ (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.,ρ = 0.378 andα = 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 7. 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 prostatespecific antigen (PSA) test in each county as an additional explanatory variable by R> 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.

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.   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.
In addition to discussing previously proposed models, we introduced a novel spatial logarithmic 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.