Heterogeneous Spatiotemporal GARCH Models

• Heterogeneous spatiotemporal GARCH models are a conditional heteroscedasticity framework that jointly captures temporal dynamics, spatial interactions, and location-specific heterogeneity.
• They incorporate diverse architectures—such as unified, vec-spARCH, and geostatistical forms—to model instantaneous spillovers and lagged variance effects in complex networks.
• Methodologies employ various estimation techniques and weight specifications, addressing regime dependence, asymmetric effects, and cross-variable interactions for improved volatility forecasting.

A heterogeneous spatiotemporal GARCH model is a conditional heteroscedasticity framework in which volatility depends jointly on temporal dynamics, spatial or network structure, and non-homogeneous mechanisms such as location-specific parameters, cross-variable effects, fixed effects, threshold regimes, or state dependence. In the recent literature, the label does not denote a single canonical recursion. Instead, it covers a family of related constructions: unified spatial and spatiotemporal ARCH/GARCH systems with instantaneous spillovers across spatial units (Otto et al., 2019), generalized spatial and spatiotemporal GARCH frameworks based on fixed-point equations (Otto et al., 2021), multivariate vec-form spatiotemporal ARCH models with contemporaneous spatial and lagged temporal variance interactions (Otto, 2022), dynamic spatial panel log-volatility models with fixed effects (Otto et al., 2022, Otto et al., 2023), geostatistical GARCH models with spatially varying coefficients and spatially correlated innovations (Aouri et al., 11 Aug 2025), and spatiotemporal E-GARCH extensions with asymmetric network spillovers (Chrisko et al., 7 Nov 2025). This suggests that the subject is best understood as a broad modeling paradigm rather than a single specification.

1. Lineage and conceptual scope

The modern literature begins from the observation that standard time-series GARCH captures volatility clustering over time but does not encode cross-sectional dependence in the conditional second moment. Early spatial ARCH-type work formalized the idea that local risk can cluster over space and time, and that the same basic representation,

$Y=\operatorname{diag}(h)^{1/2}\varepsilon,$

can be indexed over a purely spatial domain or over a joint space-time domain by taking $s=(s_s,t)'$ (Otto, 2018). The 2019 unified treatment then organized spatial ARCH, spatial GARCH, exponential spatial GARCH, and logarithmic spatial GARCH into a single framework, while explicitly stating that classical time-series GARCH appears as the one-dimensional special case (Otto et al., 2019). The 2021 generalization sharpened this perspective by writing the second-moment dynamics as a general functional equation and by emphasizing instantaneous spillovers across spatial units (Otto et al., 2021).

A second line of development moved from univariate fields to multivariate and panel settings. The dynamic spatiotemporal ARCH model of Otto et al. introduced contemporaneous spatial effects, temporal effects, spatiotemporal effects, and additive region and time fixed effects in a log-volatility equation (Otto et al., 2022). The multivariate vec-spARCH model then made same-time spatial spillovers, lagged ARCH effects, and cross-variable interactions explicit through a Kronecker-product representation (Otto, 2022). Later work added small- and large-sample comparisons of QMLE and GMM estimators in dynamic spatiotemporal ARCH panels (Otto et al., 2023), Markov-switching regime dependence in spatio-temporal log-ARCH systems (Khoo et al., 2023), and log-linear asymmetric spillovers in spatiotemporal E-GARCH models (Chrisko et al., 7 Nov 2025).

A third line reintroduced classical GARCH recursion directly at each spatial location. The heterogeneous STGARCH model keeps the familiar $\mathrm{GARCH}(p,q)$ update but lets the coefficients vary smoothly over a spatial or latent network domain, while contemporaneous dependence is introduced through a geostatistical covariance structure on the innovations rather than an explicit spatial lag in the variance recursion (Aouri et al., 11 Aug 2025). This differs materially from adjacency-matrix formulations. A plausible implication is that the literature now contains two major branches: direct spatial-lag volatility models and geostatistical models with spatially correlated innovations.

2. Canonical model architectures

The principal specifications differ in where space enters the volatility dynamics and in what is allowed to vary across locations, variables, or regimes.

Family	Representative form	Heterogeneity device
Unified spatial/spatiotemporal GARCH	$F=\alpha+W_1\gamma(Y^{(2)})+W_2F$	location-specific levels, weights, functions
Vec-spARCH	$$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$	variable-specific levels and cross-variable effects
Dynamic spatiotemporal log-ARCH	$$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$	region and time fixed effects
Heterogeneous STGARCH	$$\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)$$	spatially varying GARCH coefficients
Spatiotemporal E-GARCH	$$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$	node-specific intercepts and asymmetric spillovers

In the unified spatial and spatiotemporal GARCH framework, the central object is the transformed second-moment vector $F=(f(h(s_1)),\dots,f(h(s_n)))'$, with

$F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$

Choosing $s=(s_s,t)'$0 and $s=(s_s,t)'$1 yields spatial ARCH; keeping $s=(s_s,t)'$2 yields spatial GARCH,

$s=(s_s,t)'$3

and with $s=(s_s,t)'$4, formally

$s=(s_s,t)'$5

when the inverse exists (Otto et al., 2019, Otto et al., 2021). The $s=(s_s,t)'$6 term propagates squared shocks across locations, whereas $s=(s_s,t)'$7 is the spatial autoregressive feedback in the second moment. The literature repeatedly emphasizes that these models allow instantaneous spillovers across spatial units rather than forcing spatial effects to appear only with a time lag (Otto et al., 2021).

The multivariate vec-spARCH model is structurally different. For $s=(s_s,t)'$8, it specifies

$s=(s_s,t)'$9

or, in vec form,

$\mathrm{GARCH}(p,q)$0

Here, diagonal entries of $\mathrm{GARCH}(p,q)$1 and $\mathrm{GARCH}(p,q)$2 represent own-variable spatial and temporal ARCH effects, while off-diagonal entries represent cross-variable effects. The paper is explicit that spatial effects are contemporaneous and that Tobler’s law motivates this separation between same-time spatial interaction and lagged temporal dependence (Otto, 2022).

Dynamic spatiotemporal ARCH panel models begin from

$\mathrm{GARCH}(p,q)$3

and transform squared observations using $\mathrm{GARCH}(p,q)$4. The resulting equation,

$\mathrm{GARCH}(p,q)$5

contains a contemporaneous spatial lag $\mathrm{GARCH}(p,q)$6, a temporal lag $\mathrm{GARCH}(p,q)$7, and a spatiotemporal lag $\mathrm{GARCH}(p,q)$8, plus region and time fixed effects (Otto et al., 2022). This class is ARCH-based rather than GARCH-based in formal structure, but the papers explicitly position it as a generalized heterogeneous spatial volatility model with a GARCH-like interpretation.

The geostatistical heterogeneous STGARCH model restores a literal GARCH recursion at each location: $\mathrm{GARCH}(p,q)$9 Its distinctive feature is that $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$0 is spatially correlated, so there are no explicit contemporaneous spatial lag terms in the variance equation. Instead, dependence enters through the innovation field and propagates into $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$1 and then into $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$2 via the recursive GARCH structure (Aouri et al., 11 Aug 2025).

Finally, the spatiotemporal E-GARCH model embeds Nelson’s log-volatility asymmetry into the spatial-network setting: $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$3 with

$F=\alpha+W_1\gamma(Y^{(2)})+W_2F$4

Here $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$5 captures asymmetric instantaneous network spillovers, and $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$6 captures contemporaneous spatial volatility spillovers in log-variance (Chrisko et al., 7 Nov 2025).

3. What “heterogeneous” means in this literature

Heterogeneity is not used in a single sense. In the unified spatial GARCH literature, heterogeneity first means that the spatial units need not share the same baseline or the same neighborhood response. The intercept vector $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$7 may be location-specific, the weight matrices may vary by region or distance, and even the functions $F=\alpha+W_1\gamma(Y^{(2)})+W_2F$8 can differ across units in the general formulation (Otto et al., 2021). The 2019 unified paper further notes heterogeneous weighting schemes such as

$F=\alpha+W_1\gamma(Y^{(2)})+W_2F$9

or kernel-based weights $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$0, which permit anisotropic or distance-decay dependence (Otto et al., 2019).

In multivariate spatiotemporal ARCH, heterogeneity means that different variance channels are separately parameterized. The vec-spARCH paper states that “heterogeneous” means the model distinguishes and separately parameterizes own-variable instantaneous spatial effects, cross-variable instantaneous spatial effects, own-variable temporal ARCH effects, cross-variable temporal ARCH effects, and variable-specific unconditional variance levels (Otto, 2022). This is a channel-wise notion of heterogeneity rather than a regime notion.

In dynamic panel formulations, heterogeneity is additive in the log-volatility equation through region fixed effects $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$1 and time fixed effects $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$2. The papers explicitly allow these effects to be arbitrarily correlated with regressors (Otto et al., 2022, Otto et al., 2023). In this usage, heterogeneity refers to unobserved cross-sectional and temporal variation that is not captured by common coefficients.

Other models use regime or state mechanisms. The Markov-switching spatio-temporal ARCH model lets

$$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$3

with $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$4 following a Markov chain, so that different volatility regimes have different spatial, temporal, and spatiotemporal dependence parameters (Khoo et al., 2023). The Poisson threshold network GARCH model uses node- and time-varying autoregressive coefficients,

$$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$5

so different nodes may be in different threshold regimes at the same time (Pan et al., 2024). The spatiotemporal E-GARCH model introduces asymmetry through the leverage parameter $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$6, with $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$7 meaning negative shocks have a larger effect on volatility than positive shocks of the same magnitude (Chrisko et al., 7 Nov 2025).

A closely related, but not spatial, notion appears in the state heterogeneous GARCH-Itô model. There, intraday volatility evolves continuously in time but switches according to an exogenous low-frequency state $$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$8, and the daily integrated volatility has a regime-switching GARCH-like structure (Chun et al., 2021). This is not a spatiotemporal model, but it clarifies that heterogeneity can also be indexed by observed economic states rather than by spatial position or latent regimes. This suggests that the broader heterogeneous-volatility literature shares a common concern: replacing a single homogeneous recursion with parameter variation tied to state, location, variable, or network channel.

4. Estimation, identification, and asymptotic theory

Estimation strategies reflect the structural diversity of the models. For unified spatial and spatiotemporal ARCH/GARCH systems, nonlinear least squares is derived from the transformed relation

$$\ln vec(H_t)=vec(A)+(\Psi'\otimes W)\ln vec(Y_t^{(2)})+(\Pi'\otimes I)\ln vec(Y_{t-1}^{(2)})$$9

leading to criteria such as

$$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$0

with consistency proved under compactness, differentiability, finite-moment, and row-standardization conditions (Otto et al., 2021). The 2019 unified approach also derives maximum likelihood estimators for Type II models by transforming from $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$1 to $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$2, which introduces a Jacobian determinant and requires numerical optimization (Otto et al., 2019).

For vec-spARCH, the log-squared transformation yields a multivariate spatiotemporal autoregressive process,

$$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$3

where $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$4. Stability requires all eigenvalues of $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$5 to have modulus smaller than one. Under Assumptions 1–5, $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$6, $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$7, and $$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$8 are uniquely identifiable, and the Gaussian QMLE is consistent (Otto, 2022).

Dynamic spatiotemporal ARCH panels are estimated by GMM after orthonormal and demeaning transformations remove the fixed effects. The moment vector combines linear and quadratic moments,

$$Y_t^*=\rho_0MY_t^*+\gamma_0Y_{t-1}^*+\delta_0MY_{t-1}^*+X_t\beta_0+\mu_0+\alpha_{t0}1_n+\varepsilon_t^*$$9

and the optimal estimator minimizes $$\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)$$0. Under the stated assumptions, the GMM estimator is consistent and asymptotically normal both when $$\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)$$1 is finite and $$\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)$$2 and when $$\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)$$3 (Otto et al., 2022). The 2023 follow-up compares two QMLEs and a GMM estimator: the direct QMLE performs best for the spatial coefficient $$\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)$$4, the GMME performs best for the temporal coefficient $$\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)$$5, and the transformation QMLE can be notably biased when $$\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)$$6 (Otto et al., 2023).

The geostatistical heterogeneous STGARCH model uses local weighted quasi-likelihood at a target location $$\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)$$7,

$$\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)$$8

with

$$\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)$$9

Consistency is proved as $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$0 with fixed $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$1, $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$2, and $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$3, and asymptotic normality includes a kernel-smoothing bias of order $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$4 (Aouri et al., 11 Aug 2025).

The spatiotemporal E-GARCH model requires numerical inversion because contemporaneous spatial shock terms make the mapping from innovations to observables nontrivial. Its QML estimator uses recursive inversion of $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$5, Jacobian determinants, and inverse-Hessian standard errors; the paper recommends discarding the first five observations because initial conditions affect numerical inversion (Chrisko et al., 7 Nov 2025). The Markov-switching spatio-temporal ARCH model, by contrast, uses Hamilton filtering, Kim smoothing, and QMLE under Gaussian quasi-errors (Khoo et al., 2023).

Across the literature, stability conditions are central. Representative examples include $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$6 for simple spGARCH (Otto et al., 2021), $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$7 for local STGARCH (Aouri et al., 11 Aug 2025), and

$$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$8

for spatiotemporal E-GARCH (Chrisko et al., 7 Nov 2025). A recurring theme is that spatial feedback alters the usual time-series stationarity logic: stability depends not only on temporal parameters but also on the spatial operator and the chosen weight matrix.

5. Empirical domains and comparative evidence

Empirical work spans real estate, financial networks, stock indices, and count processes on networks. The Berlin real-estate application in the vec-spARCH paper uses monthly returns for undeveloped land, developed land, and condominiums across $$\log(h_t)=\alpha_1+\rho_0W_1g(\varepsilon_t)+\rho_1g(\varepsilon_{t-1})+\lambda_0W_2\log(h_t)+\lambda_1\log(h_{t-1})$$9 postcode regions and $F=(f(h(s_1)),\dots,f(h(s_n)))'$0 months. The main result is that temporal effects dominate spatial effects; spatial spillovers are positive but weak, and interactions between property types occur only in the temporally lagged variables (Otto, 2022). The dynamic spatiotemporal ARCH application to Berlin condominium returns from January 1995 to December 2015 reports statistically significant temporal, spatial, and spatiotemporal lags, with $F=(f(h(s_1)),\dots,f(h(s_n)))'$1, $F=(f(h(s_1)),\dots,f(h(s_n)))'$2, and $F=(f(h(s_1)),\dots,f(h(s_n)))'$3, and interprets the positive $F=(f(h(s_1)),\dots,f(h(s_n)))'$4 as instantaneous spatial volatility spillovers across adjacent postcodes (Otto et al., 2022).

Financial-network applications show more varied patterns. The heterogeneous STGARCH model is applied to 50 publicly listed firms over 276 daily observations, using a latent proxy space built from 44 standardized balance-sheet features and evaluating all $F=(f(h(s_1)),\dots,f(h(s_n)))'$5 candidate two-dimensional proxy spaces. In five-fold cross-validation, the best-performing proxy space uses other current liabilities and long-term debt, with best RMSE about 0.0112 when predicted volatility is compared with absolute log returns (Aouri et al., 11 Aug 2025). The spatiotemporal E-GARCH model is fitted to NYSE/NASDAQ, DAX, and CAC 40 stocks over the COVID period; temporal parameters $F=(f(h(s_1)),\dots,f(h(s_n)))'$6, asymmetry parameters, and the contemporaneous spatial GARCH term $F=(f(h(s_1)),\dots,f(h(s_n)))'$7 are significant, while the contemporaneous spatial shock coefficient $F=(f(h(s_1)),\dots,f(h(s_n)))'$8 is often small or insignificant (Chrisko et al., 7 Nov 2025).

Regime dependence has also been empirically important. The Markov-switching spatio-temporal ARCH model is applied to 26 Asian stock indices and 2 U.S. indices from 2014 to 2019, using a Piccolo-distance-based financial network. The two-regime model improves the BIC relative to the one-regime counterpart, and the smoothed probabilities identify regime changes around September 2014, February 2015 to June 2015, June 2015, and January 2017 (Khoo et al., 2023).

Comparative evidence is especially informative. An empirical study on 16 sectorally diversified S&P 500 stocks from 22 December 1998 to 20 October 2024 compares DCC, BEKK variants, proximity-structured BEKK, Dynamic Spatiotemporal ARCH, circular STGARCH, spGARCH-X, and STEGARCH under multiple network specifications. The Dynamic Spatiotemporal ARCH / Dynamic Network Log-ARCH model achieves the lowest RMSFE and MAFE across almost all network specifications, with the best reported configuration under $F=(f(h(s_1)),\dots,f(h(s_n)))'$9: RMSFE $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$0 and MAFE $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$1 (Chrisko et al., 2 Mar 2026). The same paper reports that STEGARCH performs poorly across all networks and that model rankings are supported by pairwise Diebold–Mariano tests. This suggests that, in forecasting exercises, parsimonious heterogeneous spatiotemporal ARCH-type structures can dominate both standard GARCH benchmarks and more flexible but numerically fragile spatiotemporal GARCH variants.

A related empirical benchmark, though not spatial, comes from the SG-Ito model for S&P 500 intraday volatility. Using 2015–2018 high-frequency data and daily pre-averaging realized volatility, the paper finds significant state heterogeneity associated with leverage effects, investor attention or trading volume, market illiquidity, stock market comovement with the Hong Kong/Chinese market, and pre-holiday effects, while post-holiday effects are weak or absent (Chun et al., 2021). Although the model is not spatial, it reinforces the broader empirical message that heterogeneous volatility mechanisms become more detectable when richer data and more structured state variables are used.

6. Methodological issues, misconceptions, and open questions

One recurring misconception is that all heterogeneous spatiotemporal GARCH models are full GARCH systems in the narrow time-series sense. The literature does not support this. Several central contributions are explicitly ARCH, log-ARCH, or E-GARCH analogues rather than standard GARCH recursions with lagged conditional variances (Otto et al., 2022, Khoo et al., 2023). Conversely, the geostatistical STGARCH model is a literal $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$2 construction, but it does not use explicit spatial lag terms in the variance recursion (Aouri et al., 11 Aug 2025). The term therefore spans formally distinct architectures.

A second unresolved issue concerns the interpretation of $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$3. The unified spatial GARCH literature notes that in simultaneous spatial settings $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$4 need not always be the exact conditional variance in the classical sequential sense; it is more generally a latent second-moment quantity that behaves like local conditional variance (Otto et al., 2019). This is not merely terminological. It affects how likelihoods are constructed, how causality is interpreted, and how one compares spatial models with ordinary time-series GARCH.

The choice of spatial or network weights remains a major practical problem. The papers repeatedly state that $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$5 is crucial and usually must be specified in advance (Otto et al., 2021). Comparative evidence shows that forecasting performance depends strongly on the network specification, with Granger-filtered and correlation-based matrices often outperforming alternatives in financial applications (Chrisko et al., 2 Mar 2026). The heterogeneous STGARCH paper addresses this by abandoning a fixed adjacency matrix in favor of a continuous proxy space estimated from firm characteristics (Aouri et al., 11 Aug 2025). This suggests that weight specification is not a secondary tuning choice but a primary modeling decision.

Several technical limitations are model-specific. Log-squared transformations create zero-value problems, so real data may require ad hoc adjustments or dedicated zero models (Otto, 2022). In dynamic spatiotemporal ARCH panels, QMLEs can require bias correction and stronger growth-rate conditions than GMM (Otto et al., 2023). In spatiotemporal E-GARCH, recursive inversion and Jacobian evaluation can be numerically delicate (Chrisko et al., 7 Nov 2025). In count-valued network GARCH, asymptotics require both $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$6 and $F=\alpha+W_1\gamma(Y^{(2)})+W_2F.$7 large together with weak dependence and decay conditions on network weights (Pan et al., 2024). In state-heterogeneous high-frequency models, future states may be unobserved, making transition probabilities an essential extension (Chun et al., 2021).

The current literature therefore presents heterogeneous spatiotemporal GARCH not as a settled single model, but as a rapidly differentiating class of volatility systems. The common principle is stable across formulations: volatility is allowed to vary across space, time, network position, or economic state, and dependence may be contemporaneous, lagged, asymmetric, or regime-specific. The main substantive payoff has been the ability to represent spatial spillovers, temporal clustering, and heterogeneous persistence within a single volatility framework, while the main methodological challenge remains the joint treatment of network specification, identification, and computational tractability.