Spatial and Spatiotemporal Volatility Models: A Review

Spatial and Spatiotemporal Volatility Models: A Review

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

{\it Keywords:} Survey, volatility, GARCH models, stochastic volatility, spatial and spatiotemporal dependence, nonlinear models.

\spacingset{1.45}

Introduction and motivation

When observations of a random process have a natural ordering, such as the temporal ordering for time series or geographical locations for spatial processes, they are typically dependent. For time series, observations close together in time will be more closely related than observations further apart. Similarly, for spatial processes, fisher1935design stated that “patches in close proximity are commonly more alike, as judged by the yield of crops, than those which are further apart”, or more commonly known as Tobler's first law of Geography: “Everything is related to everything else, but near things are more related than distant things” Tobler70. The similarity may be reflected in the mean or trend behaviour, which motivates autoregressive processes and integrated processes, but also in the process variation or volatility, motivating autoregressive heteroscedasticity models or stochastic volatility models. Spatial and spatiotemporal data models must fulfil this key property: instant and direct dependence due to spatial or temporal proximity.

Various types of models have been introduced in time series analysis and spatial analysis to describe the dependence over time and space. The most popular approaches are based on linear structures where, e.g., the value of the time series (spatial process) at a certain time point (point in space) is a linear combination of preceding or neighbouring observations and possible regressors. Such types of processes, e.g. SARIMA models, have been analysed in time series analysis in detail brockwell2009time,brockwell2002introduction. Additionally, over the last few decades, the field of econometrics has seen the development of numerous linear spatial models, which have proven to be highly practical and useful anselin1988spatial, Lee:2004, Lesage:2009, KP:2010, elhorst2010applied, Elhorst:2014.

Unfortunately, these linear approaches are only of limited use to describe a dependence behaviour in the process' volatility. This point was discussed in detail by Engle82,Bollerslev86 for temporal processes. They introduced a new type of nonlinear temporal process, which has turned out to be extremely useful for modelling returns of financial data. They are based on a multiplicative decomposition and are, thus, nonlinear in the white noise process. Their main advantage is that they can describe a time-varying behaviour of the conditional variance, the so-called volatility. In general, there are different definitions of volatility andersen2010parametric. For instance, for ARCH and GARCH models, the volatility term coincides with the conditional variance of the response variable. For other volatility time series models, such as stochastic volatility models, the conditional variance is non-explicit francq2019garch. The time-varying volatility makes these nonlinear models quite attractive in practice since the risk behaviour of the process may change over time. Linear time series, e.g., ARMA processes, do not possess this property; they have a constant volatility. For an overview of nonlinear time series models, we refer the interested reader to turkman2016non or fan2003nonlinear.

In time series analysis, generalised autoregressive conditional heteroscedasticity (GARCH) and stochastic volatility models are widely used to model time-varying volatility Francq11,francq2019garch. However, many real-world phenomena are observed at multiple geo-referenced locations and, thus, exhibit spatial or spatiotemporal dependence, which means that the volatility of neighbouring series may influence the volatility of a series. Spatial and spatiotemporal volatility models have been recently developed to account for this spatial dependence. After first being mentioned “as a byproduct” in Bera04, Otto16_arxiv introduced spatial ARCH models (published in Otto18_spARCH,otto2019stochastic), and Sato17 introduced spatial log-ARCH models. Interestingly, both models were developed simultaneously and independently of each other to be able to represent the ARCH-like dependencies in the variance of spatial processes. Moreover, Robinson:2009 and tacspinar2021bayesian introduced spatial stochastic volatility models that include an additional stochastic term in the volatility process.

Since then, various nonlinear spatial and spatiotemporal volatility models with spatial and spatiotemporal dependence in the (conditional) variance have been proposed. However, there has been little comparison among these models (e.g., OttoSchmid19_arxiv_unified compared some models in a unified framework). In this paper, we aim to provide a comprehensive review of the recent literature on spatial and spatiotemporal volatility models. As these models resemble their time series counterparts, we begin by briefly reviewing important aspects of discrete time series volatility models. Subsequently, we describe spatial and spatiotemporal volatility models considered in the literature. Besides motivating each model, we describe important estimation strategies that are applicable in the spatiotemporal context, e.g., quasi-maximum-likelihood (QML) approach, generalised methods-of-moments (GMM) approach, and the Bayesian MCMC sampling schemes. We describe possible extensions and provide directions for future research.

The remainder of the paper is structured as follows. In Section (ref), we start with discrete time series volatility models and their multivariate extensions, which are also applicable to spatiotemporal data. In Section (ref), we discuss spatial GARCH properties (i.e., instantaneous GARCH-type interactions across the cross-sectional or spatial domain) and compare different specifications that have been proposed in the literature. In Section (ref), our focus is on spatiotemporal volatility models allowing for these instantaneous spatial GARCH-type interactions in addition to the temporal GARCH structure. Further extensions and related models are discussed in the ensuing Section (ref). Finally, in Section (ref), we conclude the review with an outlook on future research and a discussion of open problems. Some of the technical details are collected in an appendix.

Time series volatility models

This section will briefly describe two classes of time series volatility models: (i) ARCH and GARCH-type models and (ii) stochastic volatility models. These models provide the basis for all spatial and spatiotemporal ARCH/GARCH and stochastic volatility models developed in the literature. It is well known that the volatility of many financial time series, such as asset returns and exchange rates, changes over time. For example, Figure (ref) displays a univariate time series of daily logarithmic returns of the Dow Jones Index from 2011 to 2023. The bottom plot shows measures of volatility using absolute returns and a 30-day moving average, which indicates temporal variations in volatility levels. These variations can be useful for identifying patterns or trends in the data and informing investment strategies or risk management decisions.

Both classes of models aim to describe the dependence in the conditional variance, i.e., volatility, which should be separated from the mean behaviour of the random process, such as temporal trends, seasonality, or autoregressive dependence. For ARCH and GARCH models, the volatility is a function of the squares of previous observations and past conditional variances (i.e., volatilities). In the case of stochastic volatility models, the volatility is modelled through a latent stochastic process, also depending on the temporally lagged squared observations. In Sections (ref) and (ref), we briefly describe important univariate and multivariate versions suggested in the literature for both classes of models. Throughout this section, we consider the discrete time series $\{ \boldsymbol{Y}_t \in \mathbb{R}^r : t \in \mathbb{Z}\}$, i.e., an $r$-dimensional random process $\boldsymbol{Y}_t$ with equidistantly ordered temporal observations.

ARCH and GARCH models

Univariate ARCH and GARCH models

In the univariate case, $r = 1$, the random process $Y_t$ is given by

where $h_t^{1/2}$ is a scaling factor and $\{ \varepsilon_t \}$ is a sequence of independent and identically distributed (i.i.d.) variables with mean $0$ and variance $1$. Moreover, $h_t$ depends on the past realisations of $Y_t$, and it coincides with the conditional variance of $Y_t$. Thus, it can be interpreted as the volatility of the series. This idea traces back to the seminal work of Engle82. To be precise, the volatility process is defined as

for an ARCH($p$) process with unknown parameters $\alpha_0, \ldots, \alpha_p$, which have to be estimated. To obtain a positive conditional variance, the constant term $\alpha_0$ is assumed to be positive and $\alpha_i \geq 0$ for all $i = 1, \ldots, p$. In practice, large lag orders $p$ are often needed to capture the volatility dynamics, e.g., those of inflation indices engle1983estimates. A GARCH process extends this model by allowing for autoregression and moving average components in the conditional variance equation, reducing the required lag orders and, thus, the computational burdens. The basic GARCH($p$,$q$) model is given by

with $\beta_i \geq 0$, $i = 1, \ldots, q$, being additional parameters to be estimated Bollerslev86. A GARCH process is weakly stationary if $\sum_{i=1}^p \alpha_i + \sum_{j=1}^q \beta_j < 1$. For a detailed discussion of the statistical properties of ARCH and GARCH models and the parameter estimation, we refer the interested reader to the textbook of francq2019garch.

ARCH and GARCH models are also particularly useful as error processes of other regression and time series models to describe time-dependent and dynamic model uncertainties. Since they have an expectation of zero, they can be directly applied to any other mean model, as Engle82 already demonstrated for an ARCH regression model with a linear regression component. Later, weiss1984arma,weiss1986asymptotic,weiss1986arch considered ARCH models for the errors of autoregressive moving average processes. These models exhibit an autoregressive dependence in both the mean and the conditional variance.

Since then, several extensions and adaptations of GARCH models have been proposed. Each of these models has its own strengths and weaknesses, and choosing the best model for a particular data set depends on the nature of the data, the research question, and other factors. First, to avoid the non-negative constraints of the model coefficients, Geweke86,Pantula86,Milhoj87 proposed a logarithmic expression of the (log-)volatility equation, i.e.,

Consequently, the process exhibits multiplicative dynamics in the volatility, while the standard ARCH and GARCH models have additive volatility dynamics. Moreover, this logarithmic GARCH (log-GARCH) model allows for a direct transformation of the log-squared observations to an autoregressive moving average process of orders $p$ and $q$. Interestingly, the log-GARCH model is invariant to power log-GARCH models where $\log h_t$ is replaced by $\log h_t^\delta$ with the power $\delta > 0$, which acts as a scaling factor of the logarithmic terms sucarrat2019log. Since $h_t = \text{exp}(\alpha_0 + \sum_{i = 1}^{p} \alpha_i \log Y_{t-i}^2 + \sum_{i = 1}^{q} \beta_i \log h_{t-i})$, the conditional variance is always positive for real-valued coefficients $\{\alpha_i:i=1,\ldots,p\}$ and $\{\beta_i:i=1,\ldots,q\}$. This is particularly interesting when exogenous covariates influence volatility. That is,

with $\delta_j$ being the corresponding regression coefficients of the $j$-th covariate $x_{j,t}$ sucarrat2019log. In the GARCH counterpart, the regressors are assumed to be almost surely positive with non-negative coefficients francq2019qml.

Second, to replicate the so-called leverage effect, which is often observed in financial return data, the model has been extended to Exponential GARCH models (EGARCH, Nelson91) that allow for asymmetric dependence by including the sign of the residuals in the log-volatility equation:

Third, higgins1992class proposed a non-linear ARCH (NARCH) model, which nests the linear ARCH model of Engle82 as a special case and converges to the log-ARCH model in the limiting case. More precisely, the volatility equation of the NARCH model is given by

with $\phi_i \geq 0$ for $i = 0,1,\ldots, p$ and $\delta > 0$. This model converges to a log-ARCH model for $\delta$ approaching zero.

Apart from these models, several other extensions have been proposed, which should only briefly be mentioned, e.g., threshold GARCH models (TGARCH, zakoian1994threshold,glosten1993relation, adding a threshold term to the GARCH equation for different volatility dynamics below and above the threshold), Glosten, Jaganathan and Runkle GARCH models (GJR-GARCH, asymmetric volatility by including both positive and negative residuals in the GARCH equation), or fractionally integrated GARCH models (FIGARCH, long memory in the volatility process by using fractional integration in the conditional variance equation). For a more detailed overview of univariate time-series ARCH and GARCH models, we refer the interested reader to the survey paper of bera1993arch and the textbook of francq2019garch.

Multivariate ARCH and GARCH models

In general, financial asset returns tend to move together over time. Figure (ref) shows exemplarily a multivariate time series consisting of daily logarithmic returns of three selected stocks in the Dow Jones Index, Procter & Gamble (PG), 3M (MMM), and International Business Machines (IBM). The top plot displays the raw data of the daily returns over time, and the bottom plot shows the 30-day moving averages of the squared daily returns, indicating the co-movement of volatility among the three stocks. This can provide insights into the correlations and dependencies between the stocks, which can be useful for portfolio and risk management. Thus, a multivariate framework is well-suited for modelling the time-varying conditional covariance matrix for all returns. Several multivariate GARCH (MGARCH) models have been proposed to account for such dependence in the volatilities.

Recall that $\boldsymbol{Y}_t$ is an $r$-dimensional vector of returns at time point $t$. It is assumed that

where $\{ \boldsymbol{\varepsilon}_t \}$ is a sequence of i.i.d. $r$-dimensional random variables with $\operatorname*{E}( \boldsymbol{\varepsilon}_t ) = {\boldsymbol 0}$ and $\text{Var}(\boldsymbol{\varepsilon}_t) = \mathbf{I}_r$. The square root has to be understood in the sense of the Cholesky factorization, i.e. $\mathbf{H}_t^{1/2}$ is the unique symmetric and positive symmetric matrix with $\mathbf{H}_t^{1/2} (\mathbf{H}_t^{1/2})^\prime = \mathbf{H}_t$. The matrix $\mathbf{H}_t$ is allowed to depend on certain parameters, and it is assumed to be a measurable function with respect to the $\sigma$-algebra generated by $\boldsymbol{Y}_v, v < t$. Thus,

i.e., $\mathbf{H}_t$ is the volatility matrix of $\boldsymbol{Y}_t$.

Most of the published papers on this topic appeared at the end of the 80s and the 90s of the 20th century. In principle, all these models only differ in modelling the matrix $\mathbf{H}_t$. Here we want to sketch some of the most relevant approaches briefly.

The first paper on MGARCH models is due to bollerslev1988capital. They introduced the so-called vector GARCH model, briefly VEC model. It is based on the idea of transforming the matrix $\mathbf{H}_t$ to a vector by using the vech-operator harville1998matrix. For a square $r \times r$ matrix $\mathbf{B}$, the operator $\text{vech}(\mathbf{B})$ is defined as the $r(r+1)/2$-dimensional vector obtained by stacking the columns of the lower triangular part of $\mathbf{B}$. Let $\boldsymbol{h}_t = \text{vech}(\mathbf{H}_t)$. Then, it holds for a VEC-GARCH(1,1) process that

where $\boldsymbol{\eta}_{t-1} = \text{vech}(\boldsymbol{\varepsilon}_{t-1} \boldsymbol{\varepsilon}_{t-1}^\prime)$, $\mathbf{A}$ and $\mathbf{B}$ are assumed to be $r(r+1)/2 \times r(r+1)/2$ parameter matrices, and $\boldsymbol{\omega}$ is a $r(r+1)/2$ parameter vector. The total number of parameters of this model is $r(r+1)(r(r+1)+1)/2$. For $r=2$ it is $21$, for $r=3$ it is $78$, and for $r=4$ it is already $210$. Thus, the number of parameters increases very fast as $r$ increases. This is the reason why in practice, the model is only used for small values of $r$, e.g., $r=2$ or $r=3$.

Another possibility to reduce the number of parameters is to assume that $\mathbf{A}$ and $\mathbf{B}$ are diagonal matrices bollerslev1988capital. This model is denoted as the diagonal vector GARCH model, briefly DVEC-GARCH(1,1) model. Using the Hadamard product harville1998matrix, it can be written as

where $\mathbf{A}^\#$, $\mathbf{B}^\#$, and $\mathbf{\Omega}^\#$ are $r \times r$-matrices implied by $\mathbf{A} = \text{diag}(\text{vech}(\mathbf{A}^\#))$, $\mathbf{B} = \text{diag}(\text{vech}(\mathbf{B}^\#))$, and $\mathbf{\Omega} = \text{diag}(\text{vech}(\mathbf{\Omega}^\#))$. Consequently, $h_{ij,t}$ depends only on its own lag and on $\varepsilon_{i,{t-1}} \varepsilon_{j,{t-1}}$. This restriction dramatically simplifies the parameter estimation; however, it is still not suitable for large-scale systems.

A further simplification of this model was discussed by ding2001large. They choose $\mathbf{A}^\#$ and $\mathbf{B}^\#$ as matrices of rank one or as multiples of the matrix consisting purely on the element $1$. Further, riskmetrics1996jp uses the exponentially weighted moving average model, which can be written as

which corresponds to a scalar VEC-GARCH(1,1) model. Riskmetrics recommends choosing the factor $\lambda$ equal to $0.94$ for daily data and $0.97$ for monthly data bollen2015should,gonzalez2007optimality.

Strong restrictions on the parameters are necessary to ensure the matrix $\mathbf{H}_t$ to be positive definite. For that reason, engle1995multivariate proposed another approach, the BEKK (Baba, Engle, Kraft, and Kroner) model. The BEKK-GARCH(1,1,$K$) process is given by

where $\mathbf{A}_i$, $\mathbf{B}_i$, and $\mathbf{\Omega}$ are $r \times r$ matrices and $\mathbf{\Omega}$ is positive definite. The BEKK-GARCH model is a special case of the VEC-GARCH approach, but the converse is not true stelzer2008relation. The number of parameters is again high. Similar to the VEC-GARCH and DVEC-GARCH, the application of the BEKK-GARCH model reduces to cases where $r$ is small.

Bollerslev90 introduced a type of MGARCH model where the conditional correlations are constant (CCC model). For the CCC-GARCH(1,1) process, it holds that

where $\mathbf{R} = ( \rho_{ij} )$ and $\mathbf{D}_t = \text{diag}(h_{11,t}^{1/2},..., h_{rr,t}^{1/2})$ and

Here, the number of parameters is $r(r+5)/2$.

A model with dynamic conditional correlations (DCC model) was proposed by Engle02 and tse2002multivariate. For the DCC-GARCH(1,1) model of Engle02 it holds that

with $\mathbf{D}_t$ as above and

and $\boldsymbol{u}_t = ( u_{i,t} )$ with $u_{i,t} = \varepsilon_{i,t}/\sqrt{h_{ii,t}}$. $\bar{\mathbf{Q}}$ denotes the unconditional covariance matrix of $\boldsymbol{u}_t$ and $\alpha$ and $\beta$ are non-negative numbers satisfying that $\alpha + \beta < 1$.

Besides these models, many further proposals have been made. The above models seem to be the most applied ones. There are also other attempts to overcome the problem of dimensionality. Factor GARCH models are one step in that direction. Here the idea is that some factors drive the behaviour of the stock returns. Such an approach was discussed by, e.g., engle1990asset, lin1992alternative and bollerslev1993common.

We refer to the overview papers by bauwens2006multivariate and silvennoinen2009multivariate, as well as the book by francq2019garch, where the presented multivariate models and many further ones are discussed in more detail.

Stochastic volatility models

Univariate stochastic volatility models

In stochastic volatility models, the volatility process is modelled through a latent stochastic process. Although it is difficult to determine the exact origin of these models as they arose from various research efforts addressing different issues, Taylor:1982, Taylor:1986 seem to be the first to consider a univariate discrete version that can be considered as an alternative to the ARCH process. To learn more about the origin and development of stochastic volatility models, refer to Eric:1996 and Shephard:2005. A standard discrete time stochastic volatility model is specified in the following way:

where $Y_t$ is the observed response variable, $\{h_t\}$ is the sequence of the unobserved log-volatility, assuming an AR(1) process with a mean parameter $\mu_h$ and an autoregressive parameter $|\phi|<1$. The model includes two independent disturbance terms denoted by $\varepsilon_t$ and $u_t$. The sequence $\{\varepsilon_t\}$ includes the independent random variables with an identical distribution with mean $0$ and variance $1$. The disturbance term $\{u_t\}$ in the log-volatility equation are independent and have identical distribution with mean $0$ and variance $\sigma^2_h$. In this specification, the sign of $Y_t$ is determined by that of $\varepsilon_t$, and the volatility clustering and fat tail properties observed in the marginal distribution of $Y_t$ are delivered by the time-varying log-volatility (see Eric:1996 on the statistical properties of stochastic volatility models). This standard model can also be obtained as a discrete-time approximation to various diffusion processes in the continuous-time asset pricing literature Hull:1987, Wiggins:1987, Melino:1990, Scott:1989.

The outcome and log-volatility equations can be modified to formulate alternative specifications. The log-volatility equation is specified as an AR(1) process and can be generalised to any ARMA process. For example, if $h_t$ follows an AR(p) process, then it will take the following form:

where $\phi_1,\hdots,\phi_p$ are unknown autoregressive parameters.

An alternative specification can be obtained by assuming the presence of infrequent jumps in the outcome equation. Adding a jump component to the outcome equation can improve the fit of the observed time series of returns because the jump component may capture outliers as well as asymmetry in the return distribution Andersen:2002, Chib:2002. The stochastic volatility model with a jump component can be specified as

where $q_t$ is the jump random variable, and $k_t$ is the jump size random variable. The jump random variable is a Bernoulli random variable with success probability $P(q_t=1)=\kappa$, and the jump size is modelled as $\log(1+k_t)\sim N(-0.5\delta^2,\delta^2)$. In this model, $\kappa$ and $\delta$ are additional unknown parameters that we need to estimate along with $\mu_h$, $\phi$, and $\sigma^2_u$.

Another variant can be defined by allowing the volatility feedback in the outcome equation Koopman:2002:

where the scalar unknown parameter $\alpha$ gives the effect of volatility on the outcome variable. Chan:2017 extended this model by considering time-varying parameters:

where $\boldsymbol{x_t}$ is the $k\times1$ vector of covariates with matching time-varying parameter vector $\boldsymbol{\beta}_t$. This model generalises the model suggested in Koopman:2002 by allowing time-varying parameters $\boldsymbol{\beta}_t$ and $\alpha_t$ in the outcome equation. Let $\boldsymbol{\gamma}_t=(\alpha_t,\boldsymbol{\beta}^{'}_t)^{'}$ be the $k\times1$ vector of time-varying parameters. Chan:2017 assumes an a random walk process for $\boldsymbol{\gamma}_t$ such that $\boldsymbol{\gamma}_t=\boldsymbol{\gamma}_{t-1}+\boldsymbol{\nu}_t$, where $\boldsymbol{\nu}_t\sim N(\boldsymbol{0},\mathbf{\Gamma})$, and $\mathbf{\Gamma}$ is the $(k+1)\times(k+1)$ covariance matrix.

Harvey:1996 consider a variant that allows for the so-called “leverage effect” via introducing correlation in the disturbance terms of the outcome and the log-volatility equations. See also Eric:2004, Yu:2005 and Omori:2007 on the different versions of this model. Shephard:2005 notes that Hull:1987 were the first to propose a continuous-time stochastic volatility model incorporating the leverage effect. This work, in turn, inspired the development of the EGARCH model proposed by Nelson (1991) Shephard:2005. The variant proposed by Harvey:1996 can be specified as

where $\varrho$ is the correlation parameter. In this specification, Yu:2005 defines the leverage effect as the negative relationship between $E(h_t|Y_t)$ and $Y_t$, and derived the following equation:

This result suggests that this specification will exhibit the leverage effect whenever $\varrho<0$.

Another variant can be obtained by assuming a scale mixture distribution for the outcome variable Eric:2004, Chib:2002. This version takes the following form:

where $\{\omega_t\}$ is the sequence of latent variables that are independent and have identical distribution. Under the assumptions that $\varepsilon_t\sim N(0,1)$ and $\omega_t\sim IG(\nu/2,\nu/2)$, where $IG$ denotes the inverse gamma distribution, it can be shown that the marginal distribution of $\omega^{1/2}_t\varepsilon_t$ (unconditional on $\omega_t$) is the standard $t$ distribution with $\nu$ degrees of freedom Geweke:1993. Omori:2007 consider the same model under the assumption that $\log(\omega_t)\sim N(-0.5\tau^2,\tau^2)$, where $\tau$ is a scalar unknown parameter with $\tau^2\sim\text{Gamma}(1,1)$.

Chan:2013 introduces a class of models that includes both the moving average and stochastic volatility components that can nest a variety of specifications as special cases. This specification takes the following form:

where $\mu_t$ is the time-varying conditional mean process, $\psi_1,\hdots,\psi_q$ are the unknown moving average parameters, and the disturbance terms $v_t$ and $u_t$ are independent of each other for all leads and lags. Let $\boldsymbol{\mu}=(\mu_1,\hdots,\mu_T)^{'}$, $\boldsymbol{h}=(h_1,\hdots,h_T)^{'}$ and $\boldsymbol{\psi}=(\psi_1,\hdots,\psi_T)^{'}$. Then, this specification gives

which indicates that the conditional variance of $Y_t$ is a moving average of $q+1$ most recent variances $e^{h_t},e^{h_{t-1}},\hdots,e^{h_{t-q}}$. Moreover, unlike the standard stochastic volatility model in (ref) and (ref), $Y_t$ is serially correlated even after conditioning on $\boldsymbol{h}$. Chan:2013 shows that

where $\psi_0=1$. Thus, the conditional covariances are also time-varying because of the presence of the log-volatility $h_t$. As stated in Chan:2013, popular specifications can be obtained from this model by choosing a suitable conditional mean process $\mu_t$. For example, some of these specifications are (i) an AR(p) model: $\mu_t=\beta_0+\beta_1Y_{t-1}+\hdots+\beta_pY_{t-p}$, (ii) a linear regression model: $\mu_t=\boldsymbol{x}^{'}_t\boldsymbol{\beta}$, where $\boldsymbol{x}_t$ is the $k\times1$ vector of covariates, (iii) an unobserved component model: $\mu_t=\tau_t$, where $\tau_t=\tau_{t-1}+\varepsilon^{\tau}_t$ and $\varepsilon^{\tau}_t\sim N(0,\sigma^2_{\tau})$, (iv) a time-varying regression model: $\mu_t=\boldsymbol{x}^{'}_t\boldsymbol{\beta}_t$, where $\boldsymbol{\beta}_t=\boldsymbol{\beta}_{t-1}+\boldsymbol{\varepsilon}^{\boldsymbol{\beta}}_t$ and $\boldsymbol{\varepsilon}^{\boldsymbol{\beta}}_t\sim N(\boldsymbol{0},\mathbf{\Sigma}_{\boldsymbol{\beta}})$.

The likelihood function of a stochastic volatility model is a mixture over the distribution of $\mathbf{h}$ and therefore requires evaluation of a high-dimensional integral. For example, the likelihood function of the model in (ref)-(ref) can be defined as $f(\boldsymbol{Y}|\boldsymbol{\theta})=\int f(\boldsymbol{Y}|\boldsymbol{h},\boldsymbol{\theta})f(\boldsymbol{h}|\boldsymbol{\theta})\text{d}\boldsymbol{h}$, where $\boldsymbol{Y}=(Y_1,\hdots, Y_T)^{'}$, $f(\boldsymbol{h}|\boldsymbol{\theta})$ is the prior distribution of $\boldsymbol{h}$ determined by (ref), and $\boldsymbol{\theta}=(\mu_h,\phi,\sigma^2_h)^{'}$.\footnote{Note that $f(\boldsymbol{Y}|\boldsymbol{\theta})$ is also called the observed-data likelihood function or the integrated likelihood function. Two other functions that can be defined are (i) the conditional likelihood function $f(\mathbf{Y}|\mathbf{h},\boldsymbol{\theta})$, and (ii) the complete-data likelihood function given by $f(\boldsymbol{Y},\boldsymbol{h}|\boldsymbol{\theta})=f(\boldsymbol{Y}|\boldsymbol{h})\times f(\boldsymbol{h}|\boldsymbol{\theta})$. The conditional and complete-data likelihood functions are readily available for all stochastic volatility models.} This feature indicates that the estimation based on the exact maximum likelihood is not readily available and poses difficulties for likelihood-based estimation procedures. In the literature, various estimation methods exist, including the generalised method of moments (GMM), the quasi-maximum likelihood method, spectral GMM based on the characteristic function, indirect inference methods, Monte Carlo maximum likelihood methods, and Markov chain Monte Carlo (MCMC) based methods. According to Asai:2006, the choice of an estimation method can be based on the following properties: (i) efficiency, (ii) estimation of volatility, (iii) optimal filtering, smoothing, and forecasting methods, (iv) computational efficiency, and (v) applicability to flexible models. Among others, see JPR:1994, Kim:1998, and Broto:2004 on the estimation methods suggested in the literature for the univariate stochastic volatility models.

Multivariate stochastic volatility models

As in the case of multivariate ARCH and GARCH models described in Section (ref), the cross-dependence among the volatility of different financial assets leads to modelling volatility in a multivariate framework, which can lead to efficient estimation. Let $\boldsymbol{h}_t=(h_{1t},\hdots,h_{rt})^{'}$ be the $r\times 1$ vector of log-volatility terms, and $\mathbf{H}^{1/2}_t=\text{diag}\left(e^{h_{1t}/2},\hdots,e^{h_{rt}/2}\right)$ be the $r\times r$ diagonal matrix with the $i$th diagonal element $e^{h_{it}/2}$. Harvey:1994 consider the following multivariate version:

where $\boldsymbol{\mu}$ and $\boldsymbol{\phi}$ are the $r\times1$ vectors of unknown coefficients, $\circ$ denotes the Hadamard product, $\mathbf{P}_{\boldsymbol{\varepsilon}}=(\rho_{ij})$ is a positive definite correlation matrix with $\rho_{ii}=1$ and $|\rho_{ij}|<1$ for $i,j=1,\hdots,r$, and $\mathbf{\Sigma}_{\boldsymbol{u}}=(\sigma_{u,ij})$ is the $r\times r$ positive definite covariance matrix. Harvey:1994 also consider the multivariate $t$ distribution for $\boldsymbol{\varepsilon}_t$. The volatility process can be generalised to a VARMA structure in the following way:

where $\mathbf{\Phi}(L)=\left(\mathbf{I}_r-\sum_{i=1}^p\boldsymbol{\phi}_i\circ L^i\right)$ and $\mathbf{\Theta}(L)=\left(\mathbf{I}_r-\sum_{i=1}^p\boldsymbol{\theta}_i\circ L^i\right)$, $\mathbf{I}_r$ is the $r\times r$ identity matrix, $L$ is the lag operator, and $\{\boldsymbol{\phi}_i\}$ and $\{\boldsymbol{\theta}_i\}$ are vectors of parameters.

Note that if the off-diagonal elements of $\mathbf{P}_{\boldsymbol{\varepsilon}}$ and $\mathbf{\Sigma}_{\boldsymbol{u}}$ are zeros, then this multivariate model simply specifies a univariate standard stochastic volatility model for each component of $\boldsymbol{Y}_t$. Although the non-zero off-diagonal elements of $\mathbf{\Sigma}_{\boldsymbol{u}}$ introduce correlation across the volatility terms, the model does not allow for the leverage effect. To introduce the leverage effect, Asai:2006 consider the following model:

where $\mathbf{L}=\text{diag}\left(\lambda_1\sigma_{\boldsymbol{u},11},\hdots,\lambda_r\sigma_{\boldsymbol{u},rr}\right)$. Thus, the model exhibits the leverage effect when $\lambda_i<0$ for $i=1,\hdots,r$.

Some other alternative versions that may not lead to the leverage effect as defined by Yu:2005 are also considered by Dani:1998 and Chan:2006. See Yu:2006 for further details on the models that can deliver the leverage and asymmetric effects. There are also alternative specifications in the literature, including parsimonious specifications based on additive and multiplicative factors structures, time-varying correlation matrix models, matrix exponential models, Cholesky decomposition-based models, Wishart models, and range-based models. The details of these models and the estimation approaches considered in the literature are surveyed in Yu:2006, Renate:2006, and Chib:2009.

Spatial volatility models

Now, suppose that the process is observed across space. In contrast to time series, where the index is a scalar value, the index of each observation is now at least two-dimensional. Consider the random process $\{ \boldsymbol{Y}(\boldsymbol{s}) \in \mathbb{R}^r : \boldsymbol{s} \in D \subseteq \mathbb{R}^d, d > 1\}$, where $D$ is the spatial domain, typically a subset of $\mathbb{R}^d$ with $d > 1$ and positive volume. For instance, if $D \subseteq \mathbb{Z}^d$, the spatial domain would be called lattice (e.g., satellite image sequences for $d = 2$, or CT images for $d = 3$), while a continuous spatial process is present if $D \subseteq \mathbb{R}^d$ (e.g., soil samples, or species distributions). Another typical example is the case when $D$ is a discrete set $\{ \boldsymbol{s}_1, \ldots, \boldsymbol{s}_n\}$ of locations or polygons (e.g., air quality measurement stations, economic country/county-level data). Furthermore, $D$ could be considered to be a spherical space, $\mathbb{S}^d = \{ \boldsymbol{x} \in \mathbb{R}^{d+1} : || \boldsymbol{x} || = c\}$ with the radius $c$, which is particularly useful for modelling global data on the Earth. Figure (ref) presents an example of a purely spatial process where the monthly log returns of condominium prices in Berlin are shown on the map. Each location, postcode region, is represented by a polygon precisely defined by geographical coordinates.

Spatial ARCH and GARCH models

Spatial ARCH models were first mentioned “as a byproduct” in Bera04. Later, Otto18_spARCH and otto2019stochastic introduced the first purely spatial ARCH model (jointly published on arXiv in Otto16_arxiv), which has an ARCH-like dependence structure in the conditional variances. The first spatial ARCH models were introduced for univariate processes and purely spatial domains, i.e., the process is observed only once for all locations in $D$. In other words, there is only one time point $t$, and we do not repeatedly observe the process over time. For the univariate case with $r = 1$ and $n$ spatial locations, a spatial ARCH model is defined by

where $h(s_i)$ is a scaling factor of the error $\varepsilon(\boldsymbol{s}_i)$ of the $i$-th location, analogously to the time-series case. The errors are supposed to be i.i.d. across all spatial locations and have a zero mean and a constant variance of 1. Further, $h(\boldsymbol{s}_i)$ depends on all adjacent realisations of the response variable $Y$, where the adjacency is defined by a so-called spatial weights matrix $\mathbf{W} = (w_{ij})_{i,j = 1, \ldots, n}$. The weights $w_{ij}$ are non-zero if $h(\boldsymbol{s}_i)$ might be influenced by $h(\boldsymbol{s}_j)$, i.e., $\boldsymbol{s}_i$ is in proximity to $\boldsymbol{s}_j$. Typically, the definition of the weights matrix depends on the geographical space and the coordinates of all locations. For example, $w_{ij}$ can be chosen as the inverse-distance between $\boldsymbol{s}_i$ and $\boldsymbol{s}_j$, or $w_{ij}$ could be equal to $1/k$ for all $k$ nearest neighbours. Then,

Due to this temporal simultaneity, the volatility $h(\boldsymbol{s})$ at locations $\boldsymbol{s}$ given all other locations is difficult to interpret because all other locations simultaneously depend on $Y(\boldsymbol{s})$. Thus, the interpretation of $h(\boldsymbol{s})$ is slightly different compared to time-series ARCH models, and $h(\boldsymbol{s}_i)$ does not coincide with the conditional variance at location $\boldsymbol{s}_i$ given all neighbouring observations. Thus, we refer to $\{ h(\boldsymbol{s}_i) \}$ as the volatility process. For parameter estimation, Otto18_spARCH proposed a quasi-maximum-likelihood estimator, which is computationally implemented in the R-package spGARCH Otto19_RJournal.

Moreover, higher-order spatial dependence can be considered by

where $\{\mathbf{W}_k = (w_{k,ij})_{i,j = 1, \ldots, n} : k = 1, \ldots, p \}$ is a set of suitable weight matrices, e.g., separating the influence for different directions (i.e., $k = 1$ corresponds to the northward direction, $k = 2$ to the eastward direction, and so on). In contrast to time-series models, higher-order spatial dependence is often directly included in the first spatial lag by choosing $w_{ij}$ to be positive also for larger lag-orders or distances between $\boldsymbol{s}_i$ and $\boldsymbol{s}_j$. For instance, typical choices of $w_{ij}$ are inverse-distance-based weights, i.e., $w_{ij} = d(\boldsymbol{s}_i, \boldsymbol{s}_j)^{-k}$ where $k$ controls weight decay across space and $d(\boldsymbol{s}_i, \boldsymbol{s}_j)$ is a suitable metric to measure the distance between $\boldsymbol{s}_i$ and $\boldsymbol{s}_j$.

Alternatively, a spatial autoregressive term of the volatilities can also be included. In this way, a spatial GARCH model of order $(1,1)$ can be defined as

which was introduced by otto2022general. This model was first discussed within a unified framework proposed by OttoSchmid19_arxiv_unified. By allowing for higher-order spatial lag terms, we obtain the following volatility equation of spatial GARCH models

In contrast to time-series models that often employ the natural one-way ordering (i.e., past observations can only influence future observations, but not vice versa), spatial models must allow for two-sided influence. Thus, there is typically no causal order between the observations and further assumptions are needed for the existence of a real-valued process. Furthermore, if the spatial locations $\boldsymbol{s}$ are one-dimensional (i.e., $d = 1$), the spatial ARCH models coincide with the time-series ARCH model by Engle82, where $\mathbf{W}$ acts like a backward-shift operator. More precisely, if the locations are ordered and equidistantly spaced, $\mathbf{W}$ would be a sparse matrix with ones on the first subdiagonal for an ARCH(1) process. Similarly, an ARCH($p$) process can obtained if $\mathbf{W}_k$ has ones on the $k$-th subdiagonal for $k = 1, \ldots, p$. Note that such models with a triangular weight matrix lead to directional spatial processes merk2021directional. In these cases, a causal ordering of the locations exists.

Below, we again focus on the special case of spatial ARCH(1) and GARCH(1,1) processes with two \textcolor{blue}{possibly different} weight matrices $\mathbf{W}_1$ and $\mathbf{W}_2$ for the ARCH and GARCH term, respectively. Like for time-series GARCH models, also spatial GARCH models require additional assumptions on the parameters and/or weight matrix to ensure the non-negativity of $h(\boldsymbol{s}_i)$ for all $i = 1, \ldots, n$. To analyse this in more detail, let $\boldsymbol{h} = (h(\boldsymbol{s}_i))_{i = 1, \ldots, n}$ and $\boldsymbol{Y}^{(2)} = (Y(\boldsymbol{s}_i)^2)_{i = 1, \ldots, n}$ be the $n$-dimensional vectors of all $h(\boldsymbol{s}_i)$ and the squares of $Y(\boldsymbol{s}_i)$, respectively. Then, (ref) can be written in a matrix notation as

For a spatial GARCH model, the volatility is specified in a matrix notation as follows

or in the reduced form, i.e.,

Furthermore, let $\mathbf{A} = \text{diag}(\varepsilon(\boldsymbol{s}_1)^2, \ldots, \varepsilon(\boldsymbol{s}_n)^2)\mathbf{W}_1$. For the existence of a real-valued process, i.e., $h(\boldsymbol{s}_i) \geq 0$ for all $i = 1, \ldots, n$, we need to ensure that (i) the inverse $\mathbf{S}(\beta_1) = \left(\mathbf{I} - \beta_1 \mathbf{W}_2\right)^{-1}$ exists, and (ii) the inverse of $\tilde{\mathbf{S}}(\alpha_1) = (\mathbf{I} - \alpha_1 \mathbf{A}^2)^{-1}$ exists and all elements are non-negative. The first condition is a typical assumption in spatial econometrics, and several choices of the weights matrix and the corresponding parameter space have been discussed in the literature, e.g., if $\mathbf{W}_2$ is a row standardised weights matrix and $\beta_1 \in [0,1]$ the inverse exists. Generally, the spectral radius of $\beta_1 \mathbf{W}_2$ must be smaller than one. The second condition is more complicated to check in practice. On the one side, it can be guaranteed by limiting the squared error terms $\{\varepsilon(\boldsymbol{s}_i)^2\}$, e.g. by assuming a truncated error distribution otto2019stochastic. On the other hand, it can be ensured by considering the characteristics of underlying spatial dynamics implied by $\mathbf{W}_1$. In several specific cases, the condition is always fulfilled, e.g., in the case of directional processes in which $\mathbf{W}_1$ can be expressed as a strictly triangular matrix merk2021directional.

By contrast, in the logarithmic setting, weaker restrictions on the parameter space are needed for the non-negativity of $\boldsymbol{h}$, but they also imply a different volatility structure. Sato17 introduced a spatial log-ARCH process, for which the (log-)volatility equation is given by

With $\log \boldsymbol{h} = (\log h(\boldsymbol{s}_1), \ldots, \log h(\boldsymbol{s}_n))'$ and $\log \boldsymbol{Y}^{(2)} = (\log Y(\boldsymbol{s}_1)^2, \ldots, \log Y(\boldsymbol{s}_n)^2)'$ being the vectors of the element-wise logarithms of $\boldsymbol{h}$ and the squared observation $\boldsymbol{Y}^{(2)}$, the model can be expressed in a matrix notation as

Using the log-squared transformation Robinson:2009, the model can be transformed into a spatial autoregressive model of the log-squared observations. The transformed errors $\{\log \varepsilon(\boldsymbol{s}_i)^2\}$ of this spatial autoregressive representation follow a log $\chi_1^2$ distribution. Thus, they no longer have zero means and are heavily left-skewed. For the existence of the process, the regular conditions of spatial autoregressive models apply, i.e., the inverse $(\mathbf{I} - \alpha_1 \mathbf{W})^{-1}$ must exist.

As an extension of the spatial log-ARCH model, su2023statistical proposed the following expression of the log volatility

with $F(\alpha) = \mathbf{I}_n - A(\alpha)$ and

The coefficients $a_l(\alpha)$ are real-valued deterministic functions with $A(0) = \mathbf{I}_n$. In this general form, different kinds of spatial dependence structures can be considered, e.g., spatial autoregressive structures with $A(\alpha) = \mathbf{I}_n - \alpha \mathbf{W}_1$ (i.e., the log-ARCH model of Sato17), spatial moving average structures with $A(\alpha) = (\mathbf{I}_n - \alpha \mathbf{W}_1)^{-1}$, or matrix exponential structures with $A(\alpha) = e^{\alpha \mathbf{W}_1}$. Besides, su2023statistical consider exogenous regressors influencing the log volatility. For this reason, the constant term $\alpha_0\boldsymbol{1}_n$ should be replaced by $\mathbf{X}\boldsymbol{\delta}$ with an $n \times (k+1)$-dimensional matrix $\mathbf{X}$ and $\boldsymbol{\delta} = (\delta_{0}, \ldots \delta_{k})'$ being a vector of linear regression coefficients.

Furthermore, the log-ARCH model can be generalised to a spatial log-GARCH model. For that reason, Takaki:2021 defined the log-volatility process as follows

Moreover, Takaki:2021 also allowed for regressive effects on the volatility, which can easily be included by adding $\boldsymbol{x}(\boldsymbol{s}_i)' \boldsymbol{\delta}$ in (ref), where $\boldsymbol{x}(\boldsymbol{s}_i) = (x_1(\boldsymbol{s}_i),\hdots,x_k(\boldsymbol{s}_i))^{'}$ is a vector of exogenous variables at location $\boldsymbol{s}_i$.

Dogan2023bayesian consider a higher-order version of the spatial log-GARCH model suggested in Takaki:2021, where the scaling factors in (ref) follow

for $i=1,2,\hdots,n$. Here, $p$ and $q$ are two finite positive integers, and $(w_{1r,ij})_{r=1}^p$ and $(w_{2l,ij})_{l=1}^q$ are non-stochastic weights matrices with zero diagonal elements. The corresponding $\{\alpha_{1r}\}_{r=1}^p$ and $\{\beta_{1l}\}_{l=1}^q$ are unknown scalar parameters. For the estimation, Dogan2023bayesian transform (ref) by taking the square of both sides and then taking its natural logarithm (i.e., the log-squared transformation). Then, the transformed outcome equation can be written as

where $Y^*(\boldsymbol{s}_i) = \log Y(\boldsymbol{s}_i)^2$, $h^{*}(\boldsymbol{s}_i) = \log h(\boldsymbol{s}_i)$, and $\operatorname{\varepsilon}^*(\boldsymbol{s}_i) = \log\operatorname{\varepsilon}(\boldsymbol{s}_i)^2$. Note that $\operatorname{\varepsilon}^*(\boldsymbol{s}_i)$ has a $\log\chi^2_1$ distribution with the density

with $\operatorname*{E}(\operatorname{\varepsilon}^*(\boldsymbol{s}_i))\approx-1.2704$ and $\text{Var}(\operatorname{\varepsilon}^*(\boldsymbol{s}_i))=\pi^2/2\approx4.9348$. This density function is highly skewed with a long left tail, as visualised in Figure (ref).

Let $\boldsymbol{Y}^*=(Y^*(\boldsymbol{s}_1),\hdots,Y^*(\boldsymbol{s}_n))^{'}$, $\boldsymbol{h}^{*}=(h^{*}(\boldsymbol{s}_1),\hdots, h^{*}(\boldsymbol{s}_n))^{'}$ and $\boldsymbol{\operatorname{\varepsilon}}^*=(\operatorname{\varepsilon}^*(\boldsymbol{s}_1),\hdots,\operatorname{\varepsilon}^*(\boldsymbol{s}_n))^{'}$. Then, the higher-order spatial GARCH model in Dogan2023bayesian can be written as

where $\mathbf{W}_{1r}=(w_{1r,ij})$ and $\mathbf{W}_{2l}=(w_{2l,ij})$ are the $n\times n$ weights matrices and $\mathbf{Z}=(\boldsymbol{z}(\boldsymbol{s}_1),\hdots,\boldsymbol{z}(\boldsymbol{s}_n))^{'}$ is the $n\times k$ matrix of exogenous variables. Let $\mathbf{S}(\boldsymbol{\beta}_1)=(\mathbf{I}_n - \sum_{l=1}^q\beta_{1l}\mathbf{W}_{2l})$, where $\mathbf{I}_n$ is the $n\times n$ identity matrix and $\boldsymbol{\beta}_1=(\beta_{11},\hdots,\beta_{1q})^{'}$. Also, let $\boldsymbol{\alpha}_1=(\alpha_{11},\hdots,\alpha_{1p})^{'}$. Under the assumption that $\left\Vert \sum_{l=1}^1\beta_{1l}\mathbf{W}_{2l}\right\Vert < 1$, $\mathbf{S}(\boldsymbol{\beta}_1)$ is invertible Horn:2013. Then, the reduced form of equation (ref) is given by

Substituting this equation into (ref) and rearranging yield

Let $\mathbf{G}(\boldsymbol{\theta})= \left(\mathbf{I}_n - \sum_{r=1}^p\alpha_{1r} \mathbf{W}_{1r} - \sum_{l=1}^1\beta_{1l}\mathbf{W}_{2l}\right)$ and $\boldsymbol{\theta}= (\boldsymbol{\alpha}_1^{'},\boldsymbol{\beta}_1^{'})^{'}$. Then, under the assumption that $\left\Vert \sum_{r=1}^p\alpha_{1r} \mathbf{W}_{1r} + \sum_{l=1}^q\beta_{1l}\mathbf{W}_{2l}\right\Vert<1$, we obtain

When the weights matrices are row normalised, $\left\Vert \sum_{r=1}^p\alpha_{1r} \mathbf{W}_{1r} + \sum_{l=1}^q\beta_{1l}\mathbf{W}_{2l}\right\Vert<1$ simplifies to $\sum_{r=1}^p|\alpha_{1r}| + \sum_{l=1}^q|\beta_{1l}| < 1$ by the triangle inequality.

Takaki:2021 propose a Gaussian pseudo maximum likelihood estimator for their spatial log-GARCH model by approximating the distribution of the transformed error terms with a normal distribution. They show that the resulting likelihood estimator attains the standard large sample properties. However, it is well known in the time series literature that the Gaussian pseudo maximum likelihood estimator obtained in this way might have poor finite sample properties because the normal approximation to the distribution of the log-squared error terms provides a poor approximation JPR:1994, Shephard:1994, Kim:1998, Koopman:1998.

Dogan2023bayesian instead propose approximating the distribution of $\operatorname{\varepsilon}^{*}(\boldsymbol{s}_i)$ using a mixture of Gaussian distributions and develop a Bayesian estimation algorithm. They assume that the distribution of $\operatorname{\varepsilon}^{*}(\boldsymbol{s}_i)$ can be approximated by the following $10$-component Gaussian mixture distribution Omori:2007:

where $\varphi(\operatorname{\varepsilon}^*(\boldsymbol{s}_i)|\mu_j,\,\sigma^2_j)$ denotes the Gaussian density with mean $\mu_j$ and variance $\sigma^2_j$, and $c_j$ is the probability of $j$-th mixture component. The parameters of the $10$-component Gaussian mixture distribution are given in Table (ref). A comparison between the $10$-component Gaussian mixture distribution and the normal distribution in approximating the distribution of $\operatorname{\varepsilon}^{*}(\boldsymbol{s}_i)$ is illustrated in Figure (ref). It is evident that the $10$-component Gaussian mixture distribution provides a very accurate approximation, whereas the normal distribution offers a poor approximation. \definecolor{LightCyan}{rgb}{0.88,1,1}

To complete the model specification, Dogan2023bayesian assume the following independent prior distributions: $\alpha_{1r}\sim \text{Uniform}(-1,\,1)$ for $r=1,\hdots,p$, $\beta_{1l}\sim \text{Uniform}\left(-1,\,1\right)$ for $l=1,\hdots,q$, and $\boldsymbol{\delta}\sim N(\boldsymbol{\mu}_{\delta},\mathbf{V}_{\delta})$, where $\text{Uniform}(-1,\,1)$ is the uniform distribution over the unit interval. A detailed description of the proposed Gibbs sampler can be found in Algorithm (ref) in the Appendix.

Spatial stochastic volatility models

In this section, we review extensions of the standard stochastic volatility models in the time series literature to spatial data. Our starting point will be the stochastic volatility model considered in Robinson:2009 and tacspinar2021bayesian, where the log-volatility terms are modelled through a first-order spatial autoregressive process. The resulting spatial stochastic volatility model shares similar properties with the standard stochastic volatility model in time series, and it is designed to capture volatility clustering in spatial data.

tacspinar2021bayesian specify the spatial stochastic volatility (SSV) model as

for $i=1,\hdots, n$, where the first equation is the outcome equation and the second equation is the log-volatility equation (or the state equation). In the outcome equation, $h(\boldsymbol{s}_i)$ denotes the latent log-volatility term, and $\varepsilon(\boldsymbol{s}_i)$ is an i.i.d. normal random variable with mean zero and unit variance. In the log-volatility equation, $\mu_h$ is the constant mean parameter, and $u(\boldsymbol{s}_i)$ is an i.i.d. normal random variable with mean zero and variance $\sigma^2_u$. The $w_{ij}$'s are the non-stochastic spatial weights such that they are zero when $i=j$. These elements represent the degree of spatial association between the log-volatility terms. The scalar parameter $\phi$ is the spatial autoregressive parameter and provides a measure of spatial correlations among $h_i$'s. This model can be considered a spatial extension of the stochastic volatility model in (ref)-(ref).

Let $\boldsymbol{h}=(h(\boldsymbol{s}_1),\hdots,h(\boldsymbol{s}_n))^{'}$ be the $n\times1$ vector of log-volatilities and $\boldsymbol{u}=(u(\boldsymbol{s}_1),\hdots,u(\boldsymbol{s}_n))^{'}$ be the $n\times1$ vector of error terms. Also, let $\mathbf{W} = (w_{ij})_{i,j = 1, \ldots, n}$ be the $n\times n$ non-stochastic matrix for the spatial weights. Then, the state equation can be written in vector form as

where $\boldsymbol{1}_n$ is the $n\times1$ vector of ones. Define $\mathbf{B}(\phi)=(\mathbf{I}_n - \phi \mathbf{W})$. Under some restrictions on the parameter space of $\phi$, $\boldsymbol{h} = \mu_h\boldsymbol{1}_n + \mathbf{B}^{-1}(\phi)\boldsymbol{u}$ exists, where $\mathbf{I}_n$ is the $n\times n$ identity matrix.\footnote{The necessary and sufficient condition for the invertibility of $\mathbf{B}(\phi)$ is that the spectral radius of $\phi \mathbf{W}$ must be less than $1$. See Lee:2004, Lesage:2009, KP:2010, and Elhorst:2014 for a discussion on the parameter space of spatial parameters.}

The conditional variance of $Y(\boldsymbol{s}_i)$ varies over the relevant space, because $\text{Var}\left(Y(\boldsymbol{s}_i)|h(\boldsymbol{s}_i)\right)=e^{h(\boldsymbol{s}_i)}$. Furthermore, $ \operatorname*{E}(Y(\boldsymbol{s}_i)Y(\boldsymbol{s}_j))=0$ for all $i\ne j$ implying that $Y(\boldsymbol{s}_i)$'s are not spatially correlated. Let $\boldsymbol{k}_i(\phi)$ denote the $i$th row vector of $\mathbf{B}(\phi)$ placed in a column vector, and let $r\in\mathbb{N}$ be an even number. Then, it follows that

where $\mu_{r}=\frac{r!}{2^{r/2}\times\left(r/2\right)!}$. Therefore, $\operatorname*{E}(Y(\boldsymbol{s}_i)^4)/\left[\operatorname*{E}(Y(\boldsymbol{s}_i)^2)\right]^2 - 3 = 3\left(e^{\sigma^2_u\Vert \boldsymbol{k}_i(\phi)\Vert^2} - 1\right)>0$. Hence, $Y(\boldsymbol{s}_i)$ has a leptokurtic symmetric distribution. Moreover,

This covariance is generally not zero unless $\phi=0$. Thus, the higher moments of $Y(\boldsymbol{s}_i)$'s are correlated, implying that $y(\boldsymbol{s}_i)$'s are spatially dependent.

For the estimation, it is more convenient to turn the spatial stochastic volatility model into a linear state-space model, and to this end, both Robinson:2009 and tacspinar2021bayesian transform the outcome equation by taking the square of both sides and then taking its natural logarithm (i.e., the log-squared transformation). Then, the outcome equation can be written as

as for spatial GARCH models. Robinson:2009 approximates the distribution of $\operatorname{\varepsilon}^{*}(\boldsymbol{s}_i)$ with the normal distribution and proposes a Gaussian pseudo maximum likelihood estimator for the estimation. The pseudo maximum likelihood estimators obtained in this manner may attain the standard large sample properties, but they tend to have poor finite sample properties, as mentioned above. Alternatively, tacspinar2021bayesian propose to approximate the distribution of $\operatorname{\varepsilon}^{*}(\boldsymbol{s}_i)$ using a mixture of Gaussian distributions \textcolor{blue}{as in (ref)} so that the resulting estimation system turns into a linear Gaussian state space model. They then use the data augmentation technique to facilitate the Bayesian estimation by treating $\boldsymbol{h}$ as an additional parameter vector. The Bayesian MCMC algorithm is described in detail in Algorithm (ref) shown in the Appendix.

Spatiotemporal volatility models

In this section, we discuss spatiotemporal volatility specifications, which may allow for instantaneous spatial effects. In the context of these models, $Y_t(\boldsymbol{s})$ is now repeatedly observed for $t = 1, \ldots, T$ and at all locations $\boldsymbol{s} \in D$. In the case of spatial econometrics models, the set of locations is assumed to be constant over time. This means that we can observe the outcome variable's realisation in the same space across multiple time periods. It is important to note that spatiotemporal models are naturally included in the purely spatial models, as time could be considered as one dimension of the points $\boldsymbol{s}$. In such cases, the weight matrix $\mathbf{W}$ has to be chosen accordingly so that future values do not influence past observations, as we will point out in Section (ref).

Compared to multivariate time series described in Sections (ref) and (ref), spatiotemporal models account for spatial, temporal, and spatiotemporal dependence. Typically, a certain (geographical) structure of the effects is assumed to interpret the parameters in a geographical sense. Moreover, this implied structure makes the models suitable for cases when $n$ is larger than $T$.

As an example of a spatiotemporal process, we can consider the log-returns of the condominium sales in Berlin across time. Figure (ref) depicts the spatiotemporal process on a map for one selected time point, June 2012 (left panels), and as time series for one selected location, postcode region 12683 (right panels). Moreover, the observed monthly log returns are shown in the top panels, and the squared log returns are shown in the bottom panels.

Another example of a spatiotemporal process from finance is the series of returns of all Dow Jones stocks across time. The similarity between the companies leads to interdependence between the series. For instance, fulle2022spatial showed that the similarity of firms regarding their balance sheet data can be used to model interactions in the log-volatilities across financial networks. To represent the closeness of the stocks, we have plotted the returns in an artificial space reflecting the distances in the volatility behaviour according to piccolo1990distance, so-called Piccolo distances. Figure (ref) shows the log returns and squared returns for one selected time point, 30 December 2022, in a spatial representation on a map with Voronoi cells separating the stocks.

Multi-index representation

In general, spatiotemporal models are already included in the spatial models when time is treated as one dimension of the locations $\boldsymbol{s}$, and the spatial weight matrices are chosen appropriately. All the above-mentioned models assume an arbitrary (finite) dimension of the underlying spatial domain. For instance, consider a spatiotemporal process on the surface of the Earth with degrees latitude and longitude, $(\text{lat}, \text{long})'$. Then, the index of the $i$-th observation, $\boldsymbol{s}_i$, would be composed of the spatial coordinate $(\text{lat}_i, \text{long}_i)'$ and the time point $t_i$ of the $i$-th observation, i.e., $\boldsymbol{s}_i = (t_i, \text{lat}_i, \text{long}_i)'$. The key difference between the time index and the spatial coordinates is that future observation cannot influence past observations. That is, the weight matrix $\mathbf{W}$ must account for the causal ordering in time, and $w_{ij}$ must be equal to zero if $t_i > t_j$. For example, let

be an $nT$ dimensional vector of a spatiotemporal random process at $n$ locations and $T$ time points. With two spatial weight matrices

where $\mathbf{L}_1$ is a $T$-dimensional shift matrix (i.e., first subdiagonal equalling one), $\otimes$ is the Kronecker product, and $\mathbf{W}$ is a regular $n$-dimensional spatial weights matrix (e.g., contiguity matrix), we obtain a spatiotemporal GARCH model with a first-order temporal lag implied by $\mathbf{W}_{\text{time}}$ and a first-order spatial lag with a constant weight matrix $\mathbf{W}$ directly from spatial GARCH models as proposed by otto2022general. Similarly, the log-GARCH models of Takaki:2021 can be constructed for spatiotemporal data. In such a way, spatiotemporal processes can be modelled using spatial volatility models, and all theoretical results can be directly applied in the spatiotemporal setup.

Notice that the above-defined matrices are sparse. Current computational algorithms for sparse matrices are highly efficient from a time and memory perspective, as they only store the indices and values of the non-zero entries in the weight matrices. Thus, spatiotemporal models can often be estimated in this multi-index representation. However, when not using sparse-element objects and operations, the computational requirements can quickly explode with an increasing dimension $n$ and $T$. Thus, it is generally meaningful if the dimension of the spatial weight matrices is as small as possible, and we will return back to the representation with an index $t$ below, i.e., $Y_t(\boldsymbol{s})$ is observed for $t = 1, \ldots, T$ at all locations $\boldsymbol{s} \in D$.

Spatiotemporal ARCH and GARCH models

The total number of coefficients in multivariate ARCH and GARCH models can increase faster than the cross-sectional dimension of these models. For example, in the BEKK specification considered by engle1995multivariate, the number of parameters has an order of $O(n^2)$, where $n$ is the cross-sectional dimension of the model. Hence, multivariate models are often not applicable in realistic spatiotemporal settings. caporin2015proximity consider parsimonious structured specifications using spatial econometrics tools. Consider the following BEKK specification:

where $\boldsymbol{Y}_t=(Y_t(\boldsymbol{s}_1,\hdots,\boldsymbol{Y}_t(\boldsymbol{s}_n))^{'}$ is the $n\times1$ vector of the outcome variable, $\mathbf{\Sigma}_t$ is the $n\times n$ matrix of covariances, $\boldsymbol{\operatorname{\varepsilon}}_t$ is the $n\times1$ vector of i.i.d. random variables that have $0$ mean and unit variance, and $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ are $n\times n$ matrices of unrestricted parameters. caporin2015proximity re-parametrise this model such that the number of parameters has an order of $O(n)$ by setting:

for $j=0,1$, where $\mathbf{W}$ is an $n\times n$ spatial weights matrix, and $\boldsymbol{s}^{(1)}$, $\boldsymbol{v}$, $\boldsymbol{\alpha}^{(j)}$ and $\boldsymbol{\beta}^{(j)}$ are $n\times1$ vectors of unknown parameters. Note that a homogeneous parameters version can be obtained by assuming that the parameter vectors do not vary over the cross-sectional dimension. In this specification, the spillover effects arising through $\mathbf{A}\boldsymbol{Y}_{t-1}\boldsymbol{Y}^{'}_{t-1}\mathbf{A}^{'}$ and $\mathbf{B}\mathbf{\Sigma}_{t-1}\mathbf{B}^{'}$ depend on the specification adopted for $\mathbf{W}$. caporin2015proximity discuss alternative specifications for $\mathbf{W}$ and consider a quasi-likelihood estimation approach for the model. Monica:2021 suggest an extended version of this model by assuming that $\mathbf{W}$ has a time-varying structure. However, notice that these models do not allow for instantaneous spatial spillovers in a GARCH sense. All spatial interactions enter the model at the first temporal lag, i.e., it must always take one time period for information to spill over to neighbouring locations.

otto2022dynamic consider a spatiotemporal ARCH model with a logarithmic representation for the volatility equation. This model takes the following form:

for $i=1,2,\hdots,n$ and $t=1,\hdots T$. Here, $\log h_{t}(\boldsymbol{s}_i)$ is considered as the log volatility term in the region $\boldsymbol{s}_i$ at time $t$, and $\operatorname{\varepsilon}_{t}(\boldsymbol{s}_i)$'s are i.i.d. random variables that have mean zero and unit variance. In (ref), $\{m_{l,ij}\}_{l=1}^p$, for $i,j=1,\hdots,n$, are the non-stochastic spatial weights, where $p$ is a finite positive integer, and $\{m_{l,ii}\}_{l=1}^p$ are zero for $i=1,\hdots,n$. The spatial, temporal, and spatiotemporal effects of the log-squared outcome variable on the log-volatility term are measured by the unknown parameters $\gamma_0$, $\{\rho_{l0}\}_{l=1}^p$, and $\{\delta_{l0}\}_{l=1}^p$, respectively. In (ref), $\mathbf{x}_{t}(\boldsymbol{s}_i)$ is a $k\times1$ vector of exogenous variables with the associated parameter vector $\boldsymbol{\beta}_0$, and the regional and time fixed effects are denoted by $\boldsymbol{\mu}_0=(\mu_{0}(\boldsymbol{s}_1),\hdots,\mu_{0}(\boldsymbol{s}_n))^{'}$ and $\boldsymbol{\alpha}_0=(\alpha_{10},\hdots,\alpha_{T0})^{'}$. Both $\boldsymbol{\mu}_0$ and $\boldsymbol{\alpha}_0$ can be correlated with the exogenous variables in an arbitrary manner.

To motivate the presence of spatial, temporal, and spatiotemporal effects in the log-volatility equation, otto2022dynamic consider a monthly dataset of the real house price returns in Berlin at the postcode level over the period from January 1995 to December 2015. Figure (ref) displays the average log-squared returns over Berlin's postcodes (the top left figure), the estimated temporal autocorrelation of the log-squared returns as a series of boxplots (the top right figure), and the estimated spatiotemporal autocorrelation in terms of Moran's $I$ across the time horizon (the bottom figure). The first figure shows a clustering pattern in the log-squared returns, indicating the presence of spatial dependence. From the ACF estimates, we can observe apparent temporal volatility clustering, while the spatiotemporal dependence is of a minor degree, irregularly fluctuating around zero.

Applying the log-squared transformation to the outcome equation in (ref) yields

where $Y^*_{t}(\boldsymbol{s}_i)=\log Y^2_{t}(\boldsymbol{s}_i)$, $h^{*}_{t}(\boldsymbol{s}_i)=\log h_{t}(\boldsymbol{s}_i)$ and $\operatorname{\varepsilon}^*_{t}(\boldsymbol{s}_i)=\log\operatorname{\varepsilon}^2_{t}(\boldsymbol{s}_i)$. Then, the model can be expressed in vector form as

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

for $t=1,\hdots,T$. The elements of $\boldsymbol{\operatorname{\varepsilon}}^{*}_t$ in (ref) are i.i.d across $i$ and $t$ but their mean may not be zero. otto2022dynamic suggest adding and subtracting $\operatorname*{E}\left(\boldsymbol{\operatorname{\varepsilon}}^{*}_t\right)$ to obtain the following estimation equation:

where $\boldsymbol{u}_t=(u_{t}(\boldsymbol{s}_1),\hdots,u_{t}(\boldsymbol{s}_n))^{'}=\boldsymbol{\operatorname{\varepsilon}}^{*}_t-\operatorname*{E}\left(\boldsymbol{\operatorname{\varepsilon}}^{*}_t\right)$, and $\mu_{\operatorname{\varepsilon}}=\operatorname*{E}\left(\operatorname{\varepsilon}^{*}_{t}(\boldsymbol{s}_i)\right)$. Let $\sigma^2_0=\operatorname*{E}(u^2_{t}(\boldsymbol{s}_i))$. Then, it follows that the elements of $\boldsymbol{u}_t$ are i.i.d across $i$ and $t$ with mean zero and variance $\sigma^2_0$. To eliminate $\boldsymbol{\mu}_0$ and $\mu_{\operatorname{\varepsilon}}\mathbf{1}_n$ from the model, they consider an orthonormal transformations. otto2022dynamic suggest a GMM approach based on a set of linear and quadratic moment functions Yu:2014. The Monte Carlo simulation results indicate that the suggested GMM estimator has good finite sample properties. For details, we refer the interested reader to otto2022dynamic directly.

When considering the same temporal ARCH parameter across spatial locations, one can see an additional benefit of spatiotemporal ARCH and GARCH models. Whereas their time-series counterparts often require long time series (without structural breaks) to reliably estimate the parameters, spatiotemporal models can also be estimated for shorter time frames because they can borrow information from all locations to estimate the temporal ARCH parameter $\gamma_0$. Simultaneously, spatial interactions between the series are included to account for the cross-sectional dependence.

We conclude this section by considering the circular spatiotemporal GARCH model considered by holleland2020stationary. This model is a circular version of the following spatiotemporal GARCH model:

where $\{V_t(\boldsymbol{s})\}$ is a sequence of i.i.d. random variables that have zero mean and unit variance, $\Delta_{1i}=\{\mathbf{v}\in\mathbb{Z}^d:\alpha_i(\mathbf{v})>0\}$, $\Delta_{2i}=\{\mathbf{v}\in\mathbb{Z}^d:\beta_i(\mathbf{v})>0\}$, $\mu$ is an unknown scalar parameter, $\boldsymbol{\alpha}=\{\alpha_i(\mathbf{v}), \mathbf{v}\in\Delta_{1i}, i=1,\hdots,p\}$ and $\boldsymbol{\beta}=\{\beta_i(\mathbf{v}), \mathbf{v}\in\Delta_{1i}, i=1,\hdots,p\}$ are unknown parameters. Let $\mathcal{R}=\mathbb{Z}^d/(\mathbf{m}\mathbb{Z}^d)$ be the quotient group of order $\mathbf{m}\in\mathbb{Z}^d_{+}$. The circular model considered in holleland2020stationary is obtained by replacing $\mathbb{Z}^d$ with $\mathcal{R}$. In this way, $\{Y_t(\boldsymbol{s})\}$ becomes a process indexed on $\mathbb{Z}\times\mathcal{R}$, and the difference $(\mathbf{v}-\boldsymbol{s})$ and the sum $(\mathbf{v}+\boldsymbol{s})$ are points in $\mathcal{R}$ with respect to modulus $\mathbf{m}$. Thus, the circular model takes the following form:

holleland2020stationary present the statistical properties of this circular model and suggest a quasi-maximum likelihood estimator for estimating the parameters.

Multivariate spatiotemporal ARCH models

When we observe vectors of several features across space and time, we can model such data as a multivariate spatiotemporal process. otto2022multivariate extended the spatiotemporal log-GARCH models for multivariate response variables. For that reason, they follow the idea of VEC-GARCH models, where the matrix of volatilities is transformed using the vech-operator bollerslev1988capital,engle1995multivariate. More precisely, the multivariate spatiotemporal log-ARCH model is given by

where $\boldsymbol{Y}_t$ is an $n \times r$-dimensional matrix of all $r$ features of the response variable at all $n$ locations. Moreover, $\mathbf{\Xi}_t$ is an $n \times r$-dimensional matrix of i.i.d random errors with mean zero and unit variance, and $\mathbf{H}_t$ is the matrix of log-volatilities, which are defined as follows

The operations $(\log)$ and $(\log, 2)$ should be understood as element-wise operations, i.e., $\mathbf{H}^{(\log)}_t = \left(\log h_{j,t}(\boldsymbol{s}_i)\right)_{i = 1, \ldots, n, j = 1, \ldots, r}$, and $\boldsymbol{Y}_t^{(\log,2)} = \left(\log Y_{j,t}^2(\boldsymbol{s}_i)\right)_{i = 1, \ldots, n, j = 1, \ldots, r}$ denotes the $n \times r$-dimensional matrix of log-squared observations. The model has three coefficient matrices: (i) the location-specific intercepts $\mathbf{A}$ ($n \times r$-dimensional), (ii) own and cross-variable spatial effects in $\mathbf{\Psi}$ ($r \times r$-dimensional), and (iii) own and cross-variable temporal effects in $\mathbf{\Pi}$ ($r \times r$-dimensional). Using the log-square transformation, the model can be transformed into a multivariate spatiotemporal autoregressive model. otto2022multivariate propose a QML estimator based on the normal approximation of the transformed errors. The QML estimator is based on the estimation approach proposed by yang2017identification for a multivariate spatial autoregressive model and Yu08 for the univariate spatiotemporal autoregressive model. Alternatively, the Gaussian mixture approximation could also be considered.

Spatiotemporal stochastic volatility models

We can utilise two approaches to define a spatiotemporal stochastic volatility model. The first approach, known as the geostatistical approach Cressie:2015, involves using known parametric covariance functions to model the correlation over time and space. This approach usually includes models obtained by extending a stationary and isotropic Gaussian process (GP), which can be defined in the following way:

where $\sigma>0$ is a scalar unknown parameter and $V_t(\boldsymbol{s})$ is a Gaussian process with mean $0$, variance $1$, and stationary covariance function $\text{Cov}\left(V_{t_1}(\boldsymbol{s}_1),V_{t_2}(\boldsymbol{s}_2)\right)=C(\mathbf{d},\lambda)$, where $\mathbf{d}=\boldsymbol{s}_1-\boldsymbol{s}_2$ and $\lambda=t_1-t_2$. In this approach, $C(\mathbf{d},\lambda)$ takes a known parametric function that satisfies certain properties such as stationarity, separability, and full symmetry. Gneiting:2007 review main covariance functions suggested in the literature.

Peter:2011 consider an extension of (ref) for temperature series, which takes the form of $Y_t(\boldsymbol{s})=\sigma_t(\boldsymbol{s}) V_t(\boldsymbol{s})$, where $\sigma_t(\boldsymbol{s})$ is modelled such that it can capture the spatially varying seasonality in the variance of $Y_t(\boldsymbol{s})$. huang2011class suggest another extended version taking the following form:

where $\sigma>0$ and $\alpha>0$ are scalar unknown parameters, $h_t(\boldsymbol{s})$ is a Gaussian process with mean $0$, variance $1$, and correlation function $\rho_{h}$, and $V_t(\boldsymbol{s})$ is another independent Gaussian process with mean $0$, variance $1$ and covariance function $\rho_V$. huang2011class assume that the correlation functions $\rho_h$ and $\rho_V$ are isotropic and stationary in time in the following sense:

This model can be considered as a spatiotemporal extension of the non-Gaussian geostatistical model proposed by Steel:2006. For the estimation of the model, huang2011class consider an approximation to the likelihood function of the model and show how an importance sampling approach and Monte Carlo integration can be used to evaluate and maximise the approximation.

In the second approach, tools from spatial econometrics are used to specify spatial, temporal, and spatiotemporal effects in the log-volatility equation. These effects are incorporated into the model specification through spatial weights matrices that specify the degree of spatial dependence in the outcome variables. Otto:2022 suggested the following specification:

where $\boldsymbol{Y}_t=\left(Y_t(\boldsymbol{s}_1),\hdots,Y_t(\boldsymbol{s}_n)\right)^{'}$ is the $n\times1$ vector of outcome variable at time $t$, $\mathbf{H}^{1/2}_t=\text{diag}(e^{\frac{1}{2}h_{t}(\boldsymbol{s}_1)},\hdots,e^{\frac{1}{2}h_{t}(\boldsymbol{s}_n)})$ is the $n\times n$ diagonal matrix containing stochastic volatility terms $h_{t}(\boldsymbol{s}_i)$'s, and $\mathbf{V}_{t}=\left(V_{t}(\boldsymbol{s}_1),\hdots,V_{t}(\boldsymbol{s}_n)\right)^{'}$ is the $n\times 1$ vector of disturbance terms whose elements are i.i.d standard normal random variables. Let $\boldsymbol{h}_t=\left(h_{t}(\boldsymbol{s}_1),\hdots,h_{t}(\boldsymbol{s}_n)\right)^{'}$ be the $n\times1$ vector of stochastic volatility terms at time $t$. Otto:2022 consider the following process for $\boldsymbol{h}_t$:

where $\boldsymbol{\mu}=\left(\mu(\boldsymbol{s}_1),\hdots,\mu(\boldsymbol{s}_n)\right)^{'}$ is the $n\times1$ vector of the time-invariant site-specific effects, and $\mathbf{U}_{t}=\left(u_{t}(\boldsymbol{s}_1),\hdots,u_{t}(\boldsymbol{s}_n)\right)^{'}$ is the $n\times1$ vector of i.i.d. disturbance terms with $u_{t}(\boldsymbol{s}_i)\sim N(0, \sigma^2)$ for all $i$ and $t$, where $\sigma^2$ is a scalar unknown parameter. The $n\times n$ spatial weights matrix $\mathbf{W}$ specifies the degree of linkages among the elements of $\boldsymbol{h}_t$. The scalar parameter $\rho_1$ captures contemporaneous spatial correlation, $\rho_2$ measures the temporal effect, i.e., the time dynamic effect, and $\rho_3$ represents the spatiotemporal effect, i.e., the spatial diffusion effect. The reduced form of the volatility equation can be expressed as

where $\mathbf{S}(\rho_1)=(\mathbf{I}_n-\rho_1\mathbf{W})$ and $\mathbf{A}(\rho_2,\rho_3)=(\rho_2\mathbf{I}_n+\rho_3\mathbf{W})$. When the cross-sectional dimension is fixed, the process for the log-volatility is stable if all eigenvalues of $\mathbf{S}^{-1}(\rho_1)\mathbf{A}(\rho_2,\rho_3)$ lie inside the unit ball. Otto:2022 show that an easy-to-check sufficient condition for the stability of the volatility equation is $|\rho_1|+|\rho_2|+|\rho_3|<1$ when $\mathbf{W}$ is row normalised. The statistical properties of this model are described in Section of Appendix.

Otto:2022 introduce a Bayesian estimation approach by assuming the following independent prior distributions: $\rho_1\sim\text{Uniform}(-1,1)$, $\rho_2\sim\text{Uniform}(-1,1)$, $\rho_3\sim\text{Uniform}(-1,1)$, $\boldsymbol{\mu}|\mathbf{b}_{\mu},\mathbf{B}_{\mu}\sim N(\mathbf{b}_{\mu},\mathbf{B}_{\mu})$, and $\sigma^2|a,b\sim\text{IG}(a,b)$. The estimation approach is based on the log-squared transformation and the finite Gaussian mixture approach described in Section (ref) for the SSV model. The log-squared transformation gives

where $Y^*_{t}(\boldsymbol{s}_i)=\log Y^2_{t}(\boldsymbol{s}_i)$ and $V^*_{t}(\boldsymbol{s}_i)=\log V^2_{t}(\boldsymbol{s}_i)$. Define $\boldsymbol{Y}^*_t=\left(Y^*_{t}(\boldsymbol{s}_1),Y^*_{t}(\boldsymbol{s}_2),\hdots,Y^*_{t}(\boldsymbol{s}_n)\right)^{'}$ and $\mathbf{V}^*_t=\left(V^*_{t}(\boldsymbol{s}_1),V^*_{t}(\boldsymbol{s}_2),\hdots,V^*_{t}(\boldsymbol{s}_n)\right)^{'}$. Then, in vector form, we have

In order to convert (ref) into a linear Gaussian state-space model, Otto:2022 approximate $p(V^*_{t}(\boldsymbol{s}_i))$ with an $m$-component Gaussian mixture distribution:

where $\phi(V^*_{t}(\boldsymbol{s}_i)|\mu_j,\,\sigma^2_j)$ denotes the Gaussian density function with mean $\mu_j$ and variance $\sigma^2_j$, $p_j$ is the probability of $j$th mixture component and $m$ is the number of components. We can equivalently write (ref) in terms of an auxiliary discrete random variable $z_{t}(\boldsymbol{s}_i)\in\{1,2,\hdots,m\}$ that serves as the mixture component indicator:

where $\mathbb{P}(z_{t}(\boldsymbol{s}_i)=j)=p_j$ is the probability that $z_{t}(\boldsymbol{s}_i)$ takes the $j$th value. Let $\mathbf{Z}_t=(z_{t}(\boldsymbol{s}_1),\hdots,z_{t}(\boldsymbol{s}_n))^{'}$, $\mathbf{d}_t=(\mu_{z_{t}(\boldsymbol{s}_1)},\hdots,\mu_{z_{t}(\boldsymbol{s}_n)})^{'}$ and $\boldsymbol{\Sigma}_t=\operatorname*{diag}(\sigma^2_{z_{t}(\boldsymbol{s}_1)},\hdots,\sigma^2_{z_{t}(\boldsymbol{s}_n)})$. Then, from (ref), we have $\mathbf{V}^*_t|\mathbf{Z}_t\sim N(\mathbf{d}_t,\,\boldsymbol{\Sigma}_t)$, which indicates that

Otto:2022 suggest the Gibbs sampler described in Algorithm (ref) in the Appendix for estimating the model parameters.

Further extensions and related models

In this section, we consider alternative specifications for the spatial and spatiotemporal effects in the volatility equations. In Section (ref), we consider some extensions of the basic SSV model introduced in Section (ref). In both spatial and spatiotemporal volatility models covered in Sections (ref) and (ref), we usually consider spatial lag terms formulated with $\mathbf{W}$ to introduce spatial and spatiotemporal effects in the volatility equations. In Sections (ref) and (ref), we consider alternative methods based on the matrix exponential and conditional autoregressive approaches to specify spatial dependence in the volatility equations.

Extensions of spatial stochastic volatility models

Following the time series literature on the stochastic volatility models described in Section (ref), the SSV model described above can be extended in several directions. To this end, we consider several different spatial stochastic volatility models depending on the spatial specification adopted for the outcome and the log-volatility equations. Let $\boldsymbol{x}(\boldsymbol{s}_i)=(x_1(\boldsymbol{s}_i),\hdots,x_k(\boldsymbol{s}_i))^{'}$ be the $k\times1$ vector of explanatory variables for $i=1,2,\hdots,n$, $\mathbf{X}=(\boldsymbol{x}(\boldsymbol{s}_1),\hdots, \boldsymbol{x}(\boldsymbol{s}_n))^{'}$ be the $n\times k$ matrix of explanatory variables, and $\mathbf{W} = (w_{ij})_{i,j = 1, \ldots, n}$ and $\mathbf{M} = (m_{ij})_{i,j = 1, \ldots, n}$ be two non-stochastic $n\times n$ spatial weights matrices that have zero diagonal elements.

The first extension we consider allows for a spatial lag of the outcome variable as well as some exogenous explanatory variables in the outcome equation:

for $i=1,\hdots,n$. The first-order spatial autoregressive process for $Y(\boldsymbol{s}_i)$'s introduces spatial correlations in the outcome variable. As before, $\operatorname{\varepsilon}(\boldsymbol{s}_i)$'s are i.i.d. standard normal random variables, and $u(\boldsymbol{s}_i)$'s are an i.i.d. normal random variables with mean zero and variance $\sigma^2_u$. The scalar spatial autoregressive parameters are $\rho$ and $\phi$, and $\mu_h$ is the constant mean parameter. We will refer to this model as SAR-SSV.

The spatial autoregressive process allows for the global transmission of a shock. On the other hand, the spatial moving process transmits a shock locally anselin1988spatial, Fingleton:2008,Fingleton:2008b, Taspinar:2013, Dogan:2015. An alternative specification to the SSV model is to use a spatial moving average process for the log-volatility terms:

where $\phi$ is the scalar spatial moving average parameter and $u(\boldsymbol{s}_i)$'s are an i.i.d normal random variables with mean zero and variance $\sigma^2_u$. We refer to this model as the SMA-SSV model.

One can also allow for both a spatial autoregressive process and a spatial moving average process in the log-volatility equation to define the following extension of the SSV model:

where $\mathbf{W}_1 = (w_{1,ij})_{i,j = 1, \ldots, n}$ and $\mathbf{W}_2 = (w_{2,ij})_{i,j = 1, \ldots, n}$ are two non-stochastic $n\times n$ spatial weights matrices that have zero diagonal elements. We refer to this model as the SARMA-SSV model.

Another variant can be defined by allowing the volatility feedback in the outcome (observation) equation, which can be considered as an analogous version of the model suggested by Koopman:2002. This model is specified as

where the scalar parameter $\alpha$ indicates the effect of the stochastic volatility on $Y(\boldsymbol{s}_i)$. We refer to this model as the SARM-SSV model.

An analogous version of the leverage effect model described in Section (ref) can be defined by introducing correlation in the error terms of the outcome and log-volatility equations. This analogous takes the following form:

where $(\operatorname{\varepsilon}(\boldsymbol{s}_i),u(\boldsymbol{s}_i))^{'}$ follows the following bivariate normal distribution

and $\varrho$ is the unknown correlation parameter. We refer to this model as the SAR-SSVL model.

The next extension considers a scale mixture of normal distribution representation for the outcome variable. This is the analogous version described in Section (ref) and takes the following form:

where the latent scale variables $\omega(\boldsymbol{s}_i)$'s are i.i.d random variables having distribution $IG(\nu/2,\nu/2)$. It is well-known that this representation implies that the marginal distribution of $\omega(\boldsymbol{s}_i)^{1/2}\operatorname{\varepsilon}(\boldsymbol{s}_i)$ (unconditional on $\omega(\boldsymbol{s}_i)$) is the $t$ distribution with $\nu$ degrees of freedom Geweke:1993. This specification can be considered as the spatial extension of the stochastic models considered in Harvey:1994, Ruiz:1994 and Eric:2004. We refer to this model as the SAR-SSVt model.

The Bayesian estimation approach described in Algorithm (ref) can also be considered for some of these alternative specifications. The main requirement of the estimation approach is that the model obtained through the log-squared transformation and the Gaussian mixture approximation should be in the form of a linear Gaussian state space model. Therefore, a similar approach can be adopted to estimate the SAR-SV, SMA-SSV, SARMA-SSV, SAR-SSVL, and SAR-SSVt models. However, the approach described for estimating the SSV model may not be extended to the SARM-SSV model.

Matrix-exponential dependence structure

An alternative way to define spatial lag terms is through a matrix exponential term defined as $e^{\alpha \mathbf{W}}=\sum_{j=0}^\infty\frac{\alpha^j}{j!}\mathbf{W}^j$, where $\alpha$ is a scalar spatial parameter. The matrix exponential terms were first considered by Lesage:2007 to specify spatial dependence in an outcome variable as an alternative to the commonly used spatial autoregressive process. This specification type introduces an exponential decay rate for the cross-sectional dependence and can provide some computational advantages. su2023statistical consider a matrix exponential term for spatial log-ARCH models, which is a special case of their logarithmic spatial heteroscedasticity model (log-SHE model). For further information, refer to Debarsy:2015 and Ye:2021, Ye:2022, Ye:2023 for the properties of spatial models defined in terms of matrix exponential terms.

Using matrix exponential terms, we can alternatively define the spatial ARCH and GARCH processes mentioned in Section (ref), respectively, in the following way:

where $\alpha_1$ and $\beta_1$ are scalar parameters. When $\alpha_1=\beta_1=0$, both processes reduce to $\boldsymbol{h}=\alpha_0\mathbf{1}_n$. Since a matrix exponential term is always invertible, the spatial GARCH process in (ref) always has the reduced form given by

where $e^{-\beta_1\mathbf{W}_2}$ is the inverse of $e^{\beta_1\mathbf{W}_2}$. Importantly, this reduced form does not invoke any restriction on the parameter space of $\beta_1$, which was not the case for the model considered in Section (ref). Similarly, we can also consider the matrix exponential terms to define alternative versions of the spatial stochastic volatility models introduced in Section 3.2. For example, the log-volatility equation in the SSV model can be specified in the following way:

where $\phi$ is the scalar spatial parameter. The reduced form of this process is $\boldsymbol{h}= \mu_he^{-\phi\mathbf{W}}\mathbf{l}_n+ e^{-\phi\mathbf{W}}\boldsymbol{u}$, which does not require any restrictions for $\phi$. Similarly, we can also introduce spatial and spatiotemporal effects in the models considered in Sections (ref) and (ref). For example, the matrix exponential version of the spatiotemporal ARCH model considered by otto2022dynamic can be formulated as

where $e^{\sum_{l=1}^p\rho_{l0}\mathbf{M}_l}$ and $e^{\sum_{l=1}^p\delta_{l0}\mathbf{M}_l}$ are higher-order terms formulated by the sequence of weights matrices $\{\mathbf{M}_l\}$. Similarly, the log-volatility equation in the spatiotemporal model suggested by Otto:2022 can alternatively be expressed as

where $\rho_1$ and $\rho_3$ are spatial and spatiotemporal parameters.

Conditional autoregressive dependence structure

Another alternative way that can be used to introduce spatial and spatiotemporal dependence in a volatility process is to use the conditional autoregressive (CAR) specification. We start with the following model considered by Besag:1991:

where $\mu$ is a scalar unknown overall mean parameter, $\operatorname{\varepsilon}(\boldsymbol{s}_i)$ is a disturbance term that has $N(0,\sigma^2_{\operatorname{\varepsilon}})$ distribution with the unknown variance parameter $\sigma^2_{\operatorname{\varepsilon}}$, $\phi(\boldsymbol{s}_i)$ is a spatially structured random effect term that has the CAR specification in (ref). In the CAR process, we assume that $b_{ij}$'s are known constants with $b_{ij}=b_{ji}$ and $b_{ii}=0$, and $\sigma^2_{\phi}$ is an unknown scalar variance parameter. The constants $b_{ij}$ can be considered as spatial weights that determine the relationship between regions $\boldsymbol{s}_i$ and $\boldsymbol{s}j$. For example, $b_{ij}$ may be set to 1 if $\boldsymbol{s}_i$ and $\boldsymbol{s}_j$ are neighbours, and $0$ otherwise. In this specification, the conditional variance of $\phi(\boldsymbol{s}_i)$ is spatially varying and depends on $\sum_{j\ne i}b_{ij}$. Besag:1974 shows that the joint distribution of $\boldsymbol{\phi}=(\phi(\boldsymbol{s}_1),\hdots,\phi(\boldsymbol{s}_n))^{'}$ can be determined as

where $\mathbf{B}$ is the $n\times n$ matrix with the $(i,j)$th element $B_{ij}=\frac{b_{ij}}{\sum_{k=1}^nb_{ik}}$ and $\mathbf{D}_{\phi}$ is the $n\times n$ diagonal matrix with the $i$th diagonal element $D_{ii}=\frac{1}{\sum_{k=1 }^nb_{ik}}$. Note that this result requires that $(\mathbf{I}_n-\mathbf{B})^{-1}\mathbf{D}_{\phi}$ is a positive definite matrix. Yan:2007 extends this model by assuming that $\operatorname{\varepsilon}(\boldsymbol{s}_i)$ has a stochastic volatility term as specified below:

where $\mu_h$ is a scalar mean parameter and $h(\boldsymbol{s}_i)$ is the log-volatility term assumed to follow the CAR process specified in (ref). In the CAR process, the weights $c_{ij}$ play the same role as $b_{ij}$ in the CAR process assumed for $\phi(\boldsymbol{s}_i)$, and $\sigma^2_h$ is a scalar variance parameter. As in the case of $\boldsymbol{\phi}$, the joint distribution of $\boldsymbol{h}$ is $\boldsymbol{h}|\sigma^2_h\sim N\left(\mathbf{0},\sigma^2_h(\mathbf{I}_n-\mathbf{C})^{-1}\mathbf{D}_{h}\right)$, where we assume that $(\mathbf{I}_n-\mathbf{C})^{-1}\mathbf{D}_{h}$ is a positive definite matrix. In order to achieve identification for $\mu$ and $\mu_h$, Yan:2007 respectively requires that $\sum_{i=1}^n\phi(\boldsymbol{s}_i) = 0$ and $\sum_{i=1}^nh(\boldsymbol{s}_i) = 0$. The posterior distribution of the model can be expressed as

where $\boldsymbol{\theta}=(\mu,\boldsymbol{\phi},\mu_h,\sigma^2_h)^{'}$ and $f(\boldsymbol{Y}|\mu,\boldsymbol{\phi},\mu_h,\boldsymbol{h})$ is the likelihood function given by

Yan:2007 assumes flat priors for the elements of $\boldsymbol{\theta}$ and shows that the conditional posterior distributions take standard forms, except for that of $\boldsymbol{h}$. Yan:2007 defines $\lambda(\boldsymbol{s}_i)=\mu_h+h(\boldsymbol{s}_i)$ and shows that the conditional posterior distribution of $\lambda(\boldsymbol{s}_i)$ can be bounded, up to a scale, by a density of normal distribution. Yan:2007 suggests using this blanket distribution in an accept-reject algorithm to generate draws for $\lambda(\boldsymbol{s}_i)$. Note that once we have draws for $\lambda(\boldsymbol{s}_i)$, we can determined $\mu_h$ and $h(\boldsymbol{s}_i)$ from the relation $\lambda_i(\boldsymbol{s}_i)=\mu_h+h(\boldsymbol{s}_i)$ via imposing the identification condition $\sum_{i=1}^nh(\boldsymbol{s}_i) = 0$.

Lesage:2008 suggested a version of the model in (ref) and (ref) that involves an alternative CAR process for the random effect term $\phi(\boldsymbol{s}_i)$. In vector form, their version can be specified as

where $\boldsymbol{\operatorname{\varepsilon}}=(\operatorname{\varepsilon}(\boldsymbol{s}_1),\hdots,\operatorname{\varepsilon}(\boldsymbol{s}_n))^{'}$ is the vector of disturbance terms, $\boldsymbol{\phi}=(\phi(\boldsymbol{s}_1),\hdots,\phi(\boldsymbol{s}_n))^{'}$ is the vector of random effects, $\sigma^2_{\phi}$ is a scalar variance parameter and $\rho$ is a scalar spatial parameter. Lesage:2008 consider this model for the knowledge spillovers arising from patent activity between European regions and specify the elements of the diagonal matrix $\mathbf{M}$ as the output gap and the elements of the symmetric matrix $\mathbf{W}$ as either based on a technological proximity index or based on an index of transport infrastructure. To ensure that $(\mathbf{I}_n-\rho\mathbf{W})^{-1}\mathbf{M}$ is positive definite, we should require that $\rho\in (1/\psi_{min},\,1/\psi_{max})$, where $\psi_{min}$ and $\psi_{max}$ are the minimum and maximum eigenvalues of $\mathbf{M}^{-1/2}\mathbf{W}\mathbf{M}^{1/2}$, respectively. In order to allow for outliers, Lesage:2008 assume that the elements of $\boldsymbol{\operatorname{\varepsilon}}$ have a scale mixture of normal distribution, i.e., $\operatorname{\varepsilon}(\boldsymbol{s}_i)|\omega(\boldsymbol{s}_i),\sigma^2_{\operatorname{\varepsilon}}\sim N\left(0,\sigma^2_{\operatorname{\varepsilon}}\omega(\boldsymbol{s}_i)\right)$. The scale mixture components $\mathbf{\omega}(\boldsymbol{s}_i)$'s are independent with $\omega(\boldsymbol{s}_i)|\nu\sim IG(\nu/2,\nu/2)$ and $\nu\sim Exp(\lambda_0)$, where $Exp(\lambda_0)$ is the exponential distribution with the rate parameter $\lambda_0$. For the remaining parameters, Lesage:2008 assume the following independent prior distributions: (i) $\sigma^2_{\phi}\sim IG(a_{\phi},b_{\phi})$, (ii) $\sigma^2_{\operatorname{\varepsilon}}\sim IG(a_{\operatorname{\varepsilon}},b_{\operatorname{\varepsilon}})$, (iii) $\mu\sim N(\mu_0, V_0)$ and (iv) $\rho\sim Beta(a_0,a_0)$, where $Beta(a_0,a_0)$ is the beta distribution defined as

Lesage:2008 show that setting $a_0=1.01$ gives a relatively uninformative prior for $\rho$ over the interval $(1/\psi_{min},\,1/\psi_{max})$. Under these priors, the conditional posterior distributions of $\boldsymbol{\phi}$, $\boldsymbol{\omega}$, $\mu$, $\sigma^2_{\phi}$ and $\sigma^2_{\operatorname{\varepsilon}}$ are in standard forms while those of $\rho$ and $\nu$ take unknown forms. In the case of $\rho$, Lesage:2008 use the univariate numerical integration over the interval $(1/\psi_{min},\,1/\psi_{max})$ to produce the conditional posterior distribution. As for $\nu$, Lesage:2008 utilised the random walk Metropolis-Hastings algorithm described in Lesage:2009 to generate posterior draws.

Conclusion and outlook

Spatial and spatiotemporal volatility models constitute a promising new class of models for modelling dependence among the volatility of neighbouring sites in spatial and spatiotemporal data. These models have been recently developed to account for the spatial and spatiotemporal effects in the volatility of an outcome variable. This paper has provided a comprehensive review of the recent literature on spatial and spatiotemporal volatility models. Compared to multivariate time-series GARCH models, which are typically not applicable in spatial settings because of the large number of cross-sectional locations, spatial and spatiotemporal volatility models incorporate a geographical structure to specify dependence in the volatilities. This structure allows modelling spatial and spatiotemporal spillovers across neighbouring locations in a GARCH-like sense. High volatilities may instantaneously spill over to neighbouring regions, thus forming spatial volatility clusters.

Besides motivating different alternative specifications and summarising estimation strategies, we discussed possible extensions and indicated future research directions. Notably, the strand of spatial and spatiotemporal volatility models can be extended in the following ways:

In summary, we believe that the class of spatial and spatiotemporal models offers intriguing opportunities for new theoretical developments and practical applications. Bearing in mind that a process' volatility is often interpreted as risk, identifying and predicting local areas of high volatilities -- risks -- is important in various fields, such as finance, economics, or environmental science.

Glossary