A Heterogeneous Spatiotemporal GARCH Model: A Predictive Framework for Volatility in Financial Networks

TL;DR

This paper presents a novel heterogeneous spatiotemporal GARCH model that captures spatial spillovers and contagion effects in financial volatility, allowing for predictions across unobserved locations using a flexible proximity measure.

Contribution

The paper introduces a new spatiotemporal GARCH framework that incorporates spatial dependence without relying on fixed adjacency matrices, using proxy-based proximity for financial networks.

Findings

The model effectively captures spatial spillovers in financial volatility.

It outperforms traditional models in predictive accuracy in cross-validation.

The framework enables volatility prediction at unobserved locations.

Abstract

We introduce a heterogeneous spatiotemporal GARCH model for geostatistical data or processes on networks, e.g., for modelling and predicting financial return volatility across firms in a latent spatial framework. The model combines classical GARCH(p, q) dynamics with spatially correlated innovations and spatially varying parameters, estimated using local likelihood methods. Spatial dependence is introduced through a geostatistical covariance structure on the innovation process, capturing contemporaneous cross-sectional correlation. This dependence propagates into the volatility dynamics via the recursive GARCH structure, allowing the model to reflect spatial spillovers and contagion effects in a parsimonious and interpretable way. In addition, this modelling framework allows for spatial volatility predictions at unobserved locations. In an empirical application, we demonstrate how the model can be applied to financial stock networks. Unlike other spatial GARCH models, our framework does not rely on a fixed adjacency matrix; instead, spatial proximity is defined in a proxy space constructed from balance sheet characteristics. Using daily log returns of 50 publicly listed firms over a one-year period, we evaluate the model's predictive performance in a cross-validation study.

Tables

Table 1. Table 1: Out-of-sample bias and RMSE of estimated nonparametric functions ω ​ ( s ) \omega(s) , α ​ ( s ) \alpha(s) , β ​ ( s ) \beta(s) , and predicted unconditional volatilities at S 2 S_{2} over M ​ C = 100 MC=100 simulations.

$n_1$	$T$	$\omega$		$\alpha$		$\beta$		Volatility
		Bias	RMSE	Bias	RMSE	Bias	RMSE	Bias	RMSE
	100	-0.0216	0.0784	0.0067	0.0774	0.0603	0.2126	0.0092	0.1205
50	200	-0.0308	0.0706	0.0042	0.0597	0.0827	0.1894	0.0042	0.1084
	300	-0.0348	0.0684	0.0047	0.0533	0.0922	0.1766	0.0037	0.1067
	100	-0.0177	0.0819	0.0073	0.0767	0.0447	0.2185	0.0020	0.1168
100	200	-0.0269	0.0718	0.0078	0.0605	0.0691	0.1850	0.0008	0.1096
	300	-0.0308	0.0651	0.0063	0.0527	0.0810	0.1671	0.0006	0.1058
	100	-0.0136	0.0819	0.0072	0.0820	0.0521	0.2182	0.0115	0.1231
150	200	-0.0240	0.0706	0.0056	0.0614	0.0698	0.1852	0.0002	0.1113
	300	-0.0270	0.0658	0.0058	0.0550	0.0771	0.1727	0.0004	0.1097

Table 2. Table 2: List of the 50 companies included in the analysis.

Companies
Allianz	JPM
Abbott	Merck
Apple	Microsoft
AT&T	MUFG
Santander	Nestlé
Bank of America	Novartis
BHP	Oracle
BNP Paribas	Pepsi
BP	Petrobras
Chevron	Pfizer
Cisco	PMI
Citigroup	Procter & Gamble
Coca-Cola	Roche
ConocoPhillips	Shell
E.ON	Samsung
Eni	Sanofi
ExxonMobil	Schlumberger
Gazprom	Siemens
General Electric	Telefónica
GlaxoSmithKline	Total
Google	Toyota
HP	Verizon
HSBC	Vodafone
IBM	Walmart
Intel	Johnson & Johnson

Equations

$$\{\begin{array}{l}{Z}_t(s)={\sigma }_t(s){\eta }_t(s)\\ {\sigma }_t^2(s)=\omega (s)+{\sum }_{i=1}^q{\alpha }_i(s)Z_{t-i}^2(s)+{\sum }_{j=1}^p{\beta }_j(s){\sigma }_{t-j}^2(s)\end{array}$$

$$\{\begin{array}{l}{\tilde{Z}}_t(s)={\tilde{\sigma }}_t(s){\eta }_t(s),\\ {\tilde{\sigma }}_t^2(s)=\omega (s_0)+{\sum }_{i=1}^q{\alpha }_i(s_0){\tilde{Z}}_{t-i}^2(s)+{\sum }_{j=1}^p{\beta }_j(s_0){\tilde{\sigma }}_{t-j}^2(s).\end{array}$$

$$\Omega =\{(\omega (s),{\alpha }_i(s),{\beta }_j(s)):\sum _{i=1}^q{\alpha }_i(s)+\sum _{j=1}^p{\beta }_j(s)\le 1-\gamma ,{\rho }_1\le \omega (s)\le {\rho }_2\text{ for }i=1,\dots ,q\text{ and }j=1,\dots ,p\},$$

$$Z_t^2(s)=\omega (s)+\sum _{i=1}^r{\delta }_i(s)Z_{t-i}^2(s)+\sum _{j=1}^p{\beta }_j(s){\zeta }_{t-j}(s)+{\zeta }_t(s),$$

$$Z_t^2(s)=\sum _{j=0}^{\infty }{\varphi }_j(s){\zeta }_{t-j}(s)+C(s)$$

$${\varphi }_j(s)={\beta }_j(s)+\sum _{i=1}^j{\delta }_i(s){\varphi }_{j-i}(s),j\ge 0$$

$$C(s)=\omega(s)\Big{(}1-\sum_{j=1}^{r}\delta_{j}(s)\Big{)}^{-1}.$$

$$Z_t^2(s)=\sum _{k=0}^{\infty }{\nu }_k(s){\zeta }_{t-k}(s)+C(s),$$

$$A(s)=\left(\begin{array}{ccccc}{\varphi }_1(s)& {\varphi }_2(s)& \cdots & {\varphi }_{p-1}(s)& {\varphi }_p(s)\\ 1& 0& \cdots & 0& 0\\ 0& 1& \ddots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \cdots & \cdots & 1& 0\end{array}\right)$$

$$\underset{s}{\sup }|{\nu }_k(s)|\le K{\lambda }^k.$$

$${\underline{Z}}_t^2\left(s_0\right)=A\left(s_0\right){\underline{Z}}_{t-1}^2\left(s_0\right)+{\underline{\zeta }}_t\left(s_0\right)+\underline{C}(s_0),$$

$${\underline{Z}}_t^2{\left(s_0\right)}^{\prime }=\left(Z_t^2\left(s_0\right),\dots ,Z_{t-p+1}^2\left(s_0\right)\right),{\zeta }_t{\left(s_0\right)}^{\prime }=\left({\zeta }_t\left(s_0\right),\dots ,0\right),\text{and}$$

$$\underline{C}(s_0)=\left(C(s_0),\dots ,0\right).$$

$$\mathbb{E}\left[Z_t^2\left(s_0\right)\mid \mathbf{Z}\right]=\sum _{k=0}^{t-p-1}{\nu }_k\left(s_0\right)\mathbb{E}\left[{\zeta }_{t-k}\left(s_0\right)\mid \mathbf{Z}\right]+{\left[A{\left(s_0\right)}^{t-p}\mathbb{E}\left\{{\underline{Z}}_p^2\left(s_0\right)\mid \mathbf{Z}\right\}\right]}_1+C(s_0)$$

$$\mathbb{E}\left[{\zeta }_{t-k}\left(s_0\right)\mid \mathbf{Z}\right]=\mathbb{E}\left[{\zeta }_{t-k}\left(s_0\right)\mid \left\{{\zeta }_{t-k}\left(s_u\right):u=1,\dots ,m\right\}\right]$$

$$\mathbb{E}\left[Z_t^2\left(s_0\right)\mid \mathbf{Z}\right]=$$

$$+{\left[A{\left(s_0\right)}^{t-p}\mathbb{E}\left\{{\underline{Z}}_p^2\left(s_0\right)\mid \mathbf{Z}\right\}\right]}_1+C(s_0).$$

$$\mathbb{E}\left[Z_t^2\left(s_0\right)\mid \mathbf{Z}\right]=\sum _{k=0}^{t-p-1}{\nu }_k\left(s_0\right)\mathbb{E}\left[{\zeta }_{t-k}\left(s_0\right)\mid \left\{{\zeta }_{t-k}\left(s_u\right):u=1,\dots ,m\right\}\right]+C(s_0)+O_p\left({\lambda }^{t-p}\right)$$

$$Z_t^2\left(s_0\mid m\right)$$

$${\underline{Z}}_t^2{\left(s_0\mid m\right)}^{\prime }$$

$${\underline{\zeta }}_t{\left(s_0\mid m\right)}^{\prime }$$

$$\sum _{k=0}^{t-p-1}{\nu }_k\left(s_0\right)\mathbb{E}\left[{\zeta }_{t-k}\left(s_0\right)\mid \left\{{\zeta }_{t-k}\left(s_u\right):u=1,\dots ,m\right\}\right]+C(s_0).$$

$$\mathbb{E}\left[{\zeta }_t\left(s_0\right)\mid {\left\{{\zeta }_t\left(s_u\right)\right\}}_{u=1}^m\right]=\sum _{u=1}^m{\gamma }_j(\underline{s}){\zeta }_t\left(s_u\right),$$

$$\mathbb{E}\left[Z_t^2\left(s_0\right)\mid \mathbf{Z}\right]=\sum _{k=0}^{t-p-1}{\nu }_k\left(s_0\right)\sum _{u=1}^m{\gamma }_u(\underline{s}){\zeta }_{(t-k)}\left(s_u\right)+C(s_0).$$

$$\underline{s}=\left(s_0,\dots ,s_m\right),\gamma (\underline{s})={\left({\gamma }_1(\underline{s}),\dots ,{\gamma }_m(\underline{s})\right)}^{\prime }=R(\underline{s}{)}^{-1}r(\underline{s})$$

$$R(\underline{s}{)}_{ij}=\mathrm{cov}\left({\zeta }_t\left(s_i\right),{\zeta }_t\left(s_j\right)\right)=c_{\zeta }\left(s_i-s_j\right)\text{ and }r(\underline{s}{)}_i=\mathrm{cov}\left({\zeta }_t\left(s_0\right),{\zeta }_t\left(s_i\right)\right)=c_{\zeta }\left(s_i-s_0\right)\text{.}$$

$$\text{Cov}({\hat{Z}}_t^2(s_0),{\hat{Z}}_{t^{\prime }}^2(s_0))=\text{Cov}\left(\sum _{k=0}^{t-p-1}{\nu }_k(s_0)\sum _{u=1}^m{\gamma }_u(s_0){\zeta }_{t-k}(s_u),\sum _{\ell =0}^{t^{\prime }-p-1}{\nu }_{\ell }(s_0)\sum _{v=1}^m{\gamma }_v(s_0){\zeta }_{t^{\prime }-\ell }(s_v)\right).$$

$$\text{Cov}({\hat{Z}}_t^2(s_0),{\hat{Z}}_{t^{\prime }}^2(s_0))=\sum _{k=0}^{t-p-1}\sum _{\ell =0}^{t^{\prime }-p-1}{\nu }_k(s_0){\nu }_{\ell }(s_0)\sum _{u=1}^m\sum _{v=1}^m{\gamma }_u(s_0){\gamma }_v(s_0)\text{Cov}({\zeta }_{t-k}(s_u),{\zeta }_{t^{\prime }-\ell }(s_v)).$$

$$\text{Cov}({\hat{Z}}_t^2(s_0),{\hat{Z}}_{t^{\prime }}^2(s_0))=\sum _{k=0}^{t-p-1}\sum _{\ell =0}^{t^{\prime }-p-1}{\nu }_k(s_0){\nu }_{\ell }(s_0)\sum _{u=1}^m\sum _{v=1}^m{\gamma }_u(s_0){\gamma }_v(s_0)R_{uv}.$$

$${\mathcal{L}}_T\left(s_0,\Theta \right)=\frac{1}{m{T}^{\prime }}\sum _{t=p+1}^T\sum _{u=1}^m{W}_b\left(s_0-s_u\right)I_t(\Theta ),$$

$$I_t(\Theta )=\frac{1}{2}\left(\log {\sigma }_t^2(s_u)+\frac{Z_t(s_u{)}^2}{{\sigma }_t^2(s_u)}\right),$$

$$\int K(x)\mathrm{d}x=1,W(s)=K(x)K(y)\text{ and }{W}_b(.)=\frac{1}{b}W\left(\frac{.}{b}\right).$$

Full text

1 Introduction

Volatility clustering is a well-known phenomenon in financial time series: large return innovations tend to be followed by periods of high volatility, while tranquil periods often coincide with reduced market activity. Generalised autoregressive conditional heteroskedasticity (GARCH) models (Engle,, 1982; Bollerslev,, 1986) provide a robust framework for capturing such temporal dependence in volatility. However, in real-world financial systems, individual assets are not isolated. Instead, they are embedded in complex financial networks shaped by common exposures, interbank lending, or structural economic similarities (see, e.g., Diebold and Yilmaz,, 2023; Chen et al.,, 2019; Mattera and Otto,, 2024). Consequently, volatility patterns often exhibit pronounced dependence across firms, sectors, or regions. This phenomenon is typically referred to interchangeably as spatial, cross-sectional, network, or topological dependence, depending on the context and the underlying representation of relationships. Similar forms of volatility propagation arise in environmental and epidemiological contexts, where uncertainty in model predictions, such as pollutant concentrations or disease incidence rates, varies dynamically across both space and time due to complex interactions between spatial covariates, temporal dynamics, and measurement errors (Huang et al.,, 2011; Otto et al.,, 2018, 2023; Rodeschini et al.,, 2024).

This motivates the development of spatiotemporal GARCH models, which incorporate spatial interaction in the volatility dynamics. Existing approaches can be grouped into four broad classes: (1) multivariate GARCH models that include spatial lags in the volatility equation but assume cross-sectional independence (in volatility) at the same time point (e.g., Borovkova and Lopuhaa,, 2012; Caporin and Paruolo,, 2015; Hølleland and Karlsen,, 2020; Billio et al.,, 2023), (2) spatiotemporal GARCH or log-GARCH models with instantaneous spatial dependence in volatility (e.g., Sato and Matsuda,, 2017; Sato and Matsuda, 2018a, ; Sato and Matsuda, 2018b, ; Otto and Schmid,, 2023; Otto et al., 2024a, ), (3) stochastic volatility models with spatially correlated innovations (Huang et al.,, 2011), and (4) the novel idea proposed in this paper where spatial dependence enters through the conditional mean equation, indirectly propagating volatility through autoregressive feedback.

Most of the spatiotemporal GARCH models have been developed for areal data or fixed sets of spatial locations (a notable exception is the stochastic volatility model of Huang et al., (2011)), limiting their ability to predict volatility at unobserved sites. We refer the interested reader to Otto et al., (2025) for a review of spatial and spatiotemporal volatility models. In contrast, geostatistical models operate on continuous spatial domains and enable interpolation (kriging) of latent variables at arbitrary locations. While kriging is widely used in mean modelling, its extension to volatility models remains underexplored. Bridging this gap requires a framework that not only accounts for spatial volatility dynamics but also enables volatility prediction at unobserved locations, a task with clear relevance in financial risk monitoring, stress testing, and systemic risk assessment, but also in other fields such as environmental science (e.g., modelling and forecasting pollution variability across monitoring networks) or epidemiology, where it supports the reliable estimation of spatial prediction intervals of disease outbreak intensities or uncertainty in infection dynamics.

Another limitation of existing models is their assumption of spatially homogeneous GARCH parameters, which may be unrealistic in heterogeneous domains. In financial networks, for instance, firms with different balance sheet structures may exhibit distinct levels of baseline volatility, persistence, and responsiveness to shocks. Similar forms of spatial heterogeneity arise in environmental applications (e.g., urban vs. rural air quality dynamics) and epidemiology (e.g., region-specific disease transmission patterns). Allowing GARCH parameters to vary smoothly across space reflects these differences and provides a more flexible and interpretable representation of volatility dynamics. Such parameter surfaces have been previously explored in nonparametric time series models (Dahlhaus and Subba Rao,, 2006), and we extend this idea to a spatiotemporal GARCH framework. Most existing spatiotemporal GARCH models rely on the assumption of spatial stationarity, which may be overly restrictive in the presence of heterogeneity, asymmetry, or directional effects (Aouri et al.,, 2025). To address this, we build on the probabilistic results of Aouri et al., (2025), who extend classical GARCH dynamics by allowing spatially varying parameters and introducing geostatistical dependence via a covariance structure on the innovation process. This captures contemporaneous cross-sectional correlations that propagate forward in time through the recursive GARCH mechanism, enabling a parsimonious yet flexible representation of spatial volatility spillovers.

In this paper, we propose a heterogeneous spatiotemporal GARCH model defined over a continuous spatial domain. Volatility is modelled using a standard GARCH($p$,$q$) recursion at each spatial location, with spatially varying parameters and geostatistical dependence in the innovation process. Estimation is based on a locally weighted quasi-likelihood approach that yields smooth parameter surfaces and enables kriging of the volatility at unobserved locations. We illustrate the proposed methodology in an empirical application to a panel of 50 global firms with daily stock returns. More precisely, we consider the spatial domain to be represented by firm-level financial indicators, such as balance sheet positions, which induce a latent economic space over which the GARCH process is defined. This setting allows us to model and predict volatility across financial networks, also for unobserved locations/firms. We demonstrate how the resulting model captures spatial heterogeneity in volatility, propagates risk through the network, and enables accurate prediction of volatility at unobserved firms. Our findings highlight the importance of both spatial dependence and heterogeneity in modelling firm-level risk dynamics.

The remainder of the paper is organised as follows. Section 2 introduces the heterogeneous spatiotemporal GARCH model. Sections 2.2 and 2.3 present the local quasi-likelihood estimation procedure and discusses kriging for volatility prediction. In Section 3, we evaluate the model’s performance through a Monte-Carlo simulation study. Section 4 provides the empirical analysis of financial network data. Section 5 concludes with a discussion and outlook.

2 Heterogeneous Spatiotemporal GARCH model

We consider a spatiotemporal process $\{Z_{t}(s):s\in\mathcal{D},\,t\in\mathcal{T}\}$, defined over a general spatial domain $\mathcal{D}\subseteq\mathcal{M}$, where $\mathcal{M}$ denotes an underlying spatial structure such as a Euclidean space, a manifold (e.g., a sphere), or a network. This formulation allows for broad applicability: for instance, $\mathbb{R}^{2}$ is suitable for most environmental or epidemiological data with geographic coordinates, the sphere $\mathbb{S}^{2}$ accommodates global climate models, and financial or transportation networks are naturally represented by graph-based domains, where distances reflect structural dissimilarity rather than physical proximity. To ensure validity across such domains, the spatial covariance function governing the innovations must respect the geometry of $\mathcal{M}$ (see, e.g., Gneiting et al.,, 2010; Porcu et al.,, 2016, 2023), and the local estimation kernels must be adapted accordingly. In the present work, without loss of generality, we focus on a Euclidean domain $\mathcal{D}=[0,1]^{2}\subset\mathbb{R}^{2}$ and assume a discrete, equidistant time index $t\in\mathcal{T}=\mathbb{Z}$. We define a spatiotemporal GARCH($p$,$q$) (STGARCH) process as follows:

[TABLE]

where $\omega(s)$, $\alpha_{i}(s)$, and $\beta_{j}(s)$ are location-specific parameters that vary smoothly over space, with the orders $q<p$. The innovation process $\{\eta_{t}(s)\}$ is assumed to be independent across time, strictly stationary and ergodic, with zero mean and unit variance, and exhibits spatial dependence through a valid covariance function on $\mathcal{D}$. In what follows, we consider a covariance function of the Matérn class (Gneiting et al.,, 2010). This geostatistical structure captures contemporaneous cross-sectional correlations in the innovation process, which subsequently propagate into future volatility via the recursive GARCH equation. Thus, the model allows for spatial heterogeneity in both the innovation structure and the volatility dynamics, while supporting prediction at unobserved locations. Under mild smoothness conditions on the parameter functions, the spatially nonstationary process $Z_{t}(s)$ can be locally approximated at any fixed location $s_{0}\in[0,1]^{2}$ by a spatially stationary GARCH($p$,$q$) process with constant parameters evaluated at $s_{0}$. In particular, define

[TABLE]

This locally stationary approximation underpins our asymptotic analysis in Section 2.4. Further details on the approximation can be found in Appendix A.2.

Remark 1**.**

Although the volatility process $\sigma_{t}^{2}(s)$ in (1) does not include instantaneous spatial interactions (i.e., there are no terms of the form $\sigma_{t}^{2}(s^{\prime})$ with $s^{\prime}\neq s$), and is, thus, temporally causal at each location $s$, it nonetheless inherits spatial dependence through the squared observations $Z_{t-i}^{2}(s)$. Since $Z_{t}(s)=\sigma_{t}(s)\eta_{t}(s)$ and the innovation process $\eta_{t}(s)$ is spatially correlated across $s\in\mathcal{D}$, the process $Z_{t}(s)^{2}$ also exhibits spatial dependence. Importantly, the range of spatial dependence in $Z_{t}(s)^{2}$ differs from that in $\eta_{t}(s)$. For instance, if $\eta_{t}(s)$ follows a zero-mean, unit-variance Gaussian process with spatial covariance function $C(h)$, it holds that $\mathrm{Cov}(\eta_{t}(s)^{2},\eta_{t}(s^{\prime})^{2})=2\cdot C(h)^{2},$ where $h=\|s-s^{\prime}\|$. Thus, the spatial correlation in the squared process $\eta_{t}(s)^{2}$ decays faster than in $\eta_{t}(s)$, since $C(h)^{2}$ drops off more quickly with increasing $h$ than $C(h)$. This effect carries over to the squared process $Z_{t}(s)^{2}$, and consequently to $\sigma_{t}^{2}(s)$ through the recursive GARCH dynamics. Therefore, even though $\sigma_{t}^{2}(s)$ is not explicitly defined via spatial interactions at the same time point, it is still a spatially dependent process, with the dependence structure shaped by the properties of the innovation field.

2.1 Model Assumptions

To develop the asymptotic theory for our spatiotemporal GARCH process defined in Equation (1), we impose a set of assumptions that ensure temporal stationarity, control the spatial variation in parameters, and guarantee the existence of necessary moments. We begin by formalising the properties of the innovation process.

Assumption 1**.**

The innovations $\eta_{t}(s)$ are independent over time, and spatially, strictly stationary.

This assumption reflects the temporal white noise nature of the innovations, while allowing for spatial dependence, captured through a valid geostatistical covariance function. Stationarity in space ensures that the structure of the innovation process is spatially invariant, an important requirement for spatial prediction and asymptotic consistency. Next, we define the parameter space and regularity conditions that ensure the existence of a strictly stationary solution to the spatiotemporal GARCH recursion.

Assumption 2**.**

Let $\Omega$ be the compact set defined by:

[TABLE]

where $0<\rho_{1}\leq\rho_{2}<\infty$ and $\gamma$ is a positive constant. For each $s\in(0,1]^{2}$, we assume that $\Theta$ is in the interior of $\Omega$, where $\Theta=(\omega(s),\alpha_{1}(s),\ldots,\alpha_{p}(s),\beta_{1}(s),\dots,\beta_{q}(s))$.

This assumption guarantees that the local GARCH dynamics are well-defined, with the stationarity condition that $\sum\alpha_{i}(s)+\sum\beta_{j}(s)$ is bounded away from the unit circle, uniformly across space. It is worth noting that this condition also implies that $\inf_{s\in\mathcal{D}}\sigma_{t}^{2}(s)>\rho_{1}$, i.e., it is uniformly bounded away from zero, to avoid degeneracy of the conditional variance. Moreover, the parameter surfaces need to be smooth across space. This smoothness is crucial for the local estimation procedure and for ensuring uniform convergence of estimators. This assumption allows us to approximate the process locally by a stationary GARCH model and facilitates the application of kernel-based smoothing methods in estimation.

Assumption 3**.**

The second-order partial derivatives of $\omega_{t}(s)$, $\alpha_{i}(s)$ and $\beta_{j}(s)$ exist and all their partial derivatives up to the second order are uniformly bounded over $[0,1]^{2}.$

Typically for stationary GARCH processes, the result can be proved if the fourth moments of the innovations are finite. We require a mildly stronger moment condition on the innovation process that is slightly stronger than the standard GARCH setting. This is necessary to apply central limit results for martingale arrays, which arise due to the locally weighted quasi-likelihood estimation (cf. Hall and Heyde,, 1980, , Theorem 3.2).

Assumption 4**.**

For some $v>0$, $E\left(\left|\eta_{i}(s)\right|^{4+v}\right)<\infty$.

The assumption 2, 3 and 4 are needed to prove asymptotic normality of the estimators defined below. In geostatistical applications, a common choice for the innovation process $\eta_{t}(s)$ is a Gaussian process with mean zero and a Matérn covariance function, which satisfies all required conditions. In this case, all moments exist and Assumption 4 is trivially fulfilled.

2.2 Prediction at Unobserved Locations

Before turning to estimation, we first outline the kriging problem for spatially nonstationary spatiotemporal GARCH models as defined in equation (1). Our goal is to predict the process $Z_{t}(s_{0})$ at an arbitrary, unobserved location $s_{0}\in\mathcal{D}$, based on the available observations $\{Z_{t}(s_{u}):t=1,\dots,T;u=1,\dots,m\}$. This prediction, known as spatial interpolation or kriging in geostatistics, requires knowledge of the location-specific parameters $\omega(s_{0}),\alpha_{i}(s_{0}),\beta_{j}(s_{0})$ that define the conditional variance $\sigma_{t}^{2}(s_{0})$. While these parameter surfaces are not known in practice, it is conceptually helpful to first develop the prediction framework under the assumption that they are available. This approach mirrors classical kriging setups, where the spatial covariance structure is assumed known in the derivation of the best linear unbiased predictor. In Section 2.3, we return to the problem of estimating the parameters $\omega(\cdot),\alpha_{i}(\cdot),\beta_{j}(\cdot)$ via a locally weighted likelihood approach. For now, we focus on the construction of the kriging predictor and its associated uncertainty.

To facilitate prediction, we first derive an ARMA-type representation of the STGARCH model, which enables straightforward computation of linear predictors. We assume that the innovations $\eta_{t}(s)$ follow a spatially stationary Gaussian process, consistent with the common geostatistical modelling approach. For each fixed $s$, the process $\{Z_{t}(s)\}$ then constitutes a stationary GARCH($p$,$q$) time series under the conditions in Assumption 2. More precisely, this implies that the squared process $Z_{t}^{2}(s)$ satisfies the following stationary ARMA representation, i.e.,

[TABLE]

where $r=\max(p,q)$, and $\delta_{i}(s)=\alpha_{i}(s)+\beta_{i}(s)$, with $\alpha_{i}(s)=0$ for $i>q$ and $\beta_{i}(s)=0$ for $i>p$. The innovation term is defined as $\zeta_{t}(s)=\sigma_{t}^{2}(s)(\eta_{t}^{2}(s)-1)$, which is uncorrelated over time and has zero mean under the assumption $\mathbb{E}[\eta_{t}^{2}(s)]=1$ (follows straightforward from the unit variance of $\eta_{t}(s)$).

Hence, we can write $Z_{t}^{2}(s)$ in the following form

[TABLE]

and the sequence of non-parametric functions $\phi_{j}(s)$ is given by

[TABLE]

where $\beta_{0}(s)=1$ , $\beta_{j}(s)=0$ for $j>p$ , $\phi_{j}(s)=0$ if $j<0$ . In addition

[TABLE]

Thus, $Z_{t}^{2}(s)$ can be rewritten as follows:

[TABLE]

where $\nu_{k}(s)=[A(s)^{k}]_{1,1}$, with

[TABLE]

and we assume that

[TABLE]

Let $\mathbf{Z}=\left[\left\{Z_{j}^{2}(s_{u}):u=1,\ldots,m\right\}_{j=1}^{T}\right]$ denote the observed squared process across all spatial locations and time points. Based on the ARMA representation in equation (7), we evaluate the conditional expectation $\mathbb{E}\left[Z_{t}^{2}\left(s_{0}\right)\mid\mathbf{Z}\right]$ for $p+1\leq t\leq T$. To this end, we consider the vector-valued process

[TABLE]

where

[TABLE]

To derive the kriging predictor for the latent volatility process at an unobserved location $s_{0}$, we consider the conditional expectation $\mathbb{E}[Z_{t}^{2}(s_{0})\mid\mathbf{Z}]$, which provides the minimum mean-square error prediction of $Z_{t}^{2}(s_{0})$ given the observed spatiotemporal process. Starting from the ARMA representation in equation (3), we can express this conditional expectation as

[TABLE]

Here, $\mathbb{E}[\zeta_{t-k}(s_{0})\mid\mathbf{Z}]$ is evaluated by conditioning on $\{\zeta_{t-k}(s_{u}):u=1,\dots,m\}$ using the spatial covariance structure of the innovation process, which enables interpolation of volatility innovations across space. Specifically, for $p+1\leq t-k\leq T$, we express

[TABLE]

and substitute it into the expression above to obtain

[TABLE]

Under Assumption 2, we have the bound $\left\|A\left(s_{0}\right)^{k}\right\|_{\text{spec }}\leq K\lambda^{k}$, where $\|\cdot\|_{\text{spec}}$ denotes the spectral norm, which means

[TABLE]

Further details are provided in Appendix A.2.

Next, we demonstrate that the first term on the right-hand side of (16) can be evaluated recursively. Suppose that the random process $\left\{Z_{t}^{2}\left(s_{0}\mid m\right)\right\}$ satisfies

[TABLE]

where $Z_{t}^{2}\left(s_{0}\mid m\right)=0$ for $t\leq p$. We define the $p$-dimensional vectors

[TABLE]

We note that (2.2) can be written as $\underline{Z}_{t}^{2}\left(s_{0}\mid m\right)=A\left(s_{0}\right)\underline{Z}_{t-1}^{2}\left(s_{0}\mid m\right)+\underline{\zeta}_{t}\left(s_{0}\mid m\right)+C(s_{0})$. Given the initial condition $\underline{Z}_{p}^{2}\left(s_{0}\mid m\right)^{\prime}=(0,\ldots,0)$, the solution for $Z_{t}^{2}\left(s_{0}\mid m\right)$ is

[TABLE]

Comparing $Z_{t}^{2}\left(s_{0}\mid m\right)$ and (16), we observe that $Z_{t}^{2}\left(s_{0}\mid m\right)$ and $\mathbb{E}\left[Z_{t}^{2}\left(s_{0}\right)\mid\mathbf{Z}\right]$ are asymptotically equivalent as $t$ becomes large. As a result, $Z_{t}^{2}\left(s_{0}\mid m\right)$ can be utilised as an estimator for $\mathbb{E}\left[Z_{t}^{2}\left(s_{0}\right)\mid\mathbf{Z}\right]$.

Evaluating $Z_{t}^{2}\left(s_{0}\mid m\right)$ requires determining $\mathbb{E}\left[\zeta_{t}\left(s_{0}\right)\mid\zeta_{t}\left(s_{1}\right),\ldots,\zeta_{t}\left(s_{m}\right)\right]$, along with the non-parametric functions $\omega(\cdot),\alpha_{j}(\cdot),\text{and}\;\beta_{i}(\cdot)$. Furthermore, the conditional expectation $\mathbb{E}\left[\zeta_{t}\left(s_{0}\right)\mid\left\{\zeta_{t}\left(s_{u}\right)\right\}_{u=1}^{m}\right]$ can only be computed in closed form if the distribution of $\left\{\zeta_{t}\left(s_{u}\right):u=0,\ldots,m\right\}$ is known. As the Gaussian process is the most widely adopted distribution in geostatistics, and facilitates closed-form derivations, we consider this special case in more details below. In the case where the innovation process $\{\eta_{t}(s)\;;\;s\in[0,1]^{2}\}$ follow a spatial Gaussian process, then $\zeta_{t}(s)$ is a spatial centred $\chi^{2}$-distributed process. We can demonstrate that

[TABLE]

[TABLE]

where

[TABLE]

with

[TABLE]

This implies that by utilising estimates of the innovations $\left\{\zeta_{t}\left(s_{u}\right)\right\}_{t}$ at the observed locations, we can estimate the covariance function $c_{\zeta}\left(\cdot\right)$, e.g., using a maximum-likelihood approach. Recall that, as noted in Remark 1, if $\eta_{t}(s)$ follows a Gaussian process with covariance $C_{\eta}(h)$, then $\eta_{t}^{2}(s)-1$ has covariance $2\,C_{\eta}(h)^{2}$, implying a shorter effective spatial range. Consequently, $\zeta_{t}(s)$ inherits this reduced range whenever $\sigma_{t}^{2}(s)$ varies smoothly over space. In practice, we recommend to back-transform the estimated $\zeta_{t}(s)$ to $\eta_{t}(s)$, i.e., $\eta_{t}(s)=\text{sign}(Z_{t}(s))\sqrt{\frac{\zeta_{t}(s)}{\sigma_{t}^{2}(s)}+1}$.

Remark 2**.**

In practice, the quantities $\omega(\cdot)$, $\alpha_{j}(\cdot)$, and $\beta_{i}(\cdot)$ required in (2.2) are unknown and must be estimated from the data. One possibility is to estimate all parameter surfaces jointly with the spatial covariance parameters of $\zeta_{t}(s)$ in a single optimisation step, which fully accounts for the spatio–temporal dependence but can be computationally demanding in high–dimensional settings. Alternatively, we suggest adopting a computationally simpler two–stage strategy: we first estimate the local spatiotemporal GARCH parameter surfaces independently at each location, then compute the residuals $\{\eta_{t}(s)\}$, and finally estimate the spatial covariance parameters from these residuals. This approach is common in geostatistical regression and retains good statistical properties while being substantially faster than the joint estimation alternative.

Finally, we derive the covariance of the kriging predictor for the squared spatiotemporal GARCH process $Z_{t}^{2}(s_{0})$ at an unobserved spatial location $s_{0}$. Based on the recursive representation and spatial kriging approximation of the innovations, the predictor is given by

[TABLE]

Since $C(s_{0})$ is nonrandom, it does not contribute to the covariance. Due to the temporal independence of $\zeta_{t}(\cdot)$, and if $t-k=t^{\prime}-\ell$ we get

[TABLE]

which simplifies to

[TABLE]

The resulting expression yields the closed-form covariance of the kriging predictor for $Z_{t}^{2}(s_{0})$, determined by the temporal ARMA coefficients $\nu_{k}(s_{0})$, the spatial kriging weights $\gamma_{u}(s_{0})$, and the spatial covariance matrix $R$ of the innovation process. Since these quantities depend on unknown parameter surfaces $\omega(\cdot)$, $\delta_{i}(\cdot)$, $\beta_{j}(\cdot)$ as well as on the spatial covariance of $\zeta_{t}(s)$, their practical implementation requires consistent estimation from the observed data. The following section outlines the estimation strategy for these components, which provides the foundation for constructing the kriging predictor in applied settings.

2.3 Estimation of Parameter Surface

To estimate the spatially varying parameters of the spatiotemporal process, we employ a nonparametric regression framework based on localised estimation. Since $Z_{t}(s_{0})$ is unobserved at arbitrary locations $s_{0}$, classical approaches such as ordinary least squares are not directly applicable. Nevertheless, assuming smooth spatial variation of the model parameters (Assumption 3), nonparametric techniques allow their estimation at unobserved locations.

Inspired by local likelihood methods in nonparametric regression Tibshirani and Hastie, (1987) and the locally stationary framework developed by Dahlhaus and Subba Rao, (2006), we adopt a weighted quasi-maximum likelihood approach. This method incorporates spatial weights, giving higher importance to observations near the prediction location, and allows for flexible estimation of spatially varying GARCH parameters while accounting for the complexities of spatiotemporal dependence.

Given the smoothness assumption in Assumption 3, we estimate the parameters at $s_{0}$ using a localised likelihood approach that weights observations from nearby locations more heavily. Specifically, we define the local log-likelihood

[TABLE]

with

[TABLE]

where $\Theta=\left(\omega\left(s_{0}\right),\alpha_{1}(s_{0})\ldots,\alpha_{p}\left(s_{0}\right),\beta_{1}(s_{0}),\dots\beta_{p}(s_{0})\right)$ and ($\alpha_{q+1}=0,\alpha_{q+2}=0,\dots,\alpha_{p}=0)$, $T^{\prime}$ is the number of effective time points, i.e., $T^{\prime}=T-p-1$.

The spatial weighting function is constructed from a kernel $K:\mathbb{R}^{2}\rightarrow\mathbb{R}$ satisfying

[TABLE]

The local parameter estimate $\hat{\Theta}(s_{0})$ is then obtained as

[TABLE]

This localised estimation procedure is computationally far more efficient than a full joint estimation of all spatial and temporal parameters (see Remark 2) and scales well to large spatial datasets.

The localised likelihood in (24) provides pointwise estimates of the smoothly varying parameter functions at any location $s_{0}$. Having established the estimation framework, we next investigate its theoretical properties, focusing on consistency and asymptotic normality under the assumptions stated in Section 2.1.

2.4 Properties of the Estimator

In this section, we analyse the theoretical properties of the weighted quasi-maximum likelihood estimator $\hat{\Theta}(s_{0})$ for the spatiotemporal GARCH model. We work within a purely temporal asymptotic framework, treating the spatial grid of fixed size $m$ as fixed and letting the time series length $T\to\infty$. Under the previously stated regularity conditions, Theorem 1 establishes the consistency of $\hat{\Theta}(s_{0})$, while Theorem 2 derives its asymptotic normality.

Theorem 1** (Consistency).**

Let $\{Z_{t}(s_{u}):\,t=1,\dots,T,\;u=1,\dots,m\}$ be a spatiotemporal GARCH process satisfying Assumptions 2.1. Let

[TABLE]

be the local weighted quasi-likelihood estimator defined before. Suppose $m$ is fixed and $b\rightarrow 0\;\text{and}\;\;bmT\rightarrow\;\infty\;\text{as}\;\;T\rightarrow\infty,$ we have

[TABLE]

The proof can be found in the Appendix. Having established the consistency, we now turn to the distributional behaviour of the estimator. The next theorem shows that $\hat{\Theta}(s_{0})$ is asymptotically normal with a bias term of order $b^{2}$ due to the kernel smoothing.

Theorem 2** (Asymptotic Normality).**

Let $\{Z_{t}(s_{u}):\,t=1,\dots,T,\;u=1,\dots,m\}$ be a spatiotemporal GARCH process satisfying Assumptions 2.1, and let $\hat{\Theta}(s_{0})$ be the local quasi-likelihood estimator defined as

[TABLE]

*Then, we have the following:

[TABLE]

where:

[TABLE]

Together, Theorems 1 and 2 provide a rigorous justification for the localised likelihood estimation procedure introduced in the previous subsection. The second term in the asymptotic expansion of Theorem 2 reflects a bias that arises from the deviation of the true process from its stationary approximation over the localised spatial segment. This bias depends on the curvature of the spatial parameter functions and quantifies how closely the process resembles a spatially stationary one in the neighbourhood of $s_{0}$. It is therefore directly controlled by the bandwidth $b$, highlighting the bias–variance trade-off introduced by spatial kernel smoothing.

3 Monte Carlo Simulation Study

We assess the performance of the proposed STGARCH estimator and predictor through a Monte Carlo simulation with $MC=100$ replications. In each replication, we independently and uniformly sample $n_{1}\in\{50,100,150\}$ training sites ($S_{1}$) and $n_{2}=50$ prediction sites ($S_{2}$) from the unit square $[0,1]^{2}$. A fixed $10\times 10$ grid ($S_{0}$) is additionally used for evaluating parameter recovery on a dense spatial grid. We consider time series lengths $T\in\{100,200,300\}$, where each spatially varying coefficient function $\omega(s)$, $\alpha(s)$, and $\beta(s)$ is represented using a B-spline basis, with coefficients randomly drawn in each replication to ensure varying parameter surfaces across runs. The STGARCH(1,1) processes are simulated with spatially varying parameters and a spatially stationary Gaussian innovation process $\eta_{t}(s)\sim\mathcal{N}(0,1)$ following an exponential covariance function with range parameter $\theta=0.5$. The range parameter is estimated via the proposed weighted quasi-maximum likelihood method applied to residuals, yielding an average estimate of $\hat{\theta}=0.4992$ across replications (with $n_{1}=150$ and $T=300$).

For each $s_{i}\in S_{1}$, the STGARCH coefficient functions are estimated using weighted quasi-maximum likelihood with a uniform kernel and bandwidth $b=n_{1}^{-1/4}$. The estimated parameters are then used to spatially predict unconditional volatilities at the locations in $S_{2}$ via kriging, i.e., out-of-sample. The size of $S_{2}$ was chosen to be $n_{2}=50$. This is implemented as a two-step procedure: first, $\omega(s)$, $\alpha(s)$, and $\beta(s)$ are estimated at each $s_{i}$ via the weighted localised QMLE, and second, the innovations $\eta_{t}(s)$ and the derived $\zeta_{t}(s)$ are computed for use in the kriging procedure. For reporting accuracy, we compute Monte Carlo (MC) averages of bias and RMSE of all parameters across simulation runs and spatial locations. That is, we define

[TABLE]

These metrics are reported separately for $\theta\in\{\omega,\alpha,\beta\}$. Table 1 reports the average bias and RMSE for the estimated nonparametric functions $\omega$, $\alpha$, and $\beta$ on $S_{2}$, as well as for the predicted unconditional volatilities.

To further illustrate the predictive performance of the proposed method, Figure 1 compares the predicted and true unconditional volatility surfaces over the spatial grid $S_{0}$. The background colour gradient represents the predicted values, with red dashed contours indicating the predicted level sets and black solid contours representing the true ones. The close alignment of these contours confirms that the estimator successfully recovers the underlying spatial volatility pattern. Complementing this spatial view, Figure 2 presents the temporal evolution of the true and predicted conditional variances at three randomly selected out-of-sample locations in $S_{2}$. The model closely follows both the timing and magnitude of volatility spikes, demonstrating its strong ability to capture dynamic features of the process across time and space.

Overall, the simulation results demonstrate that the proposed STGARCH estimation and prediction framework reliably recovers spatially varying parameters and accurately predicts both unconditional and conditional volatilities across a range of spatial sample sizes and time series lengths. Bias and RMSE values decrease as either $n_{1}$ or $T$ increases, reflecting the expected gains from richer spatial and temporal information. Spatial reconstruction closely matches the true volatility surfaces, and temporal predictions successfully capture both the timing and magnitude of volatility changes at unobserved locations. These findings confirm that the localised estimation-kriging approach offers both high predictive accuracy and computational efficiency, making it a practical tool for large-scale spatiotemporal volatility modelling.

4 Volatility Prediction across Financial Networks

To demonstrate the utility of the proposed spatiotemporal GARCH model with spatially varying coefficients and spatially correlated innovations, we apply it to the problem of volatility modelling in a network of financial assets. Our setting is non-standard from a geostatistical perspective: the spatial domain is not physical space but an abstract “financial space” constructed from firm-specific balance sheet features. This allows us to embed firms into a latent space in which spatial proximity reflects economic similarity, and where spatial dependence can be exploited for volatility prediction.

The financial dataset consists of 50 firms traded on public markets, observed over 276 trading days in a single year. The full list of the 50 companies included in the analysis, spanning various sectors and geographic regions, is provided in Table 2. For each firm, 66 balance sheet variables were extracted. After removing variables with more than 5% missing values, 44 features remained. Each feature was standardised across firms, and missing values were imputed with zero, effectively corresponding to the assumption that the mean of each feature is zero in the standardised space. This ensures that the feature space is centred and that the Euclidean distances are interpretable as deviations from the average firm. To provide a descriptive overview of the input data, we visualised the distributions of daily log-returns and absolute log-returns for all 50 firms (see Figure 3). The boxplots in Figure 0a illustrate that log-returns are approximately symmetric around zero for most firms, but exhibit substantial heterogeneity in scale and tail behaviour. This variation motivates the use of firm-specific volatility dynamics in our modelling approach. Figure 0b displays the distribution of absolute log-returns, which serve as a proxy for realised volatility. Here, the differences across firms are even more pronounced, reflecting variation in baseline volatility levels, sensitivity to market conditions, or idiosyncratic shocks. These plots underscore the necessity of a model that allows for spatially varying volatility structures, as pursued in our spatiotemporal GARCH specification.

To reconstruct the financial network structure, we followed the strategy proposed by (Fülle and Otto,, 2024), where firms are embedded in a low-dimensional proxy space based on selected balance sheet positions. Since no true spatial network is observable in this context, the model’s performance depends heavily on how well the constructed space reflects economic relationships relevant to volatility transmission. Instead of selecting a single financial proxy a priori, we systematically explored all pairwise combinations of the 44 features to define candidate two-dimensional spaces. Each of the resulting 946 spaces defines a distinct notion of “distance” among firms. For each space, we fitted the spatiotemporal GARCH model and evaluated its out-of-sample performance in predicting volatility.

Volatility predictions were assessed using five-fold cross-validation. In each fold, the model parameters were estimated on the full time series of a subset of firms and used to predict volatility on the held-out firms (i.e., we left the entire time series out). While this random allocation of firms into folds does not reflect realistic market clustering, and may lead to a mild overestimation of predictive performance Otto et al., 2024b , block cross-validation is not feasible in this setting as the underlying network structure is unknown. Moreover, for a fair comparison, we kept the fold assignments constant for all 946 pairwise proxy spaces.

Figure 4 displays the ranked RMSE values across all 946 proxy spaces. The RMSE is computed between the predicted volatility $\hat{z}^{2}_{t}(s)$ and the absolute observed log-return $|y_{t}(s)|$, providing a direct measure of how well the conditional standard deviations reflect the magnitude of return innovations. Among all candidate spaces, the best prediction performance was obtained for the space defined by the variables other current liabilities and long-term debt, with an RMSE of approximately 0.0112. This proxy space reflects the structure of firms’ liabilities—both short- and long-term—and may be interpreted as a financial exposure dimension. Economically, it is plausible that firms with similar liability profiles are jointly sensitive to changes in credit conditions or interest rates, which in turn influence their volatility responses to market shocks.

Other high-performing spaces include combinations such as total cash from operating activities with other current liabilities, and net income with total cash flows from investing activities. These pairings suggest that profitability and internal cash generation are also relevant for explaining co-movement in volatility. Spaces involving net borrowings and various forms of net income point to the role of financing strategy and retained earnings in shaping a firm’s exposure to systemic volatility.

While several of these spaces result in near-identical RMSE values, we will, for the sake of interpretability and visual clarity, focus our analysis on the optimal two-dimensional space defined by other current liabilities and long-term debt. In principle, the analysis could be extended to three- or four-dimensional spaces, potentially lowering the RMSE further and capturing additional economic structure. However, such extensions come at the cost of reduced interpretability and visualisation clarity, especially when aiming to represent spatial variation in parameter surfaces and volatility fields.

Figure 5 provides a joint representation of the temporal and spatial structure of the predicted volatilities obtained from cross validation study. The left panel shows the predicted volatility trajectories over time, with grey lines representing individual firm-level volatilities and the bold black line denoting the cross-sectional median. This illustrates both market-wide and firm-specific volatility dynamics across the year 2019, with evident spikes around mid-year and increased volatility clustering. The right panel displays the firms in the optimal two-dimensional proxy space spanned by other current liabilities and long-term debt. This reveals substantial spatial heterogeneity in the volatility dynamics: firms with similar financial characteristics exhibit similar volatility levels and uncertainty, supporting the relevance of the proxy space for modelling and predicting financial risk.

To analyse the selected space in more detail, we re-estimated the model using all firms and all 276 days. The innovation covariance parameters were estimated via maximum likelihood, yielding $\hat{\theta}=80.3695$, $\hat{\sigma}^{2}=0.0573$, and $\hat{\tau}^{2}=0.17173$. These values suggest that the spatial correlation in the innovation process decays slowly. Indeed, the maximum distance among firms in the selected space was only 5.6. The large estimate of $\theta$ implies that nearly all firms are subject to common shocks, with little attenuation based on financial distance. This is further supported by the relatively large nugget effect $\tau^{2}$, which captures idiosyncratic variation not explained by the spatial process. While the variance parameters were statistically significant, the range parameter had a relatively large standard error, yielding a z-score of 1.419.

Figure 6 displays the estimated spatial surfaces of the three GARCH parameters, $\beta(s)$, $\alpha(s)$, and $\omega(s)$. To ensure stable and reliable estimation, we restrict the analysis to the convex hull of the financial proxy space. This is crucial because local likelihood estimators are known to perform poorly near the boundary, where neighbourhoods become asymmetric and effective sample sizes are insufficient.

The resulting parameter surfaces reveal substantial spatial heterogeneity in volatility dynamics across firms. The spatial variation of the parameters reveals distinct patterns across the financial proxy space. The persistence parameter $\beta(s)$ attains its highest values in the north-eastern and northern regions, corresponding to firms with high long-term debt and moderate to high current liabilities. This suggests that highly leveraged firms exhibit more persistent volatility dynamics, possibly due to slower adjustments in their financial risk profiles. In contrast, firms in the south-western region (characterised by low debt and liabilities) show lower persistence, indicating that shocks dissipate more rapidly. The innovation parameter $\alpha(s)$ is largest in the south-western region, highlighting that more conservatively financed firms respond more strongly to recent return shocks. Lastly, the unconditional volatility level $\omega(s)$ is elevated in the northern part of the space, again associated with high-debt firms, reflecting a higher baseline risk level. These spatial patterns offer economically meaningful insights into the relationship between firm-specific financial characteristics and the temporal dynamics of their volatility. Such heterogeneity in volatility dynamics is critical for modelling systemic risk in financial networks and may guide the development of portfolio allocation strategies and stress-testing frameworks that explicitly account for firm-specific balance sheet configurations.

5 Conclusion

We introduced a heterogeneous spatiotemporal GARCH framework on continuous spatial domains that permits spatially varying volatility parameters and contemporaneous cross-sectional dependence through a geostatistical innovation field. We developed a localised quasi–likelihood estimator that yields smooth parameter surfaces and a kriging scheme for predicting volatility at unobserved locations. Theoretical results established consistency and asymptotic normality under a fixed–$m$, $T\to\infty$ regime. In simulations, the method recovered parameter surfaces with low bias and RMSE and produced accurate spatial and temporal volatility predictions. In an empirical study on a financial proxy space spanned by balance-sheet variables, the model captured economically interpretable spatial heterogeneity in baseline risk, shock sensitivity, and persistence, and delivered competitive predictive performance while remaining computationally efficient via a two–stage estimation–kriging procedure.

Several limitations suggest natural extensions. The Gaussian–process assumption for the innovation field is convenient but not essential; accommodating heavy–tailed or skewed innovations would improve robustness to extreme events, which might be particularly useful for financial data. The fixed–domain asymptotics with increasing time simplify analysis but leave open theory for double asymptotics with $m\to\infty$ and $T\to\infty$, data–driven bandwidth selection, and boundary correction. Another promising direction for future research is to allow coefficients to vary in both space and time, i.e., $\omega(s,t)$, $\alpha(s,t)$, and $\beta(s,t)$, enabling local regime changes and structural breaks to be captured explicitly; this, in turn, calls for a space–time kriging framework for volatility that propagates uncertainty across both dimensions. Extending the spatial domain beyond Euclidean settings to networks and manifolds with covariance functions valid on graphs or spheres would broaden applicability to financial, transportation, and global climate systems. On the computational side, joint estimation with sparsity or low-rank structure, or adaptive bandwidths. Finally, for the empirical application, model selection for the proxy space, multivariate co–volatility extensions, and mechanisms for time-varying network topology offer promising avenues for capturing evolving interdependencies in financial markets and beyond.

Appendix A Appendix

A.1 Proofs of Estimator Properties

In this appendix, we collect the detailed proofs of the main asymptotic results stated in Section 2.4 . In particular, we first establish the consistency of the local weighted quasi-maximum likelihood estimator (Theorem 1), and then derive its asymptotic normality (Theorem 2).

Proof.

Proof of Theorem 1.

We recall the weighted quasi-maximum likelihood estimator at $s_{0}$ by

[TABLE]

where

[TABLE]

and

[TABLE]

Moreover, let $\Theta$ be the entire parameter surface across all $s\in\mathcal{D}$.

We also define the stationary–approximate criterion

[TABLE]

where

[TABLE]

and $\widetilde{Z}_{t}(s)$ is locally stationary B.1.

Step 1. Set

[TABLE]

We show

[TABLE]

Indeed,

[TABLE]

with

[TABLE]

[TABLE]

Under Assumption 2, stationarity, and ergodicity, and by the deviation bounds in Theorem B.1 and Corollary B.1, both $I_{1}$ and $I_{2}$ are $o(1)$.

Step 2. Let

[TABLE]

By the ergodic theorem in $t$ and the fact that $\sum_{u}W_{b}(s_{0}-s_{u})\approx 1$, we have for each fixed $\Theta$,

[TABLE]

A standard equicontinuity argument on the compact $\Omega$ then yields

[TABLE]

Step 3. Combining Steps 1–2 gives

[TABLE]

Step 4. Let $\Theta_{s_{0}}=\arg\min_{\Theta}L(s_{0},\Theta)$, which is unique by Lemma 5.5 of Horváth and Kokoszka, (2003). By definition of $\hat{\Theta}$ and uniform convergence,

[TABLE]

Hence $\hat{\Theta}_{s_{0}}\xrightarrow{\mathcal{P}}\Theta_{s_{0}}$.

∎

Proof.

Proof of Theorem 2.

We use Taylor expansion around $\Theta_{s_{0}}$:

[TABLE]

Since $\nabla\mathcal{L}\left(s_{0},\hat{\Theta}\right)=0$, we get

[TABLE]

We also have

[TABLE]

which implies

[TABLE]

and hence

[TABLE]

Replacing this in Equation 25, we obtain:

[TABLE]

or equivalently,

[TABLE]

We use $\nabla\tilde{\mathcal{L}}\left(s_{0},\Theta_{s_{0}}\right)$ to show asymptotic normality:

[TABLE]

where

[TABLE]

$\nabla I_{t}(s_{0},\Theta)$ is a martingale difference, which implies that $\nabla\tilde{\mathcal{L}}\left(s_{0},\Theta_{s_{0}}\right)$ is a weighted sum of martingale differences. Using the central limit theorem and verifying the conditions of Corollary 3.1 Hall and Heyde, (1980), we conclude that

[TABLE]

where

[TABLE]

As a result, we have

[TABLE]

Now, consider

[TABLE]

Using a Taylor expansion of $\nabla I_{t}(\Theta)$ around $s_{0}$, we get:

[TABLE]

Hence,

[TABLE]

where

[TABLE]

and

[TABLE]

is the second moment of the kernel, while

[TABLE]

Hence, we obtain:

[TABLE]

As a result, we conclude:

[TABLE]

∎

A.2 The Approximated Spatially Stationary Process

In this section, we investigate the probabilistic properties of the STGARCH process. Specifically, we demonstrate that the spatially nonstationary process $Z_{t}^{2}(s)$ can be locally approximated by a spatially stationary process $\widetilde{Z}^{2}_{t}(s)$ in the neighborhood of any fixed spatial location $s_{0}$. This approximation plays a critical role in establishing the asymptotic properties of the estimator presented in Appendix A.1. We begin with the following lemma, which is essential for proving that the process $Z^{2}_{t}(s)$can be locally approximated by a spatially stationary process.

Lemma B.1**.**

Let $Z_{t}(s)$ be STGARCH(p,q) model and satisfies the previous Assumptions 1 to 4. Suppose that for some $1-\gamma<\lambda<1$:

[TABLE]

We assume for all k and $C>0$

[TABLE]

Proof.

Since there exists a $1-\gamma<\lambda<1$ and finite $C$ such that for all $k$,

[TABLE]

we have

[TABLE]

Furthermore, since $\sup_{s}\left|\nu_{k}(s)\right|\leq K\sup_{s}\left\|A(s)^{k}\right\|,$ $(1)$ follows directly. Further, since $\mathbb{E}\left(\left|\eta_{t}^{4+v}(s)\right|\right)<\infty$, $(2)$ can easily be verified. We now prove $(3)$. Using $\nu_{k}(s)=[A(s)^{k}]_{(1,1)}$ and the norm inequalities, we have

[TABLE]

By the Lipschitz continuity of $\omega(s)$ , $\alpha_{i}(s)$ and $\beta_{j}(s)$ we have that

[TABLE]

Hence, by using $\sup_{s}\left\|A(s)^{k}\right\|_{\text{spec}}\leq C\lambda^{k}$ we obtain

[TABLE]

for some $C$ independent of $s,v$. This establishes $(3)$. Finally, we prove $(4)$ we consider the partial derivatives of $\nu_{k}(s)$. From the definition of $\nu_{k}(s)$ we have

[TABLE]

we require a bound for $\left\|\frac{A(s)^{k}}{\partial x}\right\|_{2}$. Expanding $\frac{A(s)^{k}}{\partial x}$ gives

[TABLE]

By using the uniform boundedness of $\frac{\omega(s)}{\partial x},\;\frac{\alpha_{i}(s)}{\partial x}\;\text{and}\;\frac{\beta_{j}(s)}{\partial x}$ we have

[TABLE]

we have $\left|\frac{\partial\nu_{k}(x,y)}{\partial x}\right|\leq Ck\lambda^{k-1}$. By applying the same method, we obtain the bounds for all partial derivatives of $\nu_{k}(s)$ up to the second order.

∎

The theorem presented below illustrates that the spatially stationary process $\widetilde{Z}^{2}_{t}(s)$ offers a local approximation of the spatially nonstationary process $Z^{2}_{t}(s)$ within a certain neighborhood of $s_{0}$.

Theorem B.1**.**

Let $Z_{t}(s)$ be STGARCH model and satisfies Assumptions 2.1. $\widetilde{Z}^{2}_{t}(s)$ is defined as follow:

[TABLE]

Then we have

[TABLE]

where for some $1-\gamma\leq\lambda<1$,

[TABLE]

$C$* and $C_{1}$ are a finite constant independent of s, and $W_{t}(s)$ is a stationary process in both space and time with a finite variance.*

Proof.

By using the MA( $\infty$ ) representation, we can evaluate an upper bound for $\left|Z_{t}^{2}(s)-\widetilde{Z}^{2}_{t}(s)\right|$. From the definition of $\widetilde{Z}^{2}_{t}(s)$, Lemma B.1(3) and Lipschitz continuity, we have

[TABLE]

[TABLE]

Therefore,

[TABLE]

where $W_{t}(s)=C\sum_{k=0}^{\infty}(k+1)^{2}\lambda^{k}|\zeta_{t-k}(s)|+C_{1}$. ∎

Corollary B.1**.**

Let $Z_{t}(s)$ be STGARCH model and satisfies Assumptions 2.1, and

[TABLE]

if

[TABLE]

then

[TABLE]

where

[TABLE]

$N_{t}(s)$* is a stationary process in both space and time with a finite variance.*