Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
119,875 characters · 8 sections · 68 citation commands
A Dynamic Spatiotemporal and Network ARCH Model with Common Factors
\singlespacing
JEL-Classification: C11, C23, C58.\\ Keywords: Spatial ARCH, volatility clustering, volatility, spatial dependence, common factors, Bayesian estimation. \onehalfspacing
Volatility modeling plays a crucial role in understanding the dynamics of economic, financial, and other complex systems, as it captures the uncertainty and the risk associated with the temporal evolution of an underlying outcome variable. The most common example of heteroskedastic time series can be found in the financial domain, albeit evidence of dynamic volatility has been found in macroeconomics kiani2004business, house prices bollerslev2016daily,dougan2023bayesian, and also in environmental time series huang2011class,otto2023dynamic. Traditional approaches model volatility dynamics through autoregressive processes bollerslev1986generalized that often rely on the information of volatility proxies and observable exogenous variables. In the multivariate setting, multivariate GARCH models are adopted to study the relationships between volatilities and correlations among a set of time series. These models are helpful under the assumption that the volatility of a time series is transmitted to another so that the information on the volatility of each time series is used to model the volatility of the others. However, these models are often not suitable for high-dimensional setups, as they usually occur in spatial (i.e., having a large number of spatial locations) or network settings (i.e., having a large number of nodes) because long time series are required for estimation. Instead, spatiotemporal and network ARCH models assume a certain structure in the (instantaneous) cross-sectional dimension (spatial or network dimension), which can be interpreted in the standard GARCH-type sense otto2024spatial.
The dependencies between time series can be significantly strengthened by incorporating information on spatial proximity. Tobler's First Law of Geography -- “everything is related to everything else, but near things are more related than distant things” -- suggests that the volatility dynamics of time series associated with spatial units that are geographically close are likely to exhibit stronger interdependencies than those that are farther apart. This implies that spatial proximity amplifies the connections between time series, intensifying the impact of local shocks and spillovers on volatility patterns. Traditional multivariate GARCH models, such as the BEKK engle1995multivariate or the DCC engle2002dynamic, do not allow for testing the existence of this kind of spillovers between the multivariate time series under study, nor quantifying precisely its magnitude. For financial time series, these instantaneous interactions in the volatility result from contemporaneous decisions of actors on the financial market. In other words, agents observe the financial risk of all (or most) financial assets and instantaneously account for these risks in their own trading decisions. Similar findings have been found for house prices cohen2016spatial,otto2018spatiotemporal. In the same manner, for environmental applications, environmental conditions affect the variation in all locations contemporaneously, with a decreasing magnitude for an increasing geographical distance rodeschini2024scenario,otto2023dynamic. In the case of macroeconomics, for example, there is widely documented evidence of the existence of spillovers between neighboring economies in terms of many indicators benos2015proximity,dell2013spatial,otieno2024public.
Network spillovers, however, offer a broader perspective on interactions compared to spatial spillovers. Unlike spatial models, which rely on geographic proximity, network spillovers consider complex linkages between entities that are connected through various forms of relationships, such as similar economic indicators in the case of countries, industry sectors in the case of stocks, or any other source of shared characteristics, that is irrespective of the geographic location. Proximity-based networks extend, therefore, the concept of similarity by defining 'proximity' in a more abstract sense beyond geographic distance fulle2024spatial. In these networks, nodes that are proximate in terms of such similarities are likely to experience stronger volatility spillovers. By incorporating proximity-based networks into volatility models, we capture a richer and more nuanced view of how shocks propagate, offering critical insights into the broader dynamics of interconnected systems. Second, the network structure can be defined by the proximity/correlation of the observed time series at each node of the network diebold2014network,mattera2024network. Thirdly, the network may be defined due to its physical nature, such as river networks tonkin2018role or energy networks malinovskaya2023statistical.
Despite spatial and/or network interactions, other approaches have been proposed to improve the fit of volatility models. For instance, empirical evidence in many fields suggests that better estimates of volatility dynamics are obtained by including additional latent factors. barigozzi2016generalized,barigozzi2017generalized,uddin2023risk shows that in large panels of financial time series with dynamic factor structure on the levels or returns, the volatilities of the common and idiosyncratic components are exposed to the same volatility shocks. Therefore, estimating and including these latent factors in the model improves the estimation and the prediction of the volatility for the multivariate time series. In a different applicative domain, gorodnichenko2017level identified separated factors for levels and volatility in macroeconomic time series, such as inflation and industrial production. Furthermore, christoffersen2019factor provided evidence that the factor structure in daily commodity volatility is much stronger than the factor structure in returns for commodities.
In this study, we introduce a new dynamic spatiotemporal model to specify the volatility process of an outcome variable. We assume that the log volatility of the underlying outcome depends on the log-squared outcome of other entities in the network, the time lag of the log-squared outcome, and the time lag of the log-squared outcome of other entities in the network. We refer to the influence of other entities' log-squared outcomes as the spatial effect (the contemporaneous volatility spillover), the influence of the time-lag term as the temporal effect, and the influence of the time-lag of other entities' log-squared outcomes as the spatiotemporal effect (the lagged volatility spillover). In addition to these effects, we also account for the influence of observed and unobserved factors on the log-volatility process. To capture the impact of latent factors, we employ a factor structure that allows for entity-specific effects and permits arbitrary correlations between the latent and observed factors.
Our approach allows us to identify all these effects, and we propose an efficient estimation method based on the Bayesian MCMC framework. We begin by applying a log-squared transformation to derive an estimation equation that is linear in the parameters and the log-squared original error terms Taspinar:2021, Robinson:2009. Following the time series literature Kim:1998, Omari:2007, we approximate the distribution of the log-squared error terms in the transformed model using a finite mixture of normal distributions. To facilitate Bayesian estimation, we express the mixture density using an auxiliary latent random variable that serves as a mixture indicator. Finally, we augment the transformed model with the latent mixture indicators and suggest an efficient Gibbs sampler for estimation. Importantly, our approach enables us to recover the volatility estimates for each entity in the network over time as a natural by-product of the estimation process.
We recommend two Bayesian methods for determining the number of latent factors in the model. The first method is the deviance information criterion (DIC) proposed by David:2002, which we demonstrate how to apply in our context. The DIC measures predictive performance, balancing model complexity and model fit. The second approach is based on the Bayesian Lasso, as considered by Park:2008. In this approach, we assume shrinkage priors for factor loadings to automatically exclude the latent factors that have less impact on the log-volatility process. In a simulation study, we show that our Bayesian estimation approach performs well and produces point estimates that are very close to true parameter values.
In two empirical applications, we use our model to estimate the volatility processes of the US house price returns at the state level and the US stock market returns. Our results from both applications indicate the presence of economy-wide latent factors that affect the volatility of both the US house price returns and US stock market returns. For the US house returns, our findings show that the spatial, temporal, and spatiotemporal effects estimates are positive and statistically significant. For the US stock market returns, we find that the spatial and spatiotemporal effects estimates are positive and statistically significant, with relatively larger magnitudes.
The rest of the paper is structured as follows. Section (ref) discusses the proposed dynamic spatiotemporal ARCH model with common factors. The interpretation of the model in terms of network process is therein discussed. Section (ref) shows the posterior analysis of the model, entering in the technicalities of the two algorithms proposed, while the simulation study is presented in detail in Section (ref). Section (ref) provides two applications, one for house prices in the US and another one for the US stock market. While the first application adopts a purely spatial perspective, the second one shows the application of our method from a network perspective. Section (ref) concludes with final remarks and future research directions.
We consider the random process $\{Y_t(\mathbf{s}): \mathbf{s}\in \mathbf{D}_1\subseteq\mathbb{R}^d, d>1, t\in D_2\subseteq\mathbb{R}\}$, where $\mathbf{s}$ and $t$ denote the spatial location and time point, respectively. The data-type determines the structure of the spatial domain $\mathbf{D}_1$ and the time domain $D_2$, and we assume that $\mathbf{D}_1=\{\mathbf{s}_1,\hdots,\mathbf{s}_n\}$ and $D_2=\{1,2,\hdots,T\}$. Then, we consider the following data generating process for $Y_t(\mathbf{s}_i)$:
where $h_{t}(\mathbf{s}_i)$ is considered as the volatility term in location $\mathbf{s}_i$ at time $t$, and $\{\operatorname{\varepsilon}_{it}(\mathbf{s}_i)\}$ is a sequence of i.i.d normal random variables with mean zero and unit variance. The log-volatility terms follow the process specified in (ref), where $\{m_{ij}\}$, for $i,j=1,\hdots,n$, are the non-stochastic weights specifying the strength of links among spatial units, $\mathbf{x}_{t}(\mathbf{s}_i)$ is a $k\times1$ vector of exogenous variables with the associated parameter vector $\boldsymbol{\beta}$, $\mathbf{f}_t$ is the $q\times1$ vector of unobserved common factors, and $\boldsymbol{\lambda}(\mathbf{s}_i)$ is the corresponding $q\times1$ vector of factor loadings. In the log-volatility equation, the spatial, temporal, and spatiotemporal effects of the log-squared outcome variable on the log-volatility are measured by the unknown scalar parameters $\gamma$, $\rho$, and $\delta$, respectively. We assume that the initial value vector $\mathbf{Y}_0=(Y_{0}(\mathbf{s}_1),\hdots,Y_{0}(\mathbf{s}_n))^{'}$ is observable.
The factor loadings and factors can be correlated with $\mathbf{x}_{t}(\mathbf{s}_i)$ in an arbitrary manner. This specification generalizes the dynamic spatiotemporal model suggested by Otto:2023 because the common factors specification $\boldsymbol{\lambda}^{'}(\mathbf{s}_i)\mathbf{f}_t$ can include the additive time and spatial fixed effects. For example, by setting $\mathbf{f}_t=(1,\alpha_{t0},0,\hdots,0)^{'}$ and $\boldsymbol{\lambda}(\mathbf{s}_i)=(\mu_0(\mathbf{s}_i),1,0\hdots,0)^{'}$, then we will obtain $\boldsymbol{\lambda}^{'}(\mathbf{s}_i)\mathbf{f}_t=\mu_0(\mathbf{s}_i)+\alpha_{t0}$, which is the additive structure assumed in Otto:2023. Let $\boldsymbol{\Lambda}=(\boldsymbol{\lambda}_1,\hdots,\boldsymbol{\lambda}_n)^{'}$ be the $n\times q$ matrix of loadings. It is well known that the factors and loadings cannot be identified separately Bai:2003. For example, for any orthogonal $q\times q$ matrix $\mathbf{Q}$, we have an observational equivalent relation $\boldsymbol{\Lambda}\mathbf{f}_t=\boldsymbol{\Lambda}\mathbf{Q}\mathbf{Q}^{'}\mathbf{f}_t$. Thus, to identify $\boldsymbol{\Lambda}$ and $\mathbf{f}_t$, we need to impose $q^2$ restrictions to pin down $\mathbf{Q}$. Our ensuing analysis will show that the estimation of the log-volatility term and other parameters only requires the estimation of the product $\boldsymbol{\Lambda}\mathbf{f}_t$. That is, our analysis does not require the identification of $\boldsymbol{\Lambda}$ and $\mathbf{f}_t$ separately.
Define $Y^*_{t}(\mathbf{s}_i)=\log Y^2_{t}(\mathbf{s}_i)$, $h^{*}_{t}(\mathbf{s}_i)=\log h_{t}(\mathbf{s}_i)$ and $\operatorname{\varepsilon}^*_{t}(\mathbf{s}_i)=\log\operatorname{\varepsilon}^2_{t}(\mathbf{s}_i)$. Then, we apply the log-squared transformation to (ref) and obtain
Using vector notation, we can express (ref) and (ref) as
where $\mathbf{M}=(m_{ij})$ is the $n\times n$ network matrix, $\mathbf{Y}^{*}_t=(Y^{*}_{t}(\mathbf{s}_1),\hdots,Y^{*}_{t}(\mathbf{s}_n))^{'}$, $\mathbf{h}^{*}_t=(h^{*}_{t}(\mathbf{s}_1),\hdots,h^{*}_{t}(\mathbf{s}_n))^{'}$, $\boldsymbol{\operatorname{\varepsilon}}^{*}_t=(\operatorname{\varepsilon}^{*}_{t}(\mathbf{s}_1),\hdots,\operatorname{\varepsilon}^{*}_{t}(\mathbf{s}_n))^{'}$, and $\mathbf{X}_t=(\mathbf{x}_{t}(\mathbf{s}_1),\hdots,\mathbf{x}_{t}(\mathbf{s}_n))^{'}$. Then, we obtain an estimation equation by substituting (ref) into (ref):
The estimation equation in (ref) is in the form of a spatial dynamic panel data model with common factors, as considered by Shi:2017 and Bai:2021, with two notable new features. Firstly, the outcome variable, the spatial lag term ($\mathbf{M}\mathbf{Y}^{*}_t$), and the spatiotemporal lag term ($\mathbf{M}\mathbf{Y}^{*}_{t-1}$) are formulated in terms of the log-squared outcome variable. Secondly, the elements of $\boldsymbol{\operatorname{\varepsilon}}^{*}_t$ are the log-squared original disturbance terms, i.e., $\operatorname{\varepsilon}^*_{t}(\mathbf{s}_i)=\log\operatorname{\varepsilon}^2_{t}(\mathbf{s}_i)$. If we assume that $\operatorname{\varepsilon}_{t}(\mathbf{s}_i)\sim N(0,1)$, then $\operatorname{\varepsilon}^{*}_{t}(\mathbf{s}_i)\sim\log\chi^2_1$, which is the log-chi squared distribution with one degree of freedom with the following density:\footnote{In our ensuing analysis, we use $p(\cdot)$ to denote density functions and omit $\mathbf{X}_t$ in the conditional set for the sake of simplicity.}
The first two moments of $\log\chi^2_1$ are $\mathbb{E}(\operatorname{\varepsilon}^{*}_{t}(\mathbf{s}_i)) = -\vartheta - \log(2) \approx-1.2704$, where $\vartheta$ is Euler's constant, and $\text{Var}(\operatorname{\varepsilon}^{*}_{t}(\mathbf{s}_i))=\pi^2/2\approx4.9348$ Peter:2012. In Figure (ref), we compare this density with that of $N(-\vartheta - \log(2),\,\pi^2/2)$. As seen from the figure, the log-chi squared distribution exhibits significant skewness with a long left tail.
Let $\mathbf{S}(\rho)=\mathbf{I}_n-\rho\mathbf{M}$ and $\boldsymbol{\varphi}=(\gamma,\delta)^{'}$. Then, under the assumption that $\mathbf{S}^{-1}(\rho)$ exists, the reduced form of the model can be expressed as
where $\mathbf{A}(\rho,\boldsymbol{\varphi})=\mathbf{S}^{-1}(\rho)\left(\gamma\mathbf{I}_n+\delta\mathbf{M}\right)$. Thus, if all eigenvalues of $\mathbf{A}(\rho,\boldsymbol{\varphi})$ lie within the unit circle, then $\mathbf{Y}^{*}_t$ will have a stable process Hamilton:1994. By the spectral radius theorem, $\mathbf{Y}^{*}_t$ has a stable process if $\Vert\mathbf{S}(\rho)\Vert<1$ and $\Vert\mathbf{A}(\rho,\boldsymbol{\varphi})\Vert<1$, where $\Vert\cdot\Vert$ is any matrix norm. If we use the matrix row sum norm and assume that $\mathbf{M}$ is row normalized, then it can be shown that a sufficient condition for stability is $|\rho|+|\gamma|+|\delta|<1$ Yu:2008.
We note that our specification can be applied to model more general data observed across (random) networks rather than exclusively spatio-temporal data fulle2024spatial, mattera2024network. A network $G$ is defined as the tuple $G=(V,E)$, where $V$ is a set of $n$ distinct vertices/nodes that are connected by directed or undirected edges, comprised in the set $E$. In practice, the edges are often represented by a (potentially sparse) $n$-dimensional adjacency matrix. In addition, the edges could be weighted to accommodate different “lengths/distances” between the nodes. Let $V=\mathbf{D}_1$ and define $\{m_{ij}\}$ the network weights representing the strength of the connection between two nodes $i$ and $j$. These weights can either be known in advance or estimated based on cross-sectional features fulle2024spatial or temporal similarities mattera2024network. In this case, we refer to the process as Network log-ARCH with common factors. While the technical details remain the same in the case of a network structure, the choice of different network definitions influences the interpretation of the results. This idea could also be extended to random graphs with a fixed set of nodes (potentially including isolated nodes), where the weight matrices $\mathbf{M}$ are allowed to vary over time.
For posterior analysis, we assume the following independent priors:
where $\text{Uniform}(c_1,\,c_2)$ is the uniform distribution over the interval $(c_1,\,c_2)$. When $\mathbf{M}$ is row normalized, $\mathbf{S}^{-1}(\rho)$ exists for $\rho\in(-1,1)$. Therefore, we assume that $\rho\sim\text{Uniform}(-1,1)$. When $\mathbf{M}$ is not row-normalized, $\mathbf{S}^{-1}(\rho)$ exists for $\rho\in(-1/\tau,1/\tau)$, where $\tau$ is the spectral radius of $\mathbf{M}$ KP:2010. Thus, we can alternatively assume that $\rho\sim\text{Uniform}(-1/\tau,1/\tau)$. Additionally, we also require that the prior distributions of $\rho$ and $\boldsymbol{\varphi}$ are subject to the stability condition stated in Section (ref). The prior distributions assumed for $\varphi$ and $\boldsymbol{\beta}$ are commonly used priors in the literature. Finally, following Geweke:1996 and Han:2016, we assume that $\mathbf{f}_t\sim N(\mathbf{0},\mathbf{I}_q)$, where $\mathbf{I}_q$ is the $q\times q$ identity matrix.
Following the time series literature on the standard stochastic volatility models Kim:1998, Omari:2007, we use a finite mixture of normal distributions to approximate $p(\operatorname{\varepsilon}^*_{t}(\mathbf{s}_i))$:
where $\phi(\operatorname{\varepsilon}^*_{t}(\mathbf{s}_i)|\mu_j,\,\sigma^2_j)$ denotes the Gaussian density function with mean $\mu_j$ and variance $\sigma^2_j$, and $p_j$ is the probability of $j$th mixture component. The parameters of this mixture distribution are given in Table (ref). These parameters are chosen by matching the first four moments of the ten-component normal mixture distribution with that of $p(\operatorname{\varepsilon}^*_{t}(\mathbf{s}_i))$. Under (ref), our model becomes a mixture of normal distributions. Furthermore, this approach does not pose any estimation difficulties, as the parameters of the mixture distribution are pre-determined.
\definecolor{LightCyan}{rgb}{0.88,1,1}
Let $Z_{t}(\mathbf{s}_i)$ be the mixture component indicator variable that takes values from $\{1,2,\hdots,10\}$. Then, using $Z_{t}(\mathbf{s}_i)$, we can express the finite mixture distribution in (ref) as
where $\mathbb{P}(Z_{t}(\mathbf{s}_i)=j)=p_j$ is the probability that $Z_{t}(\mathbf{s}_i)$ equals $j$. Let $\mathbf{Z}_t=(Z_{t}(\mathbf{s}_1),\hdots,Z_{t}(\mathbf{s}_n))^{'}$, $\mathbf{d}_t=(\mu_{Z_{t}(\mathbf{s}_1)},\hdots,\mu_{Z_{t}(\mathbf{s}_n)})^{'}$ and $\boldsymbol{\Sigma}_t=\operatorname{Diag}(\sigma^2_{Z_{t}(\mathbf{s}_1)},\hdots,\sigma^2_{Z_{t}(\mathbf{s}_n)})$, where $\operatorname{Diag}(a_1,\hdots,a_n)$ denotes the $n\times n$ diagonal matrix with the $i$th diagonal element $a_i$. Then, the mixture component indicator variable allows us to determine $\boldsymbol{\operatorname{\varepsilon}}^*_t|\boldsymbol{Z}_t\sim N(\mathbf{d}_t,\,\boldsymbol{\Sigma}_t)$. Thus, we can express the likelihood function of (ref) as
where $\mathbf{b}_{Y_t}=\mathbf{W}_t\boldsymbol{\varphi}+\mathbf{X}_t\boldsymbol{\beta}+\boldsymbol{\Lambda}\mathbf{f}_t+\mathbf{d}_t$ with $\mathbf{W}_t=(\mathbf{Y}^{*}_{t-1},\mathbf{M}\mathbf{Y}^{*}_{t-1})$. Let $\mathbf{Y}^{*}=(\mathbf{Y}^{*'}_1,\hdots,\mathbf{Y}^{*'}_T)^{'}$, $\mathbf{Z}=(\mathbf{Z}^{'}_1,\hdots,\mathbf{Z}^{'}_T)^{'}$, and $\mathbf{F}=(\mathbf{f}^{'}_1,\hdots,\mathbf{f}^{'}_T)^{'}$. Then, the joint posterior distribution $p(\mathbf{Z},\boldsymbol{\beta},\boldsymbol{\Lambda},\mathbf{F},\boldsymbol{\varphi},\rho|\mathbf{Y}^{*})$ can be expressed as
where $p(\mathbf{Y}^{*}|\mathbf{Z},\boldsymbol{\beta},\boldsymbol{\Lambda},\mathbf{F},\boldsymbol{\varphi},\rho)=\prod_{t=1}^Tp(\mathbf{Y}^{*}_t|\mathbf{Z}_t,\boldsymbol{\beta},\boldsymbol{\Lambda},\boldsymbol{f}_t,\boldsymbol{\varphi},\rho)$ and $p(\mathbf{a})$ denotes the prior distribution of $\mathbf{a}$ for $\mathbf{a}\in \{\mathbf{Z},\boldsymbol{\beta},\boldsymbol{\Lambda},\mathbf{F},\boldsymbol{\varphi},\rho\}$. We note that only the conditional posterior distribution of $\rho$ does not take a known form because the likelihood function is non-linear in $\rho$ as shown in (ref). For this parameter, we resort to the random walk Metropolis-Hastings (MH) algorithm suggested by Lesage:2009.\footnote{It is also possible to sample $\rho$ and $\boldsymbol{\varphi}$ jointly by using the adaptive MH algorithm considered in Haario:2001 and Roberts:2009. The main feature of the adaptive MH algorithm is that historical draws are used to determine the mean and covariance of the proposal distribution. Han:2016 and Yang:2023 use this algorithm for static and dynamic spatial panel data models.} The conditional posterior distributions for the remaining parameters take known forms and can be determined by the standard Bayesian analysis used for a linear regression model. We summarize how these distributions can be determined in Algorithm (ref).
In Algorithm (ref), the factors and loadings enter the conditional posterior distributions of $\boldsymbol{\beta}$, $\boldsymbol{\varphi}$ and $\rho$ as a product term $\mathbf{\Lambda}\mathbf{f}_t$. Moreover, the volatility estimates can be recovered from $\widehat{\text{$\mathbf{h}$}}^{*}_t=\widehat{\text{$\rho$}}\mathbf{M}\mathbf{Y}^{*}_t+\widehat{\text{$\gamma$}} \mathbf{Y}^{*}_{t-1}+\widehat{\text{$\delta$}}\mathbf{M}\mathbf{Y}^{*}_{t-1}+\mathbf{X}_t\widehat{\text{$\boldsymbol{\beta}$}}+\widehat{\text{$\boldsymbol{\Lambda}\mathbf{f}$}}_t$, where $\widehat{\text{$\rho$}}$, $\widehat{\text{$\gamma$}}$, $\widehat{\text{$\delta$}}$, $\widehat{\text{$\beta$}}$ and $\widehat{\text{$\boldsymbol{\Lambda}\mathbf{f}$}}_t$ are the estimated posterior means obtained through Algorithm (ref). Thus, our estimation approach does not require the identification of $\mathbf{f}_t$ and $\boldsymbol{\Lambda}$ separately.
The Gibbs sampler in Algorithm (ref) requires that the number of factors $q$ is known. The literature suggests several ways that can be used to determine the number of factors Ng:2002, Alexei:2010. Following Han:2016, we suggest using the deviance information criterion (DIC) suggested by David:2002 to determine the number of factors. The DIC can be considered as the Bayesian version of AIC and take the following from:
where the expectation is taken with respect to the joint posterior distribution and $ \sbox{\myboxA}{$\m@th\mathbf{a}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi\in\{ \sbox{\myboxA}{$\m@th\mathbf{Z}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\beta}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\Lambda}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\mathbf{F}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\varphi}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\rho$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi\}$ denotes the posterior mean. Let $\{\mathbf{Z}^{(g)},\boldsymbol{\beta}^{(g)},\boldsymbol{\Lambda}^{(g)},\mathbf{F}^{(g)},\boldsymbol{\varphi}^{(g)},\rho^{(g)}\}_{g=1}^R$ be a sequence of posterior draws generated through Algorithm (ref). Then, we can compute $\mathbb{E}\left(\ln p(\mathbf{Y}^{*}|\mathbf{Z},\boldsymbol{\beta},\boldsymbol{\Lambda},\mathbf{F},\boldsymbol{\varphi},\rho)|\mathbf{Y}^{*}\right)$ in the DIC measure by $\mathbb{E}\left(\ln p(\mathbf{Y}^{*}|\mathbf{Z},\boldsymbol{\beta},\boldsymbol{\Lambda},\mathbf{F},\boldsymbol{\varphi},\rho)|\mathbf{Y}^{*}\right)\approx\frac{1}{R}\sum_{g=1}^R\ln p(\mathbf{Y}^{*}|\mathbf{Z}^{(g)},\boldsymbol{\beta}^{(g)},\boldsymbol{\Lambda}^{(g)},\mathbf{F}^{(g)},\boldsymbol{\varphi}^{(g)},\rho^{(g)})$. The second term $\ln p(\mathbf{Y}^{*}| \sbox{\myboxA}{$\m@th\mathbf{Z}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\beta}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\Lambda}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\mathbf{F}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\boldsymbol{\varphi}$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi, \sbox{\myboxA}{$\m@th\rho$} \setbox\myboxB\null \ht\myboxB=\ht\myboxA \dp\myboxB=\dp\myboxA \wd\myboxB=0.5\wd\myboxA \sbox\myboxB{$\m@th\overline{\copy\myboxB}$} \setlength\mylenA{\the\wd\myboxA} \addtolength\mylenA{-\the\wd\myboxB} \ifdim\wd\myboxB<\wd\myboxA \rlap{\hskip 0.5\mylenA\usebox\myboxB}{\usebox\myboxA} \else \hskip -0.5\mylenA\rlap{\usebox\myboxA}{\hskip 0.5\mylenA\usebox\myboxB} \fi)$ in (ref) is simply the log-likelihood function evaluated at the posterior means. The DIC is a measure of predictive performance. Thus, a model with the smallest DIC value performs relatively better in terms of predictive accuracy. It is important to note that the DIC also does not require the separate identification of $\mathbf{f}_t$ and $\boldsymbol{\Lambda}$, as the factors and factor loadings are incorporated as a product term in the likelihood function.
An alternative approach for determining the number of factors can be based on a Bayesian Lasso approach, which can be very useful when there are many factors. In the context of a linear regression model, Tibshirani:1996 show that the Lasso coefficient estimates can be interpreted as the Bayesian posterior mode estimates under the independent Laplace prior distributions for the coefficients. Using this connection, Park:2008 develop a Gibbs sampler by assuming independent Laplace prior distributions for the coefficients. The sampler is based on the property that a Laplace distribution can be represented as a scale mixture of normal distribution, where the scale variable has an exponential distribution Andrews:1974. Let $\lambda_m(\mathbf{s}_i)$ be the $m$th element of $\boldsymbol{\lambda}(\mathbf{s}_i)$ for $m=1,\hdots,q$. Following Park:2008, we assume the following Laplace distribution for $\boldsymbol{\lambda}(\mathbf{s}_i)$:
where $\phi\geq0$ is the shrinkage parameter. In order to express this prior distribution in terms of a scale normal distribution, we use the following identity Andrews:1974:
In (ref), if we set $c=\phi$, $h=\lambda_m(\mathbf{s}_i)$, and $a=\tau^2_m$, where $\tau^2_m$ is an hyper-parameter for $\lambda_m(\mathbf{s}_i)$, we obtain
The right hand side of (ref) suggests that $\lambda_m(\mathbf{s}_i)|\tau^2_m\sim N(0,\tau^2_m)$ and $\tau^2_m|\phi^2\sim\text{Exp}(\phi^2/2)$, where $\text{Exp}(a)$ denotes the exponential distribution with rate $a$. Let $\boldsymbol{\tau}^2=(\tau^2_1,\hdots,\tau^2_q)^{'}$ and $\mathbf{D}_{\tau}=\text{Diag}(\tau^2_1,\hdots,\tau^2_q)$. Then, following Park:2008 and Casella:2010, we assume the following prior distributions:
where $\text{Gamma}(c,d)$ denotes the gamma distribution with shape parameter $c$ and the scale parameter $d$. This approach will require two new sampling steps for $\boldsymbol{\tau}^2$ and $\phi^2$, and the adjustment of the sampling step for $\boldsymbol{\lambda}(\mathbf{s}_i)$. We summarize this shrinkage MCMC approach in Algorithm (ref).
In this section, we consider a simulation exercise to elucidate the finite sample performance of the suggested Gibbs samplers with and without regularization. To this end, we consider the data-generating process in equations (ref) and (ref). We set $\rho=0.16$, $\gamma = 0.15$, $\delta = 0.2$ and $\beta = -2$. These parameter values are chosen to ensure that the data-generating process mimics the findings from the first empirical application in the next section.
We consider one exogenous variable drawn independently from the uniform distribution on the unit interval. We generate the error terms $\varepsilon_{t}(\boldsymbol{s}_i)$'s independently from the standard normal distribution. The number of factors $q$ is set to two, and they are drawn independently from the standard normal distribution along with their loadings. The number of spatial units is set to $49$, and they are randomly assigned to a regular rectangular lattice. We consider a queen contiguity spatial weight matrix so that $m_{ij}$ is set to 1 if the $j$-th entity is adjacent to or shares a border with the $i$-th entity. Otherwise, $m_{ij}$ is set to zero. We then row normalize the weights matrix. The number of time periods takes values from $\{50, 100\}$. The length of the Markov chain is $100000$ draws, and the first $20000$ draws are discarded to dissipate the effects of the initial draws.
We first estimate the dynamic spatiotemporal ARCH model with two factors using the Gibbs samplers in Algorithms 1 (standard MCMC) and 2 (shrinkage MCMC). We also consider the performance of the Gibbs samplers when the number of factors chosen by the researcher exceeds the true number of factors. To this end, we estimate the dynamic spatiotemporal ARCH model with three factors using the Gibbs samplers in Algorithms 1 and 2. Ideally, in this case, the Bayesian shrinkage estimator in Algorithm 2 shall shrink one of the factor loadings toward zero and shall perform better than the estimator in Algorithm 1.
Table (ref) presents the results from the simulation exercise. We present the estimated posterior median and the 95% credible interval for the parameters of the spatiotemporal ARCH model. We observe that both estimation algorithms report posterior median estimates that are close to the true values. The 95% credible interval estimates indicate evidence for the precision of the suggested estimation algorithms. For $\mathbf{h}^{*}$, we report the overall average along with the 95% credible interval using the median posterior draws. The true overall average value $\mathbf{h}^{*}$ is $-1.346$ when $T=50$ and $-1.354$ when $T=100$. The overall average estimates reported by both algorithms are close to the true averages. The DIC estimates indicate that when $T$ is short relative to $n$, the DIC may prefer the specification with a higher number of factors. However, when $T$ is large relative to $n$, the DIC chooses the model with the true number of factors.
To determine the adequacy of the length of the chains and their mixing properties, some exemplary trace plots are provided in Figure (ref). For the sake of brevity, we only present the results for Algorithm 1 when $q_{\max}=2$ and $T=100$. In these trace plots, the red solid lines correspond to the estimated posterior means. We observe that the Bayesian sampler performs satisfactorily and mixes well in all cases. For $\mathbf{h}^{*}$, we calculate the average of true values over $n$ and over $T$, respectively, and plot them against the average posterior means over $n$ and over $T$, respectively. In Figure (ref), we observe that the Bayesian sampler performs satisfactorily in terms of capturing the log volatility over spatial units as well as over time.
In this section, we use our model to estimate volatility in the US house price returns at the state level and the US stock market returns. In addition to volatility estimates, our approach allows us to estimate the spatial, temporal, and spatiotemporal effects of log-squared returns on the log volatility of housing and stock markets.
For the first application, we consider the housing market in the US, which provides a natural setting to apply a spatio-temporal model baltagi2014further, holly2010spatio, valentini2013modeling. Indeed, house prices are affected by spillovers from neighboring units, as houses share many common amenities. For example, kuethe2011regional, examining housing prices in the western US states from 1998 to 2007, finds that in most states, prices are influenced by spatial spillovers. Analyzing data from 1975 to 2011, brady2014spatial finds similar evidence across all US states. Similar findings are common for other housing markets outside the US, such as Germany otto2018spatiotemporal, the Netherlands van2011modelling, and the UK blatt2023changepoint, among others. Moreover, it is well known that house prices are influenced by common factors originating from economy-wide shocks, such as changes in monetary policy, energy prices, and technology. Previous studies have shown that models incorporating common factors are better suited for modeling and predicting house prices compared to standard models bork2018housing, yang2021common.
Spatial spillovers are also present in the volatility of house prices, a common measure of market risk. To account for this, tacspinar2021bayesian suggest using a stochastic volatility model with spatial interactions for modeling house price volatility, while dougan2023bayesian, otto2023general, otto2024multivariate propose spatiotemporal GARCH-type models, among others. While the presence of latent factors is crucial for explaining housing market volatility fairchild2015understanding, models that incorporate both spatial spillovers and latent factors in volatility modeling remain largely overlooked. We contribute to the literature on housing market volatility with our proposed specification, which accounts for the spatial, temporal, and spatiotemporal effects of log-squared returns on log volatility while also allowing for the influence of common latent factors.
For this application, we consider the quarterly house price indexes for the continental US states between 1975 and 2023. The data source is the Federal Housing Finance Agency website\footnote{Data can be retrieved from the following link: https://www.fhfa.gov/}. Our dataset includes $N=49$ spatial units and $T=193$ quarterly observations. We incorporate the average house price index for the US as an observable factor in the model. Figure (ref) shows the temporal evolution of the house price index for the US, where each boxplot represents the distribution of prices at a given point in time across all US states. We observe a generally positive trend with some local peaks, while the spread of the index across states also increases over time. Additionally, the index distribution appears to be right-skewed, with an increasing degree of asymmetry over time.
To measure spatial interaction in volatility dynamics, following previous studies, we consider a row-standardized contiguity spatial weight matrix $\mathbf{M}$. That is, all directly neighboring states are equally weighted, and only direct neighbors influence the price index in a given state. To select the number of latent factors in our specification, we use the DIC and Bayesian Lasso approaches presented in Section 3. We compute the DIC criterion for different values o $q$, using both Algorithms 1 and 2. The results are presented in Table (ref). Eventually, we choose the $q$ value that yields the lowest DIC criterion, specifically selecting the model with $q=7$ latent factors.
House price volatility can be driven by several economy-wide underlying factors of various types. For instance, we may encounter several macroeconomic factors influenced by global economic conditions, such as interest rates, employment, energy prices, wage growth, and technology. However, supply and demand also play a crucial role; factors like rising construction costs or population decline can either increase or depress prices. Investor activity is another important factor that can influence housing markets. For example, changes in investor sentiment—often triggered by speculative bubbles bursting—can lead to sharp price corrections, contributing to market volatility. Therefore, incorporating a larger number of factors allows us to account for the effects of many of these influences.
Using Algorithms (ref) and (ref), we estimate our specification by setting the length of the MCMC chains to $100000$, discarding the first $20000$ draws as burn-ins. The trace plots for the parameters $\beta$, $\gamma$, $\rho$, and $\delta$ show stable behavior, indicating convergence of the algorithm. The median posterior draws serve as reliable point estimates for the model parameters. Recall that $\beta$ is the parameter for the observable factor (representing market volatility), $\gamma$ denotes the temporal effect parameter, $\rho$ is the spatial effect parameter, and $\delta$ is the spatiotemporal effect parameter.
We present detailed estimation results based on Algorithms (ref) and (ref) in Table (ref). In this table, we provide the estimated posterior median for each parameter, $\boldsymbol{\lambda}^{'}_i \boldsymbol{f}_t$, and $\log h_t(\mathbf(s)_i)$ along with their 95% credible intervals. We include estimation results for $q$ ranging from $1$ to $8$ and highlight the selected model by the DIC criterion in gray. In this table, statistically significant coefficients are those whose confidence intervals do not include zero.
Focusing on the results obtained with $q=7$, highlighted in gray, we find that all parameter estimates are statistically significant, except for the estimate of $\beta$. The estimate of $\rho$ indicates that the house price volatility of a state is positively affected by the log-squared house price returns of neighboring states. Thus, this parameter represents the contemporaneous spillover effect, underscoring how current volatility from the housing markets in the neighboring states influences the housing sector in real time. The positive and statistically significant estimate of the spatiotemporal effect ($\delta$) shows that past shocks in house markets in the neighboring states continue to affect housing volatility. In both cases, we observe statistically significant spillover effects in the US housing market. The estimates of $\rho$ and $\delta$ are $0.16$ and $0.19$, respectively. Therefore, both instantaneous and lagged spillovers have similar effects in terms of magnitude, with the spatiotemporal effect being slightly larger than the instantaneous spatial effect. The estimated temporal effect ($\gamma$) is positive and statistically significant across all models, suggesting that volatility persists over time within the housing market, with current volatility being influenced by past market risks. This persistence is key to understanding the cyclical nature of housing market volatility, where a shock today continues to have repercussions in the future. Finally, we find that the volatility of the US house price index $(\beta)$, obtained as the average of the state-specific house price indexes, is not statistically significant. This is probably because the market effect is better captured by some of the $q$ latent factors.
Next, we analyze the estimated values of the log volatility in the US housing market. Figure (ref) shows the median posterior draws of the log volatility $\widehat{\log h_{t}(\mathbf{s}_i)}$, plotted across time for all states based on $q=7$ latent factors. In the figure, each grey line corresponds to one state, while the colored points show the medians across all states. We highlight three different time periods in Figure (ref), each characterized by a shock period. The first one corresponds to the early 1990s recession in the US, between 1989 and 1991; the second one is the dot-com bubble between 2000 and 2001, while the third one corresponds to the great financial crisis of 2007--2009 and the period just before the crisis occurred.
Figure (ref) shows that spikes of volatility are not always observed during economic crises. Between 1989 and 1990, the US economy was weakening due to a restrictive monetary policy, exacerbated by the 1990 oil shock. Despite these shocks, we do not observe a clear spike in the volatility of the housing market during this time period. The opposite applies, however, to the dot-com bubble in the early 2000s. As happened in the stock market, we also observed an increase in the volatility of housing prices during this time period. Finally, for the great financial crisis, we observe a significant spike in US house price volatility because volatility follows a decreasing pattern in the period immediately before the crisis, that is, between 2006 and early 2007. Interestingly, we find increasing volatility during COVID-19, that is, between 2020 and 2022, while we observe a strong decrease in US housing market volatility after the shock.
In Figures (ref)-(ref), we provide a detailed analysis of the regional differences in the estimated volatilities during the three considered periods. In these figures, the median posterior draws of volatilities based on the selected model with $q=7$ are shown through the US states maps. The color scale in these figures corresponds to the color scale in Figure (ref).
Considering the savings and loan crisis between 1989 and 1991 (see Figure (ref)), we find that states in the West are characterized by larger volatility compared to states in the East of the country. However, the states in the Northeast, such as New Hampshire, New York, and New Jersey, show relatively larger volatility. These results confirm the existence of positive spatial spillovers in the volatility dynamics. Figure (ref) shows the spatial volatility patterns during the dot-com bubble. In this case, the spatial differences become even more evident, as states in the West show larger volatility compared to those in the center and in the East, with the exception of northeastern states that are characterized by high volatility during the whole period.
Finally, we consider the volatility patterns during the great financial crisis. The results are shown in Figure (ref). In the quarters before the crisis, that is, during 2006, we find that all the states are characterized by relatively low volatility levels. From the first quarter of 2007, we observe that southwestern states such as Arizona, California, and Nevada show the highest volatility in the entire country. Starting from the last quarter of 2007, however, Florida and the northeastern states also begin to show high volatility levels. Arizona, California, and Nevada, in particular, exhibit a considerable increase in volatility. The first and last quarters of 2008, as well as the first quarter of 2009, show the highest levels of volatility for the entire US. After this point, we observe a sudden decrease in volatility for the remaining quarters in Figure (ref).
\FloatBarrier
Overall, the detailed analysis of these periods shows peculiar spatial patterns; that is, western and northeastern states are characterized by large volatility, while central states of the US have lower volatility levels. Moreover, our model allows us to obtain individual volatility or real-estate market risk levels for each state and time point.
In this section, we consider an application of the proposed model for analyzing financial volatility in a stock market network. We study the weekly time series of the stocks included in the Dow Jones Industrial Average over the last two years. Therefore, we consider $N=30$ stocks over a time period with length $T=100$. Figure (ref) shows the weekly log returns (top figure) and squared log returns (bottom figure) of the $N=30$ stocks included in the Dow Jones Industrial Average Index during the period 2022-2023. The weekly returns are summarized as boxplots, showing the cross-sectional variation of the log returns and squared log returns at each point in time, while the solid line represents the median returns across all $N$ stocks.
Figure (ref) highlights an important cross-sectional variability in terms of returns temporal patterns, with many outliers stocks at different points in time. For better modeling the volatility patterns of these stocks, it may be useful to consider the information of similarly outlying stocks. Moreover, consistent with previous studies alizadeh2002range,hansen2006consistent, Figure (ref) also shows that the log squared returns time series is a very noisy proxy of volatility, as it is subject to significant random fluctuations and may not consistently reflect the underlying volatility patterns. We are, therefore, interested in getting a suitable estimate of the volatility in this application.
Before presenting the estimation results, we compare our specification with commonly used GARCH models and their exponential extensions, such as log-ARCH models. Compared to the standard log-ARCH model, the proposed network log-ARCH approach provides two key features for volatility estimation. First, it accounts for contemporaneous and temporally lagged volatility spillovers through the parameters $\rho$ and $\delta$, which are well acknowledged in the financial literature diebold2014network,diebold2015trans. Second, it considers the effect of latent economy-wide factors that influence volatility across the entire stock market han2006asset,barigozzi2016generalized,calzolari2021latent. The inclusion of latent factors is particularly important when considering stocks traded in the same market index, as the effect of latent factors is stock-specific in the model. Since the model allows for the inclusion of observed factors, $\mathbf{x}_t(\mathbf{s}_i)$, we consider market volatility—the log-volatility of the Dow Jones Index—as an observed covariate in the model. As a result, our volatility estimation accounts for the spatial effect, the temporal effect, the spatiotemporal effect, the economy-wide latent factors effect, and the observed market volatility.
To construct the financial network matrix $\mathbf{M}$, we follow a consistent strand of literature onnela2004clustering,tumminello2010correlation,kumar2012correlation that adopts a correlation-based distance to specify the elements of the network matrix. More specifically, we first define the correlation-based distance as
where $\rho_{ij}$ is the Pearson correlation coefficient between two stocks $i$ and $j$. Given this transformation, the individual correlation coefficients are mapped from $[-1,1]$ to $[2,0]$, creating a symmetric distance matrix. We then define $m_{ij}=d^{-1}_{ij}$, so that a higher pairwise correlation coefficient indicates a stronger connection between the two stocks $i$ and $j$. Figure (ref) shows the network and spatial representation of the $N=30$ stocks included in the Dow Jones Industrial Average Index. In particular, highly correlated stocks appear close together in the network depicted in Figure (ref), while weakly correlated stocks are positioned further apart. All edges with a weight greater than the median edge weight are displayed and colored according to their values, with darker colors representing higher weights.
Figure (ref) illustrates the existence of core-periphery patterns in the Dow Jones Industrial Average index. For instance, CVX is positioned farthest from the main cluster of stocks and is correlated with other peripheral stocks like DOW and CAT. Another group of peripheral stocks includes MRK, UNH, and JNJ, while IBM stands out as the most central stock. This is not surprising, as IBM is one of the oldest stocks and one of the largest in terms of market capitalization. The darkest links can be found in the right part of the network, which also represents its central region. This suggests that the stocks in the core of the network are the most correlated, while those in the periphery are less correlated with the core stocks. These results are particularly significant for volatility modeling, as the information from peripheral stocks can enhance the modeling of their volatility, while core information can be utilized to model the volatility dynamics of stocks located in the core of the network.
To apply the network log-ARCH model with common factors, it is essential to choose the number of factors, $q$. We estimate various models using Algorithm 1, varying the number of factors, and select the $q$ value that minimizes the DIC criterion. The DIC results for different $q$ values are presented in Table (ref). We compare the results obtained using Algorithm (ref) (standard MCMC) with those from the Bayesian Lasso approach (Algorithm (ref), shrinkage MCMC), which automatically selects the relevant factors.
In Table (ref), we observe that the Bayesian Lasso approach provides a larger DIC compared to Algorithm (ref). Therefore, following the results in Table (ref), we choose the Network log-ARCH model with $q=7$ common factors estimated using Algorithm (ref). The selection of a relatively large number of latent factors is not surprising when analyzing volatility dynamics in the stock market. For instance, mittnik2015stock identified five distinct categories of volatility factors: equity market, interest rates and bond market, exchange rates, liquidity, and macroeconomic variables. To these categories, we could also consider additional ones, for example, those based on behavioral variables akin2024behavioral. The selection of $q=7$, therefore, appears consistent with the findings of previous studies.
The parameters $\beta$, $\gamma$, $\rho$, and $\delta$ correspond to the observable factor (market volatility), temporal effect, spatial effect (contemporaneous volatility network spillover), and spatiotemporal effect (lagged volatility spillover), respectively. The trace plots indicate stable behavior, suggesting that the MCMC chain produced by the algorithm converged. The estimated posterior medians for these parameters serve as our point estimates, and the results, along with their 95% credible intervals, are presented in Table (ref). The table displays the estimates obtained using Algorithm 1 with different numbers of factors $q$ and Bayesian LASSO. Statistically significant coefficients are those with the confidence interval that does not include zero.
In Table (ref), we present the estimation results for $q$ values ranging from $1$ to $8$, with a focus on the results for $q=7$, which is highlighted in gray as the best model. In this case, both the parameters associated with the observable market volatility factor ($\beta$) and the own-nodes temporal lag ($\gamma$) are not statistically significant. However, we find that the instantaneous volatility spillover effect parameter ($\rho$) and the lagged cross-nodes volatility spillover effect parameter ($\delta$) are positive and statistically significant. This is notable because only the network information, through spatial and spatiotemporal effects, remains important once we account for the presence of latent factors. This also suggests that the past behavior of the entire network, adjusted for contemporaneous effects, is more influential than the lagged log return of the individual node (own past). Furthermore, the insignificance of the market volatility proxy suggests that this information is captured by the latent factors estimated in the model. Additionally, both the instantaneous and lagged spillover effects are positive, indicating that an increase in the volatility of adjacent nodes in the network is associated with a rise in volatility. The magnitude of the instantaneous effect $(0.58)$ is also larger than that of the lagged effect $(0.26)$. Also, note that the estimates obtained from the Bayesian Lasso approach using Algorithm (ref) with $q_{\max}=8$ are similar to those shown in the gray highlighted row.
The results obtained with other models are qualitatively very close to those obtained using $q=7$ latent factors. However, we think it is worth highlighting that in the case of $q=1$ and $q=2$, we find the observable market factor to be statistically significant. This suggests that up to two factors are not enough to remove the effect of market volatility or, said differently, none of the estimated two factors contain information about market volatility. All the estimated models, starting from $q=3$, instead show no significant $\beta$ coefficient. Considering also the DIC values shown in Table (ref), we find that these two simpler models with $q=1$ and $q=2$ are those with the lowest performance.
Finally, we show the estimated value of volatilities and the resulting network structure, highlighting similarities in the financial network. Figure (ref) shows the median posterior draws of the log volatilities obtained with Algorithm 1, assuming $q=7$ factors. The volatilities are plotted across time (left) and space/network (right). For the time series plot, the estimated log volatilities are aggregated in the form of boxplots, with the median volatility displayed as a solid central line. The interquartile ranges of the boxplots highlight a large variability in the volatility patterns for the stocks included in the DJIA index. For some periods, outlier stocks with exceptionally high (or low) volatility can be identified.
\FloatBarrier
On the right of Figure (ref), we display the estimated volatility patterns of the considered US stocks. Stocks with similar volatility dynamics are placed in close proximity, while stocks with less correlated volatility dynamics are plotted further apart. The color of the labels is chosen according to the median log volatility across time, whereby the color scale of the left and right plots matches and can be used as a reference to relate both graphics. The larger the label size of the stock, the larger its interquartile range of the log volatilities, i.e., the risk/volatility profile is changing more over time compared to all other stocks. The figure highlights the existence of two groups of stocks: those with high and low volatility levels. In particular, BA, CRM, and INTC are the stocks with the highest volatility values and are those closer in the network. The same applies to AXP, DIS, and CVX, albeit with a lower level of volatility compared to the previous ones. On the other side of the network, we find many stocks with lower levels of variability, namely AMGN, WMT, and MRK. HON and VZ, the stocks with the lowest volatility, are neighbors and located in the center of the network structure.
In summary, Figure (ref) confirms the existence of positive volatility spillovers because stocks with large volatility are closer to those with similarly large volatility levels, and the same applies to stocks with lower variability. This evidence corroborates the findings of previous studies on the existence of relevant spillovers in the stock market and on the usefulness of using network methods for modeling financial volatility.
This paper introduced a dynamic spatiotemporal and network log-ARCH model that incorporates spatial, temporal, and spatiotemporal spillover effects with common factors. A key innovation of this model is the inclusion of volatility-specific time-varying latent factors that evolve over time, capturing complex volatility dynamics more effectively. Unlike conventional methods, this approach integrates both spatial dependencies—how volatility in one location affects another—and temporal evolution, providing a richer, more detailed understanding of volatility patterns across time and space. The model also has a network interpretation, linking different regions or assets through their volatility dependencies, further enriching the understanding of how shocks propagate across space and time.
We suggest two Gibbs samplers to estimate the parameters of the model. In addition to parameter estimates, our Bayesian estimation approach also allows us to estimate the volatility of the outcome variable. We suggest using the DIC criterion and a Bayesian Lasso approach to choose the number of factors that yield better predictive performance. We provide simulation results demonstrating that our proposed algorithms perform satisfactorily in recovering the true parameter values and volatility.
The proposed model is particularly well-suited for analyzing markets where spatial dependencies play a crucial role, such as the US housing market, but it can also be applied for modeling more general types of network processes. Two practical applications are proposed: one analyzing house price volatility and the other focusing on stock market data. In both cases, the model's ability to capture spatial and spatiotemporal dynamics offers a more comprehensive view of volatility transmission mechanisms, highlighting the model's broad applicability and its potential to enhance decision-making in financial and real estate markets.
In the first application, the US housing market provides a natural context for applying spatiotemporal models, given the established evidence of spatial spillovers in house prices. The second application extends this model to stock market data, where volatility spillovers are known to occur across different stocks and market indices. The empirical findings reveal that positive spillovers are present in both the housing and stock markets. However, network spillovers are much more pronounced in the stock market, where assets are more tightly interconnected, and volatility spreads more rapidly across the network. These findings highlight the model's ability to capture varying degrees of interconnectedness across different markets, making it a versatile tool for understanding complex volatility dynamics.
While the proposed approach is useful for volatility modeling, less is known about its forecasting performance. This is particularly interesting in both empirical applications we have considered in the paper. On one hand, in the case of house prices, forecasting volatility is an important indicator of economic stability, and cycles in housing market volatility can be used for predicting crises. On the other hand, in the case of financial markets, forecasting volatility is a useful tool for portfolio selection. Moreover, accurate volatility forecasts can also be used by policymakers to track conditions in the stock market and may help prevent crises. Therefore, a future study on the forecasting performance of our model, comparing its performance with other suggested models in the literature, can be explored. Moreover, future research may extend the proposed method to accommodate more complex factorial structures, such as the existence of a cluster-level structure.