A Dynamic Spatiotemporal Stochastic Volatility Model with an Application to Environmental Risks

This article introduces a dynamic spatiotemporal stochastic volatility (SV) model with explicit terms for the spatial, temporal, and spatiotemporal spillover effects. Moreover, the model includes time-invariant site-specific constant log-volatility terms. Thus, this formulation allows to distinguish between spatial and temporal interactions, while each location may have a different volatility level. We study the statistical properties of an outcome variable under this process and show that it introduces spatial dependence in the outcome variable. Further, we present a Bayesian estimation procedure based on the Markov Chain Monte Carlo (MCMC) approach using a suitable data transformation. After providing simulation evidence on the proposed Bayesian estimator's performance, we apply the model in a highly relevant field, namely environmental risk modeling. Even though there are only a few empirical studies on environmental risks, previous literature undoubtedly demonstrated the importance of climate variation studies. For example, for local air quality in Northern Italy in 2021, we show pronounced spatial and temporal spillovers and larger uncertainties/risks during the winter season compared to the summer season.

random process. In particular, for small-scale spatial units, the variance of the process is increased (known as Arbia's law of geography, Arbia and Espa, 1996). In addition to the instantaneous spatial correlations, we also have to account for the natural temporal correlations, which usually occur if we repeatedly observe a random process over time. The closer two observations are to each other in time, the more strongly they can correlate in general. In this paper, we introduce a new model for spatial, temporal, and spatiotemporal correlations in the log-volatilities allowing for additional random errors in the mean and volatility equation. Furthermore, we apply this model to environmental data and show for the first time how it can be used to analyze environmental risk factors such as air pollution.
There are generally two ways to account for spatial/cross-sectional correlations in spatial statistics. Firstly, it can be modeled in the covariance matrix of the process, where each entry is supposed to follow a certain (non-)parametric covariance function depending on the distance between their locations. This idea is typically known as the geostatistical approach (Cressie and Wikle, 2015;Zimmerman, 2019). The selection of a suitable parametric covariance function with certain properties such as stationarity, separability and full symmetry is one of the main modeling issues of this approach (Porcu et al., 2016;Huang et al., 2011). See Gneiting et al. (2007) for a review on the spatiotemporal covariance functions suggested in the literature. Secondly, the observations on an outcome variable can be explicitly correlated with the adjacent observations, where the adjacency is defined fairly generally by a spatial weights matrix. This second approach is referred to as spatial autoregression in spatial econometrics, where the spatial lags of variables are used to model spatial correlations. See LeSage and Pace (2009); Anselin (1988); Elhorst (2014); Lee (2004); Kelejian and Prucha (2010) on the specification and estimation issues in spatial econometrics. Both approaches can be equivalent under certain conditions (see e.g. Ver Hoef et al., 2018a, for simultaneous and conditionally autoregressive models).
In this paper, we consider the second approach. Our model consists of an outcome and logvolatility equation with separate independent error terms, whereby the log-volatility process introduces spatial dependence in the outcome variable. Specifically, the log-volatility equation allows for spatial, temporal, and spatiotemporal correlations, as well as time-invariant site-specific effects (unobserved heterogeneity). Also, assuming that the error terms in both equations have normal distributions, it is possible to show that the outcome variable has a leptokurtic symmetric distribution under our suggested model. To introduce a Bayesian estimation approach, we use a transformation approach such that the outcome equation becomes linear in the log-volatility terms. We use a Gaussian mixture distribution to approximate the distribution of the transformed error terms in the outcome equation. This approximation turns our model into a linear state-space model, where the log-volatility equation becomes the state equation. Following recent developments in the precision-based algorithms (Chan and Jeliazkov, 2009;Chan, 2017), we suggest a Gibbs sampler that consists of five steps for the estimation. We provide simulation evidence showing that the suggested sampler can perform satisfactorily.
Theoretically, our paper is related to the spatial econometric literature that addresses the pres-ence of cross-sectional correlations in higher moments of spatial data. This strand of the literature considers the spatial extensions of generalized autoregressive conditional heteroskedasticity (GARCH) and stochastic volatility (SV) models to account for volatility clustering patterns observed over space (Otto et al., 2018;Hølleland and Karlsen, 2020;Sato and Matsuda, 2021;Taşpınar et al., 2021;Robinson, 2009;Yan, 2007). Our model can be considered as the longitudinal data (panel data) extension of the cross-sectional spatial SV models suggested by Yan (2007);Taşpınar et al. (2021); Robinson (2009). While these studies allow for the presence of spatial dependence in the log-volatility equations, they do not include temporal, spatiotemporal and unobserved heterogeneity terms in the log-volatility equations. Our suggested process is also related to the separable and non-separable space-time filters considered in the spatial panel data models for modeling spatiotemporal interactions. For example, LeSage (2012, 2011) consider a separable space-time filter for the outcome variable, while Lee and Yu (2015) and Wang and Lee (2018) consider non-separable space-time filters for the disturbance terms of spatial panel data models.
In contrast to these studies, we consider a general space-time filter that also allows for unobserved heterogeneity, i.e., the site-specific effects, for the log-volatility of an outcome variable.
In an empirical application, we use our suggested model to assess environmental risk stemming from the variations in the log-volatility of air quality predictions. Spatial and spatiotemporal interactions in the local climate and environmental risks were addressed in comparably few studies in previous empirical research, even though it has been shown that an increased variation in environmental processes can be harmful (Tewksbury et al., 2008;Paaijmans et al., 2013;Vasseur et al., 2014;Iaco et al., 2012). Most previous studies focused on correlations in the (conditional) mean levels of ecological processes (Ver Hoef et al., 2018b;Wilby et al., 2009, e.g.). In our case, we model the log-volatility of fine dust concentrations of particles having a diameter less than 10µm in Lombardy, Northern Italy. Following the literature on the ecological processes (Ver Hoef et al., 2018b;Wilby et al., 2009, e.g.), we first model the variation in the conditional mean of our outcome variable through a conventional spatial panel data model that allows for the unobserved site and time heterogeneity. Our results from this initial model indicate that there is strong and moderate spatial correlations in the outcome variable in the model with only site fixed effects and the model with both site and time fixed effects, respectively. In the next step, we use the errors from this initial model as an outcome variable in our suggested specification and aim to model the variations in its log-volatility terms. The estimation results from our suggested Bayesian approach indicate that the spatial and temporal effects in the log-volatilities are moderate while the spatiotemporal effects appear to be of minor importance. We were also able to detect a noticeable variation in air pollution risk across the year and identify measurement stations that are associated with higher risks. These are mostly located in valleys in the Alpine regions.
The rest of this paper proceeds in the following way. In Section 2, we introduce our suggested model specification, including conditions ensuring the stability of the model and prior specifications.
In Section 3, we investigate the statistical properties of the suggested model. In Section 4, we provide the details on the posterior analysis of our model, and state an algorithm for the estimation.
Section 5 provides a simulation study on the performance of the suggested Gibbs sampler. Then, in Section 6, such stochastic volatility model is applied to environmental risks for the first time.
More precisely, our focus in on local air quality modeling. Finally, in Section 7, we provide our concluding remarks.

Model Specification
Suppose that we observe the spatiotemporal process across a constant set of n locations in a geographical domain at T equidistant time points. These locations can be measurement stations (i.e., marked point data), atmospheric and remotely-sensed data, or sets of municipalities, counties, states (i.e., areal data). Moreover, these locations does not have to seen in a strict geographical sense, but can also be vertices in a network. Let Y t = (y 1t , y 2t , . . . , y nt ) be the n × 1 vector of the outcome variable at time t for t = 1, 2, . . . , T . We assume the following data generating process (DGP) for Y t : . . , e 1 2 hnt ) is the n × n diagonal matrix containing stochastic volatility terms h it 's, which are specified subsequently, and V t = (v 1t , . . . , v nt ) is the n × 1 vector of disturbance terms. We assume that v it 's are i.i.d standard normal random variables.
Let h t = (h 1t , . . . , h nt ) be the n×1 vector of stochastic volatility at time t. We assume the following process for h t : (2.2) for t = 1, 2, . . . T , where µ = (µ 1 , . . . , µ n ) is the n × 1 vector of constant means, i.e., the timeinvariant site-specific effects, and U t = (u 1t , . . . , u nt ) is the n × 1 vector of i.i.d. disturbance terms such that u it ∼ N (0, σ 2 ) for all i and t. In (2.2), W is the n × n spatial weights matrix that has zero diagonal elements. This matrix specifies how volatility terms are related over space. In a network setting, this matrix is equivalently specified as an adjacency matrix. The scalar parameter ρ 1 captures contemporaneous spatial correlation, ρ 2 measures the temporal effect, i.e., the time dynamic effect, and ρ 3 represents the spatiotemporal effect, i.e., the spatial diffusion effect. Let L be the matrix of time-lag operator such that Lh t = h t−1 . Then, (2.2) can be written as where ((I n − ρ 1 W) − (ρ 2 L + ρ 3 WL)) is called the general space-time filter LeSage, 2011, 2012;Lee and Yu, 2015). Under the assumption that ρ 3 = −ρ 1 ρ 2 , this general filter is separable and decomposes into a product of the space filter (I n −ρ 1 W) and the time filter (I n −ρ 2 L).
In our analysis, we do not impose this restrictive assumption. Let S(ρ 1 ) = (I n − ρ 1 W). Then, under the assumption that S(ρ 1 ) is invertible, the reduced form of volatility equation is where A(ρ 2 , ρ 3 ) = (ρ 2 I n + ρ 3 W). When the cross-sectional dimension is fixed, the process for the log-volatility is stable if all eigenvalues of S −1 (ρ 1 )A(ρ 2 , ρ 3 ) lie inside the unit ball (Hamilton, 1994, Proposition 10.1). Let ϑ i (W) be the ith eigenvalue of W for i = 1, . . . , n. We assume that the parameter space of ρ 1 , ρ 2 and ρ 3 are chosen such that the following conditions hold: 1 The first condition is the sufficient condition for the invertibility of S(ρ 1 ) (Kelejian and Prucha, 2010, Lemma 1). By the spectral radius theorem (Horn and Johnson, 2012, Theorem 5.6.9), we can use any matrix norm to define relatively restrictive conditions that ensure the conditions in (2.5).
Finally, to complete the model in (2.1), we assume the following prior distributions for the posterior analysis: where Uniform(c 1 , c 2 ) denotes the uniform distribution over the interval (c 1 , c 2 ) and IG(a, b) denotes the inverse gamma distribution with the shape parameter a and the scale parameter b. The priors for ρ 1 , ρ 2 and ρ 3 are subject to the stability conditions sated in (2.5).

Statistical Properties
The outcome equation of our model can be written as Thus, the conditional variance of y it given h it is Var(y it |h it ) = e h it , indicating that the conditional variance is both time and space varying. Following the time series literature, we refer to h it as the log-volatility since h it = log (Var(y it |h it )). In order to determine the unconditional moments of y it , we need to determine the distribution of h = (h 1 , . . . , h T ) . Let ρ = (ρ 1 , ρ 2 , ρ 3 ) , and define the Then, we can express the log-volatility equation as We assume that the spatial dynamic process for h t has been operating for a long time so that we can express S(ρ 1 )(h 1 − µ) in the following way: Note that when ρ = 0, Ω reduces to σ 2 I nT , and thus the result in (3.7) becomes h|µ, σ 2 ∼ N (l T ⊗ µ, σ 2 I nT ). Consider the following partition of Ω: .
of Ω st for i, j = 1, 2, . . . , n. Let r ∈ N be a natural even number. Then, the even moments of y it can be expressed as our specification suggests that y it has a leptokurtic symmetric distribution. Next, we consider the covariance between y r it and y r js : This result indicates that our specification introduces spatial dependence in the outcome variable, since Cov y r it , y r js = 0 in general. Note that Cov y r it , y r js = 0 when ρ = 0 holds, because exp r 2 4 Ω ij,ts − 1 = 0.

Posterior Analysis
To introduce a Bayesian MCMC estimation approach, we first transform our model such that the resulting outcome equation is linear in h t . We then determine the conditional likelihood function of the transformed model by approximating to the distribution of transformed disturbance term with a Gaussian mixture distribution (Kim et al., 1998;Chib et al., 2002;Omori et al., 2007). The conditional likelihood function of the transformed model facilitates the sampling steps for h t and the auxiliary mixture component indicator defined subsequently. We also provide the conditional likelihood function of the original model, which we use to determine the sampling steps of other parameters in our model.
We square both sides of (3.1) and then take the logarithm to obtain The density of v * it is highly skewed with a long tail on the left and can be expressed as In order to convert (4.3) into a linear Gaussian state-space model, we approximate p(v * it ) with an m-component Gaussian mixture distribution: where φ(v * it |µ j , σ 2 j ) denotes the Gaussian density function with mean µ j and variance σ 2 j , p j is the probability of jth mixture component and m is the number of components. In particular, we use the ten-component Gaussian mixture distribution suggested by Omori et al. (2007) to approximate . We provide the parameter values of the ten-component Gaussian mixture distribution in Table 1. The parameters in this table are chosen by matching the first four moments of the ten component Gaussian mixture distribution with that of p(v * it ). This approach has two advantages. First, the Gaussian mixture distribution with the pre-determined parameter values in Table 1 provides a well enough approximation to p(v * it ) (Omori et al., 2007). Second, this approach does not pose any estimation difficulties since the mixture parameters in Table 1 are pre-determined.
We can equivalently write (4.4) in terms of an auxiliary discrete random variable z it ∈ where P(z it = j) = p j is the probability that z it takes the jth value. Let Z t = (z 1t , . . . , z nt ) , d t = (µ z 1t , . . . , µ znt ) and Σ t = Diag(σ 2 z 1t , . . . , σ 2 znt ). Then, from (4.5), we have V * t |Z t ∼ N (d t , Σ t ), which indicates that our model in (4.3) is now conditionally linear Gaussian given the component indicator variable. Thus, from (4.3), we have which facilitates the sampling steps for h t and Z t in our suggested Gibbs sampler given in Algorithm 1. The sampling steps for the remaining parameters requires the following conditional distribution: (4.7) We are now in a position to design a Gibbs sampler by using our results on ( Then, our suggested Gibbs sampler for generating draws from p(h, Z, µ, ρ, σ 2 |Y) consists of the steps given in Algorithm 1.

Remark 1. In
Step 1, the result in (4.10) indicates that the mixture components are conditionally independent given y * it . Thus, each component is a discrete random variable taking integer values in the interval [1, 10] with the conditional posterior probability P(z it = j|y * it ). The conditional posterior results in Steps 2, 3 and 4 are obtained from a standard Bayesian analysis as in a linear regression model. In the AM algorithm described in Step 5, the proposal distribution f g ρ|ρ (0) , . . . , ρ (g−1) has two parts. The first part is N ρ (g−1) , (0.1) 2 3 × I 3 , and is used when the number of iterations is less than or equal to g 0 . The second part consists of two normal distributions. The first component is specified as N ρ (g−1) , c(2.38) 2 3 × Cov ρ (0) , . . . , ρ (g−1) , where the covariance matrix is determined from the historical MCMC draws of ρ. The second component is N ρ (g−1) , (0.1) 2 candidate values generated from f g ρ|ρ (0) , . . . , ρ (g−1) are subject to the stability conditions given in (2.4). Finally, we adjust the tuning parameter c during the estimation to achieve an acceptance rate that falls between 40 percent and 60 percent.

Simulations
In this section, we provide simulation evidence to assess sampling properties of the suggested Bayesian algorithm. The data generating process follows (2.1)  For the spatial weights matrix W, we consider row-standardized rook and queen contiguity weights matrices. To this end, we first generate a vector containing a random permutation of the integers from 1 to n without repeating elements. Then, we reshape this vector into an k × m rectangular lattice, where m = n/k. In the case of rook contiguity, we set w ij = 1 if the jth observation is adjacent (left/right/above or below) to the ith observation on the lattice. In the case of queen contiguity, we set w ij = 1 if the jth observation is adjacent to, or shares a border with the ith observation. We set k = 7, and row-normalize all spatial weights matrices. For the prior distributions, we consider the following: σ 2 ∼ IG(3, 2) and µ ∼ N (0, 10I n ). The length of the Markov chain is 22000 draws, and the first 2000 draws are discarded to dissipate the effects of the initial values.
To determine the adequacy of the length the chains and their mixing properties, some exemplary trace plots are provided in Figure 1 and 3. For the sake of brevity, we only present the results for the queen contiguity case. In these trace plots, the red solid lines correspond to the estimated posterior means. We observe that the Bayesian estimator performs satisfactorily and seems to mix well in all cases. Note also that for both ρ and σ 2 , the 95% credible intervals contain the true values chosen in the experiments.
For the n × 1 vector µ and the nT × 1 vector h, we provide evidence on the performance of our Bayesian estimator in Figure 2 and 4. In these plots, the true values are represented with solid lines and the estimates are presented with dashed lines. In the first panel, we observe that the estimated posterior means for the components of µ are in general close to the true values. Here, the shaded region refers to the 95% credible interval. For 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. The second panel presents the case where the average is taken over n, and the last panel is the case where the average is taken over T . We observe that the Bayesian estimator performs satisfactorily in terms of capturing the log-volatility over cross-sections as well as over time.

Air Quality Modeling and Environmental Risks
In this section, we demonstrate the usage of the dynamic stochastic volatility model using an empirical example from environmental science. Like for applications in financial economics, the volatility of the process can be interpreted as risks (i.e., environmental risks in our case). Previous ecological studies mostly focused on changes in the mean behavior (cf. Wilby et al., 2009), but also an increased variation might be harmful, e.g., Vasseur et al. (2014) showed that an increased temperature variation (i.e., changes in the variance) poses a greater risk than global warming (i.e., changes in the mean). Similar results on this topic were found by Paaijmans et al. (2013) and  Tewksbury et al. (2008). Moreover, there is also a connection between climate changes and the financial market, such that climate variations might have an impact on the risk of financial markets (see also Giglio et al. 2021;Hong et al. 2020).
In this article, we analyze the variation of fine dust concentrations of particles having a diameter less than 10 µm, P M 10 , in Lombardy, Northern Italy. The region is surrounded by the Alps from the west, so that the wind circulation is reduced and Lombardy becomes one of the regions in Europe with the lowest air quality (see Fassò et al., 2022). In the following empirical analysis, we use the daily P M 10 concentrations from 1.1.2021 to 31.12.2021 from the official monitoring stations of the regional environmental authority, ARPA Lombardia (Maranzano, 2022). The data are open-source provided by the Agrimonia project (Fassò et al., 2022a,b). In total, there are n = 103 measurement stations and T = 365 daily observations. To provide a first overview of the dataset, we depicted the median concentrations in µg/m 3 across space and time as a time-series plot ( Figure 5, left) and displayed on map ( Figure 5, right), respectively.
For our analysis, we first estimated a spatial panel model to describe the mean variations. The (6.1) The outcome variable C t is the n-dimensional vector of P M 10 concentrations at time point t, W is spatial weights matrix which we specify below in more detail, X t is an n × p matrix of exogenous regressors, and E t denotes the model error terms. Spatial interactions are included via the spatial autoregressive term ψWC t with an unknown autoregressive parameter ψ. Moreover, spatial fixed effects s are present for each station (even if non-significant), because there are different types of stations included in the data set, e.g., urban traffic stations located at major roads in the cities, or rural background stations located in the Alps. The spatial fixed effects describe the station-specific P M 10 concentrations and thereby serve as the model intercept. We refer to this model as Model A. Moreover, we consider the same model with spatial and temporal fixed effects as an alternative model B, that is, The temporal fixed effects a t 1 n will remove any additional (station invariant) seasonal variation across time, which is not explained by the exogenous regressors. More precisely, we included p = 4 covariates that are known to influence PM concentrations, namely the maximum height of the planetary boundary layer (PBL), the relative humidity, the air temperature, and the pressure level.
All covariates are available as daily observations for each measurement station. Moreover, they were standardized to compare the size of the effects. The spatial weights matrix has been chosen as row-standardized binary contiguity matrix, where all locations within 15 miles are considered as adjacent stations. On average, each stations has 5.2621 neighbors leading to a sparsity level of W of 94.95 %. We estimated all parameters {s, β, γ, ψ, σ 2 E } using the maximum likelihood approach implemented in the spatial econometrics MATLAB toolbox (cf. Bivand and Piras 2015;LeSage 1999). The resulting estimates including their asymptotic 95% confidence intervals are summarized in Table 2.
The estimated parameters are in line with our expectations and we see a good general fit of the model (R 2 = 0.7829 and 0.8049 with only spatial and spatiotemporal fixed effects, respectively).
However, for this paper, our focus is not on the mean variations but on the model errors E t , i.e., the model uncertainties. Hence, we do not go into further detail on the interpretation of the mean  for each measurement station, we see that the error variance is varying across space and time (see Figure 6). There are time periods of increased variations, so-called volatility clusters, e.g. at the end of the year. Moreover, we observe a similar clustering across space. Measurement stations with a higher volatility are located in close proximity, e.g., North-East of Milan or around Brescia. This provides motivation for estimating a dynamic stochastic volatility model for E t in the second step, that is, as defined in (2.1) and (2.2). Table 3 reports the median of the posterior draws for each parameter including the corresponding 95% credible intervals (500,000 posterior draws, and 100,000 burn-in draws). The space-time interaction parameters ρ 1 , ρ 2 , and ρ 3 represent the uncertainty/risk spillovers. In general, we ob-

Conclusion
We have introduced a novel spatiotemporal statistical model for stochastic volatilities, which are spatially, temporally, and spatiotemporally correlated and have different levels for each location reflecting the heterogeneity of the process across space. The model includes two different error terms for the mean and the log-volatility equations. To estimate the parameters, we suggested a  Bayesian MCMC approach. To this end, we applied a log-square transformation to transform our non-linear state-space model into a linear one. Further, we used a Gaussian mixture distribution to approximate the distribution of the transformed error terms in the space equation, transforming our model into a linear Gaussian state-space model. We analyzed the estimation performance for different parameter settings and spatial interactions in a simulation study, and showed that the suggested Bayesian sampler performs satisfactorily.
Moreover, we applied the dynamic spatiotemporal stochastic volatility model in a completely new empirical framework, namely, in the field of modeling environmental and climate risks. First, statistical modeling of the volatility process of climate variables have not been done extensively yet, even though there is a large scientific consensus that an increased variability of environmental processes, e.g., the temperature variability, is harmful for the environment. Second, stochastic volatility models were predominantly applied to financial data, because of the straightforward interpretation of the log conditional volatilities as the (return) risk of financial assets. We have transferred to this idea to environmental risks and showed how the volatility of P M 10 predictions are correlated across space and time. In particular, measurement stations with an increased model uncertainty could be identified in this way. In addition, we showed in our application that there are significant temporal and spatial spillovers in these environmental risks. Also, the temporally lagged spatial spillovers, i.e., the spatiotemporal correlations, appear to play a minor role.
Both from a theoretical and applied perspective, spatiotemporal stochastic volatility models and