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.


Introduction
When analyzing geo-referenced data, statistical models have to account for instantaneous spatial correlations due to the geographical proximity between the observations.This is commonly known as Tobler's first law of geography "everything is related to everything else, but near things are more related than distant things" (Tobler, 1970).This observation was already noted by Ronald A. Fisher in 1935 as follows, "the widely verified fact that patches in close proximity are commonly more alike, as judged by the yield of crops, than those which are further apart" (Fisher, 1935).
Even though the similarity is typically considered to be in the (conditional) mean level at each location, there might also be spatial correlations in the (conditional) variance or variation of the 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, Parent andLeSage (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 : (2.1) for t = 1, 2, . . .T , where H ) 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 (( is called the general space-time filter (Parent andLeSage, 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: 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).
Let • be any matrix norm.Then, the sufficient conditions for (2.5) are (i) ρ 1 W < 1 and (ii) where last equality follows since we assume that ρ 1 W < 1.If we choose the matrix row sum norm • ∞ and assume that W is row normalized, then (2.6) reduces to (|ρ Thus, a further restrictive sufficient condition for the stability of (2.4) is will impose these restrictions during the sampling steps for ρ 1 , ρ 2 and ρ 3 in our suggested Gibbs sampler.
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 for all i = 1, . . ., n and t = 1, . . ., T .
(3.1) 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 nT × nT matrix J(ρ) as 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: where where 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 Ω: where each Ω st for s, t = 1, 2 . . ., T is an n × n sub-matrix of Ω.Let Ω ij,st be the (i, j)th element 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 :

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 where y * it = log y 2 it and v * it = log v 2 it .The density of v * it is highly skewed with a long tail on the left and can be expressed as Then, in vector form, we have 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 p(v * it ).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 ∈ , and P(z it = j) = p j , j = 1, 2, . . ., m, (4.5) where P(z it = j) = p j is the probability that z it takes the jth value.Let Z t = (z 1t , . . ., z nt ) , 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: 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 3 × I 3 .The 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.
Remark 2. In the sampling step for ρ, P(ρ (g−1) , ρ) is calculated at each pass of the sampler, and therefore, p(h|ρ, µ, σ 2 ) is evaluated twice at each pass of the sampler.In other words,

Simulations
In this section, we provide simulation evidence to assess sampling properties of the suggested Bayesian algorithm.The data generating process follows (2.1) and (2.2).More specifically, the elements of V t and U t are drawn independently from the standard normal distribution for t = 1, 2, . . ., T .Therefore, the value of σ 2 is set to 0.25 in all experiments.The elements of µ are drawn independently from the normal distribution with mean 3.3 and standard deviation 0.35.To initialize the process, we use (3.7), and the series expression for K(ρ) is truncated at 15.We consider two sets of values for ρ, {(0.6, 0.35, −0.025), (0.3, 0.65, −0.025)}.These parameter values are chosen to ensure that the data generating process mimics the findings from our empirical application in the next section.The number of spatial units n is set to 98, and the number of time periods T is fixed at 50.
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.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 Figure 2: Estimates of µ and h (dashed lines) and their corresponding true data-generating values (solid lines) for ρ = (0.6, 0.35, −0.025) and W is a queen contiguity matrix 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.2021to 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 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 observe moderate instantaneous spatial interactions (i.e., ρ 1 = 0.5943 for Model A and ρ 1 = 0.3037 for Model B) and temporal autoregressive interactions (i.e., ρ 2 = 0.3781 for Model A and ρ 2 = 0.6699 for Model B), while spatiotemporal effects appear to be of minor importance (i.e., ρ 3 = −0.0278and ρ 3 = −0.0236for Models A and B, respectively).The posterior mean estimates for both the spatial and temporal autoregressive parameters are positive.That is, if the environmental risk (logvolatility h it ) is high it is likely to spillover to the neighboring regions and to future time points.For this reason, spatial and temporal volatility clusters are formed.The negative signs of the posterior mean estimate for the spatiotemporal lag are opposing this behavior, but they are close to zero.
Moreover, we observe that the spatial spillovers are dominating the temporal ones for Model A without temporal fixed effects.That is, the missing temporal variation of mean model could be Moreover, we analyzed the conditional volatilities for each station across the time horizon.For Model B, the results are displayed in Figure 7 as a time-series plot (left) and on a map (right).More precisely, we computed the posterior medians of h it for all locations and time points leading to an estimate of the n × T matrix of the logarithmic conditional volatilities.Then, for visualization, we computed the median and the 5% and 95% quantiles of these estimates across all spatial locations (see time series plot) and across all time points (see map).From the time series plot, we see that there is a clear annual pattern -the model uncertainties are lower in summer and highest in winter.These varying conditional volatility levels can be interpreted as environmental risks, which are the highest the winter season.In the same manner, measurement stations with a greater uncertainty/risk can be identified from the map.They are mostly located in valleys in the mountain areas.

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 environmental risk modeling are important topics for future research.The current model does not allow for temporally varying constant terms in the conditional volatilities, which are constant across space.Regarding the latter case, other environmental processes, like temperature variations, soil droughts, or atmospheric ozone concentration and optical depths are important processes, where a deep understanding of the spatiotemporal interactions in the variabilities is essential, both for obtaining accurate prediction intervals and the planning of interventions (and their impact on the environment).
ts − 1 .(3.10)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.
.7) We are now in a position to design a Gibbs sampler by using our results on (i) the mixture component indicators in (4.5), (ii) the conditional likelihood function of transformed model in (4.6), (iii) the conditional likelihood function of Y t in (4.7) and (iv) the distribution of h in (3.7).Let Y = (Y 1 , . . ., Y T ) and Z = (Z 1 , . . ., Z T ) .The joint posterior distribution p(h, Z, µ, ρ, σ 2 |Y) can be expressed as

Figure 4 :
Figure 4: Estimates of µ and h (dashed lines) and their corresponding true data-generating values (solid lines) for ρ = (0.30, 0.65, −0.025) and W is a queen contiguity matrix

Figure 5 :
Figure 5: Median P M 10 concentrations (µg/m 3 ) and the 5% and 95% quantiles across space displayed as time series (left) and median concentrations across time shown on a map (right).

Figure 6 :
Figure 6: Standard deviation of residuals E t

Figure 7 :
Figure 7: Median posterior draws of h t (Model B) and the 5% and 95% quantiles across space displayed as time series (left) and median posterior draws of h t across time shown on a map (right).

Table 2 :
Estimated parameters of the mean model

Table 3 :
Estimated parameters of the dynamic stochastic volatility model, medians and 95% posterior credible intervals.