A Multivariate Spatial and Spatiotemporal ARCH Model

This paper introduces a multivariate spatiotemporal autoregressive conditional heteroscedasticity (ARCH) model based on a vec-representation. The model includes instantaneous spatial autoregressive spill-over effects in the conditional variance, as they are usually present in spatial econometric applications. Furthermore, spatial and temporal cross-variable effects are explicitly modelled. We transform the model to a multivariate spatiotemporal autoregressive model using a log-squared transformation and derive a consistent quasi-maximum-likelihood estimator (QMLE). For finite samples and different error distributions, the performance of the QMLE is analysed in a series of Monte-Carlo simulations. In addition, we illustrate the practical usage of the new model with a real-world example. We analyse the monthly real-estate price returns for three different property types in Berlin from 2002 to 2014. We find weak (instantaneous) spatial interactions, while the temporal autoregressive structure in the market risks is of higher importance. Interactions between the different property types only occur in the temporally lagged variables. Thus, we see mainly temporal volatility clusters and weak spatial volatility spill-overs.


Introduction
In general, spatiotemporal processes can be represented as multivariate time series.However, when analysing spatial and spatiotemporal data, one has to account for one key difference compared to multivariate time series.Due to their spatial nature, geographical proximity between the observations induces instantaneous interactions between them.This is commonly known as Tobler's first law of geography: "everything is related to everything, but near things are more related than distant things" (Tobler 1970).This does not only apply to the mean behaviour of the data, but also their variance.Thus, spatiotemporal models should always allow for instantaneous spatial interactions.
In this paper, we introduce a multivariate spatial and spatiotemporal autoregressive conditional heteroscedasticity (spatial ARCH, briefly spARCH) model.Using a vector representation, we extend the spatial ARCH models of Otto et al. (2018); Sato and Matsuda (2021); Otto and Schmid (2019) to multivariate and spatiotemporal data.In that sense, the approach follows the same logic as classical time-series vec-ARCH models (cf.Engle and Kroner 1995).
Hence, we call the new multivariate, spatiotemporal ARCH process vec-spARCH.All these approaches trace back to the seminal papers of Engle (1982) and Bollerslev (1986).In contrast to previous multivariate spatiotemporal GARCH models (e.g., Borovkova and Lopuhaa 2012), we allow for instantaneous dependence over space at the same time point, which is important for spatial models.Thus, our multivariate vec-spARCH model can also be applied for purely spatial data, when there are spatial volatility clusters (i.e., clustered regions of high/low volatilities).Moreover, it is worth noting that, using this vec-representation, the above-mentioned spatial GARCH proposed by Otto and Schmid (2019) nests multivariate models, such that their results can also directly be applied.Alternative models that allow instantaneous spatial interactions in the variance are spatial stochastic volatility models, as proposed by Taspinar et al. (2021).
The vec-spARCH process distinguishes between three different effects: (1a) instantaneous spatial effects of the same variables, (1b) instantaneous cross-variable spatial effects, (2a) temporal autoregressive effects of the same variables, (2b) cross-variable temporal autoregressive effects, and (3) variable-specific unconditional volatility levels.Each of the effects is described by a parameter matrix or vector, for which we derive a quasi maximum-likelihood (QML) estimator.For estimation, a logarithmic transformation of the vec-spARCH is applied (cf.Robinson 2009), such that the model can be represented as multivariate spatiotemporal autoregressive model of the transformed quantity.The asymptotic consistency of the QML estimator has been shown by Yang and Lee (2017) for a multivariate spatial autoregressive model (i.e., without temporal dimension) and by Yu et al. (2008) for a univariate spatiotemporal autoregressive process.Under certain regularity conditions, which are commonly used in spatial econometrics, we show the identifiability and consistency of the estimators.
In practice, spatiotemporal ARCH models are particularly important, because an ARCH error process can also account for variation due to latent factors.In particular, for small spatial units, it is often difficult to quantify influential factors with the same spatial resolution.For instance, the average income of households in small spatial units, e.g.postal-code areas, does not necessarily reflects the economic power of these particular units, because people's daily cycles usually span across multiple small spatial units.That is, people do not necessarily live where they work or spend most of their time.In such cases, spatial and spatiotemporal ARCH models are important error distributions of any model to account for unobservable factors.
The remainder of the paper is structured as follows.Firstly, we introduce the multivariate modelling framework and discuss how the model can transformed to a regular univariate spatiotemporal process.Further, we derive the Gaussian logarithmic likelihood and show the asymptotic consistency for the QML estimator under several regularity assumption that are often met in spatial econometrics.Secondly, we analyse the finite-sample performance of the proposed estimator for several model specifications and two different error distributions, namely standard normal and heavy-tailed errors (t 3 -distributed).Thirdly, a real-world application is presented, for which we show that it is important to account for instantaneous spatial interactions and cross-variable correlations.To be precise, Berlin real-estate prices of three different property types are analysed and we find weak spatial interactions, even though they are dominated by the temporal effects.These interdependencies are more pronounced when the spatial units and time intervals are small.Finally, Section 5 concludes the paper with a summary and brief outlook to future research and potential fields of application.

Multivariate Spatiotemporal ARCH Model
In spatial statistics/econometrics, autoregressive spill-over effects are instantaneous.That is, no time lag is required for shocks to affect neighbouring locations.Instead, we assume that the conditional variance can vary over space depending on the realised variance at adjacent locations.This results in spatial clusters of high and low variances.For previous univariate or multivariate spatiotemporal GARCH models, such as proposed by Borovkova and Lopuhaa (2012); Hølleland and Karlsen (2020), spatial spill-overs could only occur after one time instance.In other words, the conditional variance at each locations depends on the past squared observations at the same location and its neighbours, but not on their neighbouring locations at the same time point.This is the fundamental difference between multivariate time series models covering spatiotemporal data and approaches from spatial statistics or econometrics.

Model specification
Assume that {Y t (s) ∈ R p : s ∈ D s ⊂ R q , t ∈ Z} is a p-variate spatiotemporal stochastic process in a q-dimensional space D s with positive volume (cf.Cressie and Wikle 2011).More precisely, the process is observed at n different locations s 1 , . . ., s n , i.e., at each location s i and time point t we observe a vector Y t (s i ) = (Y 1,t (s i ), . . ., Y p,t (s i )) .Moreover, let Y j,t = (Y j,t (s 1 ), . . ., Y j,t (s n )) the vector of the j-th characteristic at all locations and Y t = (Y 1,t , . . ., Y p,t ) be an n × p matrix of all observations (Y j,t (s i )) i=1,...,n,j=1,...,p at time point t.Suppose that the process is observed for t = 1, . . ., T .It is worth mentioning that a multivariate spatial log-ARCH model is present if T = 1 and a classical time-series log-ARCH models are also nested if D s is a singleton (i.e., n = 1).
Univariate spatial and spatiotemporal ARCH models have been introduced by Otto et al.
(2018) and Sato and Matsuda (2017).Moreover, Otto and Schmid (2019) generalised the model in a unified framework nesting spatial GARCH, E-GARCH, and Log-GARCH models.In this paper, we follow the idea of the symmetric spatial log-GARCH model of Sato and Matsuda (2021), which includes elements of GARCH and E-GARCH models, but does not coincide with one or the other even if D s consist of only a single location (i.e., the classical time series case).More precisely, the link function between the spatial volatility term is logarithmic like for E-GARCH models, while the volatility term depends on some transformation of the squared observed process (similar to GARCH models).In contrast to time-series models, in which the temporal lag is clearly defined by the past observations and future observations are not allowed to influence the current observation, there are complex interdependencies in spatial settings and there is no causal relation between the observations anymore.For instance, with two locations s 1 and s 2 (i.e., n = 2), location s 1 would influence s 2 at each time point and vice versa.In a time series context, this would correspond to a simultaneous influence from future and past values.Thus, for direct generalisation of GARCH or E-GARCH models like in Otto et al. (2018); Otto and Schmid (2019), difficult assumptions for the existence or invertibility of the process are required in the general case.In addition, existing software could directly be used with some adaptations for the spatiotemporal case (see Otto 2019).
The multivariate spatiotemporal ARCH model (vec-spARCH) is given by where ln Y (2) t denotes the matrix of squared observations ln Y 2 j,t (s i ) i=1,...,n,j=1,...,p , and ln H t is the matrix of all ln h j,t (s i ) with i = 1, ..., n rows and j = 1, ..., p columns.This matrix is the spatial equivalent of the conditional volatility (see Otto et al. 2019).Moreover, the matrix of disturbances is denoted by Ξ t = (ε 1,t , . . ., ε p,t ) with independent and identically distributed random vectors ε j,t with E(ε j,t ) = 0 and Cov(ε j,t ) = I.The weight n × n matrix W defines the spatial dependence structure, i.e., which locations are considered to be adjacent.Moreover, the cross-variable spatial effects are represented by the off-diagonal elements of Ψ, and the temporally lagged cross-variable effects are given by the off-diagonal elements of Π.Both matrices have dimension p×p.In addition, the own-variable spatial and temporal autoregressive ARCH effects are summarised by the diagonal entries of Ψ and Π, respectively.
Analogue to multivariate vec-ARCH time-series model of Engle and Kroner (1995), we can rewrite (1) to get the vectorised form (2) The Kronecker product is denoted by ⊗.Interestingly, using such vec-representation, one can see that the multivariate ARCH model is a special case of (univariate) np-dimensional spatial GARCH models with a weight matrix Ψ ⊗ W. Thus, also spatial GARCH and E-GARCH models can be constructed in the same manner and all results of Otto and Schmid (2019) can directly be applied.However, this will not be the focus of this paper.
Moreover, the multivariate spatiotemporal ARCH model can be written as multivariate spatiotemporal autoregressive process by applying a log-squared transformation, ln Y (2) Then, we get that ln Y (2) , the model can be rewritten as ln Y (2) Hence, the vec-spARCH model coincides with a multivariate spatiotemporal autoregressive process of the log-squared transformed process ln Y (2) t .For the multivariate but purely spatial case, Yang and Lee (2017) has derived conditions for identification and the consistency and asymptotic normality of a QML estimator.Furthermore, Yu et al. (2008) derive asymptotic results for a QML estimator of spatiotemporal but univariate process when both n and T are large.We combine these two results to propose a QML estimator for the spatiotemporal, multivariate ARCH model.

Assuming a standard normal error matrix
is the expectation of a log-Gamma distribution, i.e., E ln Ξ (2) t = γ − log(2) ≈ −1.27.Then, A can be determined from Ã, which facilitates the interpretation.With S np = I − Ψ ⊗ W, we can derive the sample log-likelihood for the spatiotemporal case with T time points, i.e., ln where σ 2 u is the variance of the transformed errors U t , which is known quantity in our case (otherwise A would not be identifiable).Furthermore, for standard normal Ξ t , we get σ 2 u = ψ(1/2) ≈ 4.93, where ψ denotes the trigamma function.It is worth mentioning that we derived the log-likelihood for multivariate Gaussian errors U t , which are in fact skewed because of the logarithmic transformation.In the following Section 2.2, however, we suppose much weaker conditions for the moments of U t , which are fulfilled in the case of standard normal Ξ t , for instance.Furthermore, let Ÿ t = ln vec(Y (2) t ) for an easier notation.With E( Ÿ t ) = S −1 np0 vec( Ã0 ) + Π 0 ⊗ Ÿ t−1 , we get the expected log-likelihood as

Assumptions and parameter space
Below, we discuss important model assumptions that are also needed to derive the asymptotic consistency of the QML estimators.
Assumption 1. Suppose that each element of Ξ t is not equal to zero with probability one for all t = 1, . . ., T .
To be able to apply the log-squared transformation of the observed process, we must assume that the response is not equal to zero with probability one.This is the case for any continuous error process Ξ t .In practice, due to missing values, there is sometimes an excess of zeros.In this case, often a small number is added to the zero values, such that the logarithmic transformation gets feasible (see, e.g., Francq and Zakoian 2011).If there are zero values with a probability larger than zero, Sucarrat and Escribano (2018) proposed an expectationmaximisation algorithm for estimation in the time-series case.This would be an interesting extension for future research.Further, we need some basic assumptions on the transformed error process U t to apply the results of Yang and Lee (2017) and Yu et al. (2008).
Moreover, the parameter space needs to be compact, as formulated in the following assumption, to prove the uniform convergence of the log-likelihood function.
Assumption 3. The parameter spaces for Ã, Ψ and Π are compact sets and all parameters in their interior generate a stable process.Moreover, the data-generating parameters Ã0 , Ψ 0 and Π 0 are in the interior of corresponding parameter space.
The key question of this assumption is the stability of the process.One could rewrite the model as Hence, the stability of the process does not only depend on the temporal parameter matrix Π but also on the weight matrix W (via S −1 np ).If the above series converges, we get a stable and stationary process.
Proposition 1.If all eigenvalues of S −1 np (Π ⊗ I) are smaller than one, the multivariate spatiotemporal ARCH process is stable across time.
Note that each stable spatiotemporal ARCH process is also weakly stationary.Furthermore, the boundary region of Ψ where |S np | = 0 can be problematic in practice.However, as long as the true parameter Ψ 0 is bounded away from this region, the maximisation algorithm will not get to these boundaries with a large probability (see also Yang and Lee 2017).
Assumption 4. The row and column sums of W in absolute values are uniformly bounded in n.Moreover, S np is invertible for all possible matrices Ψ in the parameter space and S −1 np is uniformly bounded in absolute row and column sums.
Assumption 4 is classical in spatial statistics to obtain a stable process across space (cf.Yang and Lee 2017; Kelejian and Prucha 1998;Lee 2004).Here, we could adopt the assumption as formulated in Yang and Lee (2017) for multivariate spatial autoregressive models.Together with Proposition 1, we obtain a stable process across space and time.In practice, the spatial weight matrices often standardised to meet these regularity conditions, e.g. the most widely adopted row-wise standardisation.
Assumption 5. Let n be a nondecreasing function of T and T → ∞.

Consistency of the QML estimator
Due to the presence of endogenous variables, i.e., the instantaneous spatial interactions, the identification of spatial models is generally more difficult than in the strict time-series case, where all spatiotemporal interaction may only occur after one time lag.Thus, we initially focus on the identification of the parameters which is needed for the asymptotic consistency of the QML estimator in the following Theorem 1.Since the identification is inherent with the spatial dimension of the model, we could follow the same strategy as in Yang and Lee (2017) for multivariate spatial autoregressive models.The identification is based on the information inequality, as proposed by Rothenberg (1971).
For the identification, we make use of the fact that the spatial dependence is constant across time and the temporal dependence is constant for all spatial locations.If either of them varies in space or time, further identifying information would be needed.Moreover, in contrast to Yang and Lee (2017), the errors are uncorrelated by definition and the error variance is supposed to be known.The assumption of an uncorrelated error process is essential for GARCH models for identification of the parameters in the conditional volatility equation, i.e., the so-called GARCH effects.Moreover, the assumption of a known error variance σ 2 u is of course restrictive (see also Francq and Zakoian 2011;Brockwell and Davis 2006) and it is often difficult to choose an appropriate value.For time series, ex-post scale adjustments have been proposed to circumvent this assumption (see Bauwens and Sucarrat 2010;Sucarrat et al. 2016).In this paper, however, we follow the classical approach and point to future research for these ex-post scale adjustments.Moreover, for the purely spatial case with T = 1, Ã0 must be constant across space to be identifiable.
To estimate the parameters, we propose a QML estimator based on the log-likelihood function given by (3).That is, the parameters ϑ = (vec( Ã) , vec(Ψ) , vec(Π) ) can be estimated by θnT = arg max ϑ∈Θ ln L(A, Ψ, Π|Y 0 ), where Θ is the parameter space that fulfils Assumption 3. It is worth noting that we need to condition on the observed vector at t = 0, Y 0 , because of the temporal autoregressive structure.
The asymptotic consistency of this QML estimator is summarised in the following theorem.

Monte Carlo Simulations
In the following section, we present the results of a series of simulations on the consistency of the parameters for finite samples.To give a first visual impression, we display a bivariate spatial ARCH process (T = 1, n = 900, Rook's continuity matrix) with and without spatial cross-correlation in Figure 1.For both examples, the spatial ARCH effects are equal to 0.5, a moderate level of spatial dependence.Therefore, spatial volatilities clusters can be seen in both cases.They are indicated by a higher variance, that is, more intensely coloured pixels, whereas clusters of low variance are close to zero indicated by evenly grey coloured pixels.Now, for the case with a cross-correlation of 0.35 (top panels), these clusters are aligned across the variables, while they are located at different positions in the lower panels with zero cross-correlation.
In our Monte Carlo simulation study, we simulated three different bivariate models (A, B, C) with two different error distributions (standard normal and t 3 ) with 1000 replications.For each combination, we successively increased the size of the spatial field n ∈ {25, 49, 100} and the length of the time series T ∈ {30, 100, 200}.We simulated the process on a two-dimensional grid as visualised in Figure 1 and the spatial weight matrix was chosen as row-standardised Queen's contiguity matrix.The data-generating parameters of the three considered models are as follows: (A) Spatiotemporal model with a weak spatial cross-correlation: For the first simulated model, i.e., Model A with standard normal errors, the parameter estimates are depicted as a series of boxplots for the three increasing sizes (n, T ) ∈ {(25, 30) , (49, 100) , (100, 200) } in Figure 2. In all cases, the asymptotic consistency of the QML estimator can be seen, because the boxplots are getting more centred around zero.Moreover, we see the typical bias of the QML estimators for small spatial fields, which rapidly vanishes with an increasing sample size.The same behaviour can be observed for all other settings and error distributions.The average bias and the root-mean-square errors (RMSE) are reported in Tables 1 and 2, respectively.Both the absolute values of the bias and the RMSE are approaching zero if n and T are increasing.
4 Real-World Example: Berlin Real-Estate Prices In the following section, we will show the application of the process to a real example.For this purpose, we model the changes in the average sales prices of undeveloped land, developed land and condominiums in Berlin.The data are average monthly average prices per square metre of   land or living space in each post-code region from 2002 to 2014.The average prices across all spatial locations are depicted in Figure 3 as time series process.
First of all, it must be noted that there are typically geographical dependencies in the housing market, unlike for other financial markets where trading can take place regardless of location.One of the most important factors in a purchase decision is the location of the property, whereby prices are also influenced by the surrounding neighbourhood.This dependency is in turn influenced by road connections, infrastructure or public transport.Furthermore, the price in the past plays also a role, as is typical for all time series.The temporal proximity creates a causal statistical dependence that decreases the further one looks into the past.These dependencies are observed both in the price process and in the risks in terms of price changes.
This motivates the application of the proposed multivariate spatiotemporal ARCH process to property sales returns.More precisely, we analysed the logarithmic, monthly returns of the average sales prices in each category for all n = 190 postcode regions in Berlin.The length of the time series is accordingly T = 156 and the process is p = 3-dimensional.To display the log-return process, Figure 4 shows the average log-returns across all locations in the temporal domain.Especially for the developed and undeveloped land, there were much fewer sales, such that the average returns are more volatile.In the case of no transactions in certain months and areas, we assumed that the average sales price did not change and, thus, the log-returns are zero.More precisely, we randomly simulated a normally distributed return with mean zero and standard deviation 0.0001 to not have positive probability for zero returns.In future, a more detailed analysis including a zero-transaction model would be interesting, especially for smaller time granularities and spatial locations.
The estimated parameters of the multivariate spatiotemporal ARCH process are reported in Table 3 along with their standard errors.The unconditional variance levels were assumed to be constant across space, but vary with the property types.Bearing in mind that we have modelled monthly returns, we observe interesting results.First, the spatial dependence is dominated by the temporal dependence that appears to be more important.Second, spatial spill-overs are positive (i.e., we observed clusters of higher variances/risks), but they are only significant for developed land.When increasing the temporal intervals from monthly to quarterly data, these spatial interactions will disappear.The same holds when grouping the spatial locations to   larger areas.This highlights the importance of spatial GARCH models for small spatial units and time granularities (as it is also well-known in finance).Third, cross-variable spill-overs are only significant at the first temporal lag (i.e., after one time period).More precisely, we see significant interactions only between developed and undeveloped land, but not for condominium prices.It is important to bear in mind that the spatial and temporal ARCH effects will also cover changes in the variance due to latent variables.Fourth, the unconditional variance varies across the property types with developed land experiencing the highest variance, followed by the condominium and undeveloped land.Note that undeveloped land usually does not have and will not get building permission.
able to model purely spatial data without the need of observations over time.In the empirical application, it gets obvious that the log-returns of several types of real estate are spatially autocorrelated.This indicate local clusters of increased volatilities and market risks -even though the temporal dependence appears to be more important.Thus, we could show that there are temporally and spatially varying volatilities.Furthermore, we found significant cross-variable dependence in the first temporal lags, but no significant instantaneous cross-variable interactions.This again motivates the application of a multivariate spatiotemporal ARCH model in such studies.
For this new model, we discussed the parameter estimation using a quasi-maximumlikelihood (QML) approach.For this reason, the process is reformulated in a vec-representation and a log-squared transformation is applied to obtain multivariate spatiotemporal autoregressive process.We showed the consistency of the QML estimator under regular assumptions for the error process when the spatial and temporal dimensions increase.In the finite-sample case, we could see rapidly decreasing root-mean-square errors (RMSEs) in a series of simulations with different model specifications and error distributions.All our simulations could be performed in a reasonable amount of time using a standard computer.The required computational resources are usually the bottleneck of the QML approach due to the computation of the log-determinant of the Jacobian matrix.
There are many further directions for future research and potential fields of applications.
First, we did only considered logarithmic structures in the volatility models, but no classical ARCH structures.However, since the multivariate spatiotemporal could be transformed to purely spatial models using the vec-representation, previous results of spatial ARCH and GARCH models could be applied.Furthermore, all these spatial econometric models rely on a (correctly) specified spatial weight matrix, which is, however, mostly unknown in practice.Hence, estimation methods for the entire spatial dependence structures (i.e., each spatial weight) are desirable from a practical perspective.Penalised methods seem to be promising in this case, because many links can be considered to be zero.
model without temporal dependence, but the same spatial dependence like for Model A: A 0 = 1 n 1 p , Ψ 0 = Spatiotemporal model with pronounced cross-correlation and weak spatial correlation, same temporal autocorrelation like for Model A: A 0 = 1 n 1 p , Ψ 0 =

Figure 2 :
Figure 2: Estimation performance of the QML estimator for Model A with standard normal errors.Top row: unconditional variance level a, centre: spatial coefficient matrix Ψ, bottom: temporal coefficient matrix Π.For each plot, we show the difference between the parameter estimate and the true data-generating parameter.

Figure 3 :
Figure 3: Monthly average prices (Euro/m 2 ) of undeveloped, developed land, and condominium across 190 post-code regions in Berlin from January 2002 to December 2014.

Table 1 :
Average bias of the QML estimates for all considered settings.

Table 2 :
Root-mean-square errors of the QML estimates for all considered settings.

Table 3 :
QML estimates and standard errors of the empirical example.Spatial ARCH effects are highlighted in light green, temporal ARCH effects in dark green.Significant effects are marked by an asterisk (* t-value > 1.9, ** t-value > 2).
Apart from applications in econometrics, also environmental and climate processes would be interesting and potential fields for application of the multivariate ARCH model.The process parameters can be interpreted as local risk measures, which is highly relevant in environmental studies.Furthermore, when considering our model as error process, spatially and temporally varying measurement or modelling uncertainties can be reflected in the statistical model.For example, this might be of interest for GNSS positioning in urban environments.Other fields, where local risks and cross-variable interactions are highly relevant, are epidemiological and medical studies.