Quantile-based spatiotemporal risk assessment of exceedances

Stochastic Environmental Research and Risk Assessment https://doi.org/10.1007/s00477-018-1562-9 (0123456789().,-volV)(0123456789().,-volV) ORIGINAL PAPER Quantile-based spatiotemporal risk assessment of exceedances J. L. Romero1 • A. E. Madrid2 • J. M. Angulo1 Ó Springer-Verlag GmbH Germany, part of Springer Nature 2018 Abstract Structural characteristics of random field excursion sets defined by threshold exceedances provide meaningful indicators for the description of extremal behaviour in the spatiotemporal dynamics of environmental systems, and for risk assessment. In this paper a conditional approach for analysis at global and regional scales is introduced, performed by implementation of risk measures under proper model-based integration of available knowledge. Specifically, quantile-based measures, such as Value-at-Risk and Average Value-at-Risk, are applied based on the empirical distributions derived from conditional simulation for different threshold exceedance indicators, allowing the construction of meaningful dynamic risk maps. Significant aspects of the application of this methodology, regarding the nature and the properties (e.g. local variability, dependence range, marginal distributions) of the underlying random field, as well as in relation to the increasing value of the reference threshold, are discussed and illustrated based on simulation under a variety of scenarios. Keywords Conditional simulation
Quantile-based risk measures
Space–time random fields
Threshold exceedance indicators 1 Introduction The dynamics of a wide variety of physical phenomena (e.g. in Geophysics, Environmental Sciences, etc.) can be appropriately represented in terms of random field models from which different probabilistic and statistical aspects can be addressed (see, for example, Adler 1981; Christakos 1992, 2000; Piterbarg 1996; Vanmarcke 2010; Yaglom 1987a, b). In particular, analysis of extremal behaviour constitutes one of the main objectives in many applications oriented to risk assessment. A common approach in this & J. M. Angulo [email protected] J. L. Romero [email protected] A. E. Madrid [email protected] 1 Department of Statistics and Operations Research, University of Granada, Campus Fuente Nueva s/n, 18071 Granada, Spain 2 Department of Sciences and Informatics, University Centre of Defence at the Spanish Air Force Academy, MDE-UPCT, C/ Coronel López Peña s/n, 30720 Santiago de la Ribera, Murcia, Spain context consists in describing structural aspects of excursion sets defined by threshold exceedances (hereafter, ‘excursion sets’), that is the set of locations within the domain where the values of a given random field realization reach or are above a certain reference threshold (see, for example, Adler 2008; Adler and Taylor 2007; Adler et al. 2010, 2013; Azäis and Wschebor 2009). In fact, different functionals related to geometrical characteristics of excursion sets, such as the exceedance area, the excess volume or the number of connected components, among others, provide suitable random indicators useful for quantification of risk (see, for example, Adler and Taylor 2007; Angulo and Madrid 2010, 2014; Chiu et al. 2013; Christakos and Hristopulos 1996, 1997; Madrid et al. 2012, 2016; Yang and Christakos 2015). In the last two decades, the development of a wellfounded theory of measures of risk has arisen as a new scientific discipline, mainly motivated and also impulsed by areas of applications such as Finance and Insurance, although with an increasingly broadening interest, due to its potential applicability, in many other areas of knowledge such as Environmental and Health Sciences, Geophysics, Engineering, etc. (see, for example, Föllmer and Schied 2016, Klüppelberg et al. 2014). Among the variety of risk measure families introduced in this framework, 123 Stochastic Environmental Research and Risk Assessment quantile-based risk measures such as Value-at-Risk (VaR) and Average Value-at-Risk (AVaR) have received special attention because of their direct interpretation and easy computational implementation, besides the compliance of certain axioms. There is a very vast literature on applications of risk measures under this approach in diverse fields; in the geostatistical context, see, for example, Bernardi et al. (2018), Li et al. (2018). This theory quickly advances in different challenging directions from the complexity perspective, including systemic risk measures (e.g. Kleinow et al. 2017; Souza et al. 2016), risk transfer optimization (e.g. Filipović and Kupper 2008; Haier et al. 2016), etc. In this paper, a general and flexible methodology for spatial and spatiotemporal risk assessment is proposed based on a conditional approach. More specifically, for a given random field model and available observations, empirical distributions of different threshold exceedance indicators are obtained from conditional simulation, as a basis for evaluation of measures of risk. In particular, the study is focused on the global and local assessment of exceedance area and excess volume indicators, from which different risk maps are derived. For these so-called firstorder indicators, the compound cumulative distribution function (e.g. Craigmile et al. 2005; French and Sain 2013; Lahiri et al. 1999; Wright et al. 2003; Zhang et al. 2008) plays a key role, both formally and regarding the practical threshold specification (a formal definition is given in Sect. 2.3). The paper is structured as follows. Section 2 provides a basic introduction to essential concepts related to quantilebased risk measures and elements related to the threshold exceedance indicators considered. Section 3 describes in detail the methodological aspects, first formally and then from a computational point of view. In Sect. 4, based on simulations under different spatial (Cauchy, Student and Fisher–Snedecor random fields) and spatiotemporal (blurgenerated autoregressive model) scenarios is performed, different aspects of particular relevance are illustrated. Finally, in Sect. 5 a summary is given with further reference to some directions of continuing research in this context. 2 Preliminary elements This section introduces technical notation, definitions and some preliminary relation between the elements involved in the methodology proposed. First, some aspects related to risk measures are considered. 2.1 Risk measures Definition 1 Given a probability space ðX; F ; PÞ, and MB ðX; F ; PÞ, the space of real-valued Borel-measurable functions defined on ðX; F ; PÞ, a risk measure is a functional  q : MðX; F ; PÞ ! R satisfying certain suitable conditions, where MðX; F ; PÞ is  denotes the a linear subspace of MB ðX; F ; PÞ, and R compactification of R given by R [ f 1; þ1g. Often, MðX; F ; PÞ is taken to be an Lp ðX; F ; PÞ space, with p 2 ½1; 1Š. In this work, we assume that the argument random variables represent loss in some sense. As mentioned in Sect. 1, mainly based on financial and actuarial problems, some fundamental classes such as coherent and convex risk measures (see, for example, Föllmer and Schied 2016) are attributed significance, each one with some of the properties or axioms listed below: (i) (ii) (iii) (iv) (v) Translation invariance: qðX þ aÞ ¼ qðXÞ þ a, for a 2 R. Monotonicity: X  Y P a:s: ) qðXÞ  qðYÞ. Positive homogeneity: qðaXÞ ¼ aqðXÞ, for a [ 0. Sub-additivity: qðX þ YÞ  qðXÞ þ qðYÞ. Convexity: qðaX þ ð1 aÞYÞ  aqðXÞ þ ð1 aÞqðYÞ, for a 2 ½0; 1Š. (A detailed analysis of axioms (i)–(v) in a financial and actuarial framework is developed in Föllmer and Schied 2016). A risk measure is said to be coherent if it satisfies (i), (ii), (iii) and (iv), whereas is said to be convex if it satisfies (i), (ii) and (v). Every coherent risk measure is a convex risk measure, but not conversely. Coherent and convex risk measures are important, in particular, because they can be expressed in terms of mathematical expectations by means of dual representation theorems. (From a conceptual point of view, and in the context of financial risk assessment, coherency was formally introduced and justified as a desirable set of properties in the seminal paper by Artzner et al. 1999, while convexity was proposed by Föllmer and Schied 2002, and parallelly by Frittelli and Gianin 2002, for a less restrictive instrument regarding risk diversification.) Risk measures formulated in terms of distribution moments and/or quantiles are usually considered because of their easy interpretation in terms of the probability distribution. Here, in particular, the definition of the wellknown quantile-based risk measures Value-at-Risk (VaR) and Average Value-at-Risk (AVaR) is introduced. Definition 2 Let X be a loss random variable, with cumulative distribution function (CDF) FX . Let a 2 ð0; 1Š. The Value-at-Risk at confidence level 1 a of X, 123 Stochastic Environmental Research and Risk Assessment VaR1 a ðXÞ, is defined as the lower ð1 distribution of X: VaR1 a ðXÞ ¼ inffx 2 R : FX ðxÞ  1 aÞ quantile of the ag: It follows that if X is a continuous random variable, then VaR1 a ðXÞ ¼ FX 1 ð1 aÞ ; where ‘ ’ stands for the ‘lim inf’ of the set. Definition 3 Under the previous assumptions on X and a, the Average Value-at-Risk at confidence level 1 a of X, AVaR1 a ðXÞ, is defined as Z 1 1 VaRp ðXÞdp: AVaR1 a ðXÞ ¼ a 1 a In the case where X has a continuous probability distribution, AVaR1 a ðXÞ is a conditional expectation: AVaR1 a ðXÞ ¼ E½XjX [ VaR1 a ðXފ: Hence, in this case, while VaR1 a ðXÞ represents the minimal loss that will occur in the 100a% worst scenarios, AVaR1 a ðXÞ represents the average loss that will occur in the 100a% worst scenarios. In general, VaR1 a ð
Þ is not a coherent risk measure, because of lack of subadditivity except for special cases, whereas AVaR1 a ð
Þ is always a coherent risk measure. In fact, VaR does not take into account the tail behaviour above the specified confidence level quantile, whilst AVaR does it. A generalized formulation of quantile-based risk measures is given by the so-called spectral risk measures, defined as Z 1 M/ ðXÞ ¼ VaRp ðXÞ/ðpÞdp; 0 where / is a weight function such that /ð
Þ  0 and R1 0 /ðpÞdp ¼ 1. M/ ð
Þ is a coherent risk measure for monotonically non-decreasing weight functions (e.g. AVaR1 a ð
Þ). 2.2 Random field threshold exceedances Some aspects related to notation and definition of random field threshold exceedance indicators considered for the methodology proposed are introduced below. Hereafter, X denotes a spatial random field (RF), i.e. X ¼ fXðsÞ 2 MB ðX; F ; PÞ : s 2 S  Rd g; D  S is a bounded subdomain with non-null Lebesgue measure, kðDÞ [ 0: For simplicity, all the random variables X(s) are assumed to have a continuous distribution. Let u 2 R be a pre-fixed threshold level. Definition 4 Au ðX; DÞ ¼ fs 2 D : XðsÞ  ug is the excursion set of X in D over level u 2 R. While for d ¼ 1 this is the disjoint union of closed sets some of which may be degenerate, for d [ 1 the excursion set involves a higher geometrical complexity. In particular, aspects related to the fragmented structure of the set can be analyzed in terms of its connected components. A general treatment of the geometry of random field excursion sets based on intrinsic volumes is approached, for instance, by Adler and Taylor (2007). X  Definition 5 Iu ðsÞ ¼ 1fXðsÞug : s 2 S is the u-exceedance indicator RF. Definition 6 fXu ðsÞ ¼ maxf0; XðsÞ u-excess RF. ug : s 2 Sg is the As mentioned before, risk measures are evaluated on random variables representing loss in some sense. In this work, indicators related to exceedance areas and excess volumes are considered (see Table 1 for notation), under the assumption that suitable regularity and measurability conditions hold for their existence, as well as for the existence of their expectations (see, for instance, Adler 1981; Adler and Taylor 2007; Azäis and Wschebor 2009, etc.) For a given threshold u and a subdomain D, kðAu ðX; DÞÞ expresses the total area with the values of the random magnitude represented by the RF lying above the threshold u; kD ðAu ðX; DÞÞ represents the exceedance area ratio with respect to the subdomain D; E½kðAu ðX; DÞފ and E½kD ðAu ðX; DÞފ denote their respective expected values. Further, VðAu ðX; DÞÞ expresses the total excess volume of the RF above the threshold u; V D ðAu ðX; DÞÞ represents the excess volume ratio with respect to the subdomain D; as before, E½VðAu ðX; DÞފ and E½V D ðAu ðX; DÞފ denote their corresponding expected values. In Sect. 3, some other indicators of practical interest are analyzed, such as the volume ratio with respect to the area of the excursion set, V Au ðAu ðX; DÞÞ, and the average volume ratio per (non-degenerate) connected component of the excursion set, Vcc ðAu ðX; DÞÞ (this indicator is only considered when the D number of non-degenerate connected components, Ncc , is finite or, otherwise, for connected components with a certain minimum area size). It should be noted that the specific indicators here considered are invariant with respect to changes in the values of the random field variables corresponding to subsets of locations with zero Lebesgue measure. In particular, the possibility of degenerate excursion sets consisting of only isolated points is not excluded; however, in such cases the exceedance area, as well as the excess volume, will be 0. This can be interpreted in the sense that, under such indicators, risk is evaluated on the basis of non-degenerate occurrences of exceedances with respect to the dimension of the domain. (Nevertheless, as an aspect of interest, for regularity conditions under which 123 Stochastic Environmental Research and Risk Assessment Table 1 Some first-order random indicators of RF threshold exceedances Notation Definition kðAu ðX; DÞÞ R Au ðX;DÞ R [u0  u Au0 ðX;DÞ Indicator of exceedance area Absolute D Relative (w.r.t. D) ds kðAu ðX;DÞÞ kðDÞ k ðAu ðX; DÞÞ Indicator of excess volume Absolute VðAu ðX; DÞÞ Relative w.r.t. D V D ðAu ðX; DÞÞ Relative w.r.t. Au V Au ðAu ðX; DÞÞ D Relative per connected component Ci (average, for i ¼ 1; . . .; Ncc ) Vcc ðAu ðX; DÞÞ the probability of having degenerate tangence points to a given threshold becomes null, the reader is referred to multiparameter versions of Bulinskaya’s theorem in Adler and Taylor 2007; Azäis and Wschebor 2009, among other sources). In general, the distributions of these indicators are unknown, and analytical expressions or numerical approximations for the considered risk measures are unfeasible. This supports, as mentioned before, the importance of developing suitable methodologies that allow to compute these measures under general conditions. 2.3 Some preliminary results dsdu0 VðAu ðX;DÞÞ kðDÞ VðAu ðX;DÞÞ kðAu ðX;DÞÞ 1 D Ncc PNccD nVðAu ðX;Ci ÞÞo i¼1 kðAu ðX;Ci ÞÞ Proposition 1 Given a spatial RF X on S 2 Rd satisfying suitable regularity and measurability conditions, a bounded subdomain D  S; and a fixed threshold u 2 R; the following relations hold: R (i) kðAu ðX; DÞÞ ¼ D IuX ðsÞds (ii) E½kðAu ðX; DÞފ ¼ kðDÞSD X ðuÞ (iii) E½kD ðAu ðX; DÞފ ¼ SD ðuÞ R þ1 X R (iv) V ðAu ðX; DÞÞ ¼ u kðAu0 ðX; DÞÞdu0 ¼ D Xu ðsÞds R þ1 (v) E½VðAu ðX; DÞފ ¼ kðDÞ u SDX ðu0 Þdu0 ¼ kðDÞE½XuD Š R þ1 0 0 D (vi) E½V D ðAu ðX; DÞފ ¼ u SD X ðu Þdu ¼ E½Xu Š where XuD ¼ maxf0; X D ug. Proof In this subsection, some preliminary formalizations and results concerning the indicators defined in Sect. 2.2 and their expected values are introduced. Let FXD ðxÞ be the compound cumulative distribution function (CDF) defined as Z D FXðsÞ ðxÞkD ðdsÞ; FX ðxÞ ¼ ð1Þ (i) (ii) D ds where FXðsÞ denotes the CDF of X(s), and kD ðdsÞ :¼ kðDÞ represents the normalized Lebesgue measure on D. Let X D be the corresponding compound r.v., whose values are obtained by first randomly (i.e. according to the normalized Lebesgue measure) selecting a location s within the domain D and then assigning the corresponding X(s) realization. Let SD FXD ð
Þ, i.e. the decumulative distribution X ð
Þ ¼ 1 function (DDF) of X D . In what follows, suitable regularity and measurability conditions will be referred to hold for random field X in order to guarantee that the indicators considered are well defined, as well as for the existence of their expectations (see, for instance, Adler 1981, sec. 3.2; Adler and Taylor 2007, sec. 6.2; Azäis and Wschebor 2009, sec. 7.1; Vanmarcke 2010, sec. 3.6). 123 (iii) (iv) (v) Since IuX ðsÞ ¼ 1 if s 2 Au ðX; DÞ, and IuX ðsÞ ¼ 0 if s 62 Au ðX; DÞ, it follows that kðAu ðX; DÞÞ ¼ R R X Au ðX;DÞ ds ¼ D Iu ðsÞds. R By Fubini’s theorem, E½kðAu ðX; DÞފ ¼ D P  R R  ½XðsÞ  uŠds ¼ D SXðsÞ ds ¼ D 1 FXðsÞ ðuÞ ds ¼ R kðDÞ D FXðsÞ ðuÞds, with SXðsÞ being the DDF of X(s). Hence, by Eq. (1), E½kðAu ðX; DÞފ ¼ kðDÞð1 FXD ðuÞÞ ¼ kðDÞSD X ðuÞ.   kðAu ðX; DÞÞ 1 D E½k ðAu ðX; DÞފ ¼ E E½k ¼ kðDÞ kðDÞ ðAu ðX; DÞފ ¼ SD ðuÞ. RX VðAu ðX; DÞÞ ¼ D½u;þ1Þ IuX0 ðsÞdsdu0 : Hence, by R þ1 Fubini’s theorem, VðAu ðX; DÞÞ ¼ u kðAu0 ðX; R DÞÞ du0 ¼ D Xu ðsÞds. From the proof of (iv), E½VðAu ðX; DÞފ ¼ R 0 0 D½u;þ1Þ P½XðsÞ  u Šdsdu . Hence, correspondR þ1 0 0 ingly, E½VðAu ðX; DÞފ ¼ kðDÞ u SD X ðu Þdu and R þ1 0 E½VðAu ðX; DÞފ ¼ kðDÞ u ðu uÞFXD ðu0 Þdu0 ¼ kðDÞE½XuD Š. Stochastic Environmental Research and Risk Assessment (vi) h   VðAu ðX; DÞÞ 1 E½V ðAu ðX; DÞފ ¼ E ¼ kðDÞ kðDÞ R þ1 0 0 D E½VðAu ðX; DÞފ ¼ u SD X ðu Þdu ¼ E½Xu Š. D This proposition shows the relevance of the compound CDF for the first-order indicators considered, playing a key role in the methodology proposed in the following section. 3 Methodology In this section, formal and computational aspects of the proposed methodology are described. The approach is based on a conditional risk assessment, performed at global and local (subregional) scales, with the latter being used, in particular, for risk mapping. Specifically, risk measures are evaluated based on the empirical distributions of selected indicators, derived from simulation of the underlying random field conditional to the available spatial (or spatiotemporal) observations. As mentioned, two special cases of quantile-based risk measures, VaR and AVaR, are considered here for illustration. Formally, let X be the random field representing the phenomenon under study, and let u be the reference threshold level considered for risk assessment. As mentioned before, commonly u may be prescribed as a certain critical level of a direct meaning regarding the physical magnitude and contextual aspects. However, in many cases risk assessment is performed under specification of different scenarios for threshold levels corresponding to percentage (or probabilistic) global exceedance considerations. A natural approach, then, consists in considering the compound CDF FXD for reference. More precisely, as given in Proposition 1 (iii), for a given threshold level u, the value FXD ðuÞ ¼ 1 E½kD ðAu ðX; DÞފ represents the expected area proportion within D with no exceedance over u (i.e, XðsÞ\u). In particular, for any a 2 ð0; 1Þ, the threshold defined by the corresponding ð1 aÞ quantile of the compound CDF FXD , i.e. u1 a :¼ FX D1 ð1 aÞ , satisfies E½kD ðAu1 a ðX; DÞފ ¼ a: (Conversely, for a given threshold u, the corresponding exceedance area proportion is given by au ¼ 1 FXD ðuÞ, with an analogous interpretation). In this approach, this form of threshold level specification is applied to the conditional random field derived from X with respect to the available observed data X ¼ fXðs1 Þ; . . .; Xðsn Þg, denoted as XjX. The corresponding excursion set is, in this case, Au1 a ðXjX; DÞ. In practice, the threshold u1 a for a given a value is estimated from the empirical compound CDF obtained by aggregation of the conditional simulation samples (see Sect. 3.2, step 3). XjX;D For simplicity, let us denote by I 1 a any first-order indicator (such as the exceedance area or the excess volume) which is defined as a function of the realizations on ; ðXjXÞu1 a Þ (equivalently, D of the joint random field ðIuXjX 1 a as a function of the random excursion set family   Au1 a0 ðXjX; DÞ : a0  a ). Now, for a conveniently large number M of simulated realizations of XjX, denoted here as ðXjXÞ½mŠ (m ¼ 1; . . .; M), the corresponding excursion sets n o Au1 a0 ððXjXÞ½mŠ ; DÞ and indicator values m¼1;...;M n o ½mŠ ðXjXÞ ;D I1 a can be identified. Evaluation of difm¼1;...;M ferent risk measures is then performed based on the empirical CDF, F^I , obtained from the latter. For instance, assume that a threshold level specification is given as u1 a1 , for a certain prescribed exceedance area proportion a1 . Then, the corresponding ‘empirical’ VaR and AVaR values XjX;D a1 of the I 1 indicator can be derived, for any confidence XjX;D XjX;D level 1 a2 , as VaR1 a2 ðI 1 a1 Þ and AVaR1 a2 ðI 1 a1 Þ. In a spatiotemporal scenario, this approach can be applied to the spatial cross-sections, thus allowing, by conditional simulation, the derivation of dynamic predictive risk assessment according to different temporal horizons. 3.1 Local analysis and risk mapping For a fixed threshold level u, possibly determined as u1 a1 for a certain specified global exceedance area proportion a1 , a local analysis can be performed regarding the spatial distribution of the exceedances, by applying the above described procedure restricted to subregions, D0 , of the domain. For simplicity, among different approaches for implementation, here we assume that D is a rectangular domain and D0 is defined by the restrictions of D based on a rectangular sliding window of a certain size, which is moved in the directions of the coordinate axes at regular steps with a certain degree of overlapping. The evaluation of the selected indicators and risk measures is then performed on each particular D0 subdomain, allowing the construction of ‘local risk’ maps on D at different resolution scales. The spatiotemporal extension can be similarly carried out. 3.2 Further computational aspects In practice, simulation of the underlying random field, conditional on a set of observations at given locations, is 123 Stochastic Environmental Research and Risk Assessment performed based on a regular grid. The nodes are identified as the centroids of cells covering the spatial domain. Computation of the empirical compound CDF based on the simulated realizations is obtained by discretization of the integral (1), as follows: (1) Based on a grid of N1  N2 cells for a (rectangular) domain D, generate M independent replicates of the considered random field XjX (hereafter denoted as X): f X½ði; jÞ; mŠ : i ¼ 1; . . .; N1 ; j ¼ 1; . . .; N2 g; m ¼ 1; . . .; M; where X[(i, j); m] represents the value obtained for the conditional random field (XjX) in the m th simulated realization at the location corresponding to the grid node (i, j). Let mði;jÞ ðxÞ be the number of values X½ði; jÞ; mŠ  x at the node (i, j), and P 1 P N2 mðxÞ :¼ Ni¼1 j¼1 mði;jÞ ðxÞ, i.e. the accumulated number of values X½ði; jÞ; mŠ  x in the aggregated sample f X½ði; jÞ; mŠ : m ¼ 1; . . .; M; i ¼ 1; . . .; N1 ; j ¼ 1; . . .; N2 g; which can be seen as N1 [ N2 [ fX½ði; jÞ; mŠ : m ¼ 1; . . .; M g; i¼1 j¼1 i.e. the union of the set of replicates at each point (i, j), or as M [ fX½ði; jÞ; mŠ : i ¼ 1; . . .; N1 ; j ¼ 1; . . .; N2 g; m¼1 (2) i.e. the union of the m sample realization replicates. Empirical CDF at each node (i, j): assign, for each x 2 R, F^ði;jÞ ðxÞ ¼ (3) 8 0 > > <m > M > : 1 if x\X½ði; jÞ; ð1ފ if X½ði; jÞ; ðmފ  x\X½ði; jÞ; ðm þ 1ފ; m ¼ 1; . . .; M 1 if x  X½ði; jÞ; ðmފ where fX½ði; jÞ; ðmފ : m ¼ 1; . . .; M g is the ordered sample at (i, j). Empirical compound CDF: assign, for each x 2 R, F^XD ðxÞ ¼ N1 X N2 X i¼1 j¼1 1 F^ði;jÞ ðxÞ N1 N2 N2 N1 X X mði;jÞ ðxÞ 1 ¼ M N1 N2 i¼1 j¼1 ¼ 123 mðxÞ : MN1 N2 Excursion sets, and corresponding evaluation of indicators, are approximated on a cell aggregation basis. Here, the package RandomFields of the R Statistical Computing software is used for conditional simulation. The MATLAB environment is used for the remaining computational and graphical aspects. 4 Illustration based on simulations In this section, significant aspects of the quantile-based risk measures VaR and AVaR applied to some indicators of interest, as mentioned in Sects. 2 and 3, are analyzed under different scenarios in the purely spatial and spatiotemporal frameworks by means of a simulation study. For simplicity, hereafter X implicitly stands, when appropriate, for the conditional RF XjX on given observations. A percentile level 1 a1 for the empirical compound CDF FXD , defining an excursion set over threshold u1 a1 , and a confidence level 1 a2 for the empirical CDF of each specific indicator, are considered. Global and local measurements of these risk measures are illustrated and compared. Moreover, in all the cases considered, and for illustration purposes, it must be understood that the random field of interest is the discrete version of the underlying continuous-domain random field consisting of the variables corresponding to a regular grid of a pre-fixed resolution. In this sense, we must emphasize that aspects related to regularity conditions are not involved, except for the fact that the results obtained are intrinsically related to the size of the grid cells (see Adler 1981, Ch. 8, for a discussion in reference to this issue). 4.1 Spatial analysis Three classes of well-known random field models are considered: generalized Cauchy, Student and Fisher–Snedecor. From the global and local measurements, significant interpretations and comparisons, depending on the structural characteristics of the underlying model, are derived from the analysis performed. The Cauchy class model (see Gneiting and Schlather 2004) is defined by the homogeneous and isotropic covariance function CðhÞ ¼ r2 ð1 þ jhja Þ b=a ; a 2 ð0; 2Š and b [ 0: This class allows a separately characterization of local variability and dependence ranges. For a random field on Rd the parameter a determines the fractal dimension of realizations, D ¼ d þ 1 a=2. If b 2 ð0; 1Þ, the process has long memory, with Hurst coefficient H ¼ 1 b=2. Stochastic Environmental Research and Risk Assessment The Student and Fisher–Snedecor class models (see Leonenko and Olenko 2014) are defined as g1 ðxÞ Tn ðxÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ðStudentÞ 1 2 ðxÞ þ


þ g2 ðxÞ g 2 nþ1 n  2  1 g1 ðxÞ þ


þ g2m ðxÞ m  ðFisher SnedecorÞ Fm;n ðxÞ ¼ 1  2 2 n gmþ1 ðxÞ þ


þ gmþn ðxÞ where g1 ; . . .; gmþn are independent copies of a meansquare continuous homogeneous isotropic zero-mean Gaussian random field. For this illustration, the Gaussian Cauchy class model under four scenarios, ða; bÞ ¼ ð0:5; 0:1Þ, (2, 0.1), (0.5, 0.9) and (2, 0.9), with variance r2 ¼ 0:1, and the Student and Fisher–Snedecor classes based on the Cauchy class model under these four scenarios, are considered. For each case, 200 realizations are simulated on the square ½ 100; 100Š2 , both unconditionally and conditionally. More specifically, for the conditional approach, the parameter values are fixed and an unconditional simulated realization of the random field X is generated at 50 pre-fixed locations randomly chosen within the domain. Then, adopting these values as the observed data, realizations are generated by conditional simulation using the same parameter values. Results for various indicators including areas, volumes, and volumes relative to areas (total and per connected component), are derived and compared between the different scenarios, at global and local scales. Figure 1 displays, as an example, simulated realizations based on the four scenarios considered for the Cauchy class, allowing the visualization of the dependence structure and fractality features according to the parameter values. For instance, both aspects jointly depict a particularly contrastive behaviour for the cases ða; bÞ ¼ ð0:5; 0:9Þ and (2, 0.1). The effect of conditioning can be seen comparing to the corresponding plot, for each case, in Fig. 2. Furthermore, as illustrated in Figs. 3 and 4, long-range dependence determined by smaller b provides lower variability in the neighborhood of the conditioning points. The Fig. 1 Simulated realizations on the square ½ 100; 100Š2 from Cauchy model with r2 ¼ 0:1. a (a,b) = (0.5,0.1). b (a,b) = (2,0.1). c (a,b) = (0.5,0.9). d (a,b) = (2,0.9) 123 Stochastic Environmental Research and Risk Assessment Fig. 2 Simulated realizations on the square ½ 100; 100Š2 from Cauchy model with r2 ¼ 0:1, conditional to given values at 50 pre-fixed points for each model. a (a,b) = (0.5,0.1). b (a,b) = (2,0.1). c (a,b) = (0.5,0.9). d (a,b) = (2,0.9) heavy-tailed behaviour of Student and Fisher–Snedecor models is clearly observed in Figs. 5 and 6, respectively, where the different scales of variation can be also noticed; to help visualization, nonlinearly rescaled representations of the generated values are also displayed, using the transformations gðxÞ ¼ sgnðxÞ lnð1 þ jxjÞ for Student case, and hðxÞ ¼ lnðxÞ for Fisher–Snedecor case. 4.1.1 Global measures Based on the empirical compound CDF FXD , obtained from the 200 simulated realizations, for the three models and under the different parameter values, the empirical distributions of various characteristics of excursion sets are analyzed in terms of related VaR and AVaR measures. Specifically, the aim is to assess the sensitivity of these risk measures for varying thresholds, u1 a1 , and varying confidence levels, 1 a2 , under different scenarios. 123 In Fig. 7 three indicators are analyzed for the Cauchy model with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1, namely the absolute excess volume V, the excess volume relative to the exceedance area V Au , and the average relative excess volume per connected component Vcc . Plot (a) shows the rate of decay of V for increasing threshold u1 a1 , for different confidence levels 1 a2 . As expected, for fixed u1 a1 , AVaR1 a2 ðVðAu1 a1 ðX; DÞÞÞ increases as a2 tends to 0. For this particular model, a similar behaviour occurs for V Au , as displayed in plot (b), which indicates that the excess volume decreases faster than the exceedance area for increasing threshold. However, as seen in plot (c), due to the effect of fragmentation, the indicator Vcc changes to increase for higher threshold u1 a1 values. Regarding indicator V Au , it is important to emphasize that, in the case of heavy-tailed scenarios, the excess volume may decrease slower than the exceedance area for increasing threshold, as shown in Fig. 8. Stochastic Environmental Research and Risk Assessment Fig. 3 Surface plots for the map of variances from Cauchy model. a (a,b) = (0.5,0.9). b (a,b) = (2,0.1) Fig. 4 Level plots for the map of variances (top plots) and corresponding level plot of differences of variances between scenarios (bottom plot) from Cauchy model. a (a,b) = (0.5,0.9). b (a,b) = (2,0.1). c Differences between a and b 123 Stochastic Environmental Research and Risk Assessment Fig. 5 Simulated realizations on the square ½ 100; 100Š2 from Student model, (top plots) at original scale and (bottom plots) using the scale transformation gðxÞ ¼ sgnðxÞ lnð1 þ jxjÞ, for ða; bÞ ¼ ð0:5; 0:9Þ and (2, 0.1), with r2 ¼ 0:1. a (a,b) = (0.5,0.9). b (a,b) = (2,0.1). c (a,b) = (0.5,0.9). d (a,b) = (2,0.1) 4.1.2 Local measures A similar study is performed locally, based on sliding windows, for assessment of regional variations of risk. The scale of analysis and degree of interpolation are determined by the size of the window and the sliding step. Here we use a square window of size 33  33, and step equal to 11 cells. Risk surfaces and maps are drawn for different characteristics, models and parameter values. For illustration, in all the cases, an excursion set is determined for a fixed percentile 1 a1 ¼ 0:90, based on the compound CDF FXD . Then, the VaR and AVaR risk measures for a fixed confidence level 1 a2 ¼ 0:95 are calculated for the local 0 0 indicators kD and V Au (with A0u ¼ Au ðX; DÞ \ D0 ) at each subdomain D0 . Figures 9 and 10 show the combined effect of conditioning and model structure for two contrastive Cauchy scenarios. In Fig. 11, as expected, the AVaR surface is above the corresponding VaR surface, displayed in particular for 0 indicator V Au ; however, the conditioning effect makes the relative differences between both measures to vary locally. 123 As noticed in Fig. 12, the heavy-tailed behaviour of the Student and Fisher–Snedecor random field models is reflected in a higher level of local heterogeneity in the relative excess volume compared to the corresponding exceedance area risk maps. 4.2 Spatiotemporal analysis The methodology is also applied to conditional evaluation of risk in a spatiotemporal context. For illustration, the blur-generated non-separable space–time model introduced by Brown et al. (2000) (see also Angulo and Madrid 2014, for an extended version with dynamic spatial deformation) is used: Yðs; tÞ ¼ aY½hŠðs; t 1Þ þ Zðs; tÞ; ð2Þ where Y½hŠðs; t 1Þ ¼ h  Yðs; t 1Þ ¼ Z hðs; s0 ÞYðs0 ; t 1Þds0 ; Rd ð3Þ Stochastic Environmental Research and Risk Assessment Fig. 6 Simulated realizations on the square ½ 100; 100Š2 from Fisher–Snedecor model, (top plots) at original scale and (bottom plots) using the scale transformation hðxÞ ¼ lnðxÞ, for ða; bÞ ¼ ð0:5; 0:9Þ and (2, 0.1), with r2 ¼ 0:1. a (a,b) = (0.5,0.9). b (a,b) = (2,0.1). c (a,b) = (0.5,0.9). d (a,b) = (2,0.1) with Yðs; 0Þ ¼ Y0 being a spatial random field on Rd , Z a spatiotemporal random field on Rd  R assumed to be white in time, h a blurring kernel, and a a constant such that 0\a\1. Here, in the case d ¼ 2, a Cauchy class model for Z with ða; bÞ ¼ ð2; 0:1Þ and r2 ¼ 0:1, a homogeneous Gaussian blurring kernel with variance r2h ¼ 10, and a ¼ 0:9 are considered. Given a fixed complete realization at t ¼ 0 obtained in stationary regime, 200 replicates are conditionally simulated for evolution at t ¼ 1; 2; 3. In this case, a square window of size 21  21, and step equal to 6 cell units is used for local analysis. Predictive risk assessment is then peformed for different characteristics, at horizons 1, 2 and 3. Specifically, in Figs. 13, 14 and 15, corresponding risk maps obtained for 0 0 0 the indicators kD , V D and V Au based on AVaR are displayed. As expected, the levels of risk increase for larger horizon, with this effect varying locally according to the spatiotemporal dependence structure of the underlying process. Also, these maps clearly reflect the dynamic behaviour in relation to the conditioning observation at t ¼ 0. 5 Conclusion In this paper, a flexible methodological approach for spatial and spatiotemporal risk assessment of threshold exceedances is introduced. Using the empirical distribution, derived from conditional simulation on given observations, of selected indicators related to geometrical characteristics of the excursion sets, quantile-based risk measures such as Value-at-Risk and Average Value-at-Risk are evaluated at global or subregional scales. The latter allow the construction of (static or dynamic) risk maps at different resolution levels. In particular, first-order indicators such as absolute or relative versions of exceedance area and excess volume are formally shown to have a meaningful formalization in terms of the compound cumulative (and its reciprocal decumulative) distribution function, which further provides a global quantile correspondence and interpretation regarding the determination of the threshold exceedance level of interest for study. Computational aspects of the implementation of the procedure are described. Different issues related to the characteristics of the underlying random field, in the spatial case, are 123 Stochastic Environmental Research and Risk Assessment Fig. 7 Level plots of AVaR1 a2 ðI ðAu1 a1 ðX; DÞÞÞ for a I ¼ V, b I ¼ V Au and c I ¼ Vcc based on varying (u1 a1 ; 1 a2 ) values, from Cauchy model with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1 123 Fig. 8 Level plots of AVaR1 a2 ðV Au ðA1 a1 ðX; DÞÞÞ based on varying (u1 a1 ; 1 a2 ) values from a Cauchy, b Student and c Fisher– Snedecor models with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1 Stochastic Environmental Research and Risk Assessment 0 0 Fig. 9 Local measurements of VaR for indicators a VaR0:95 ðkD Þ and b VaR0:95 ðV Au Þ from Cauchy class model with ða; bÞ ¼ ð0:5; 0:9Þ, r2 ¼ 0:1 0 0 Fig. 10 Local measurements of VaR for indicators a VaR0:95 ðkD Þ and b VaR0:95 ðV Au Þ from Cauchy class model with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1 0 Fig. 11 Local measurements of a VaR and b AVaR for indicator V Au from Cauchy class model with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1 123 Stochastic Environmental Research and Risk Assessment 0 0 Fig. 12 Local measurements of AVaR0:95 ðkD Þ (left) and AVaR0:95 ðV Au Þ (right) from a Cauchy, b Student and c Fisher–Snedecor models with ða; bÞ ¼ ð2; 0:1Þ, r2 ¼ 0:1 123 Stochastic Environmental Research and Risk Assessment 0 Fig. 13 Local measurements of AVaR0:95 ðkD Þ from model (2)–(3) at horizons (left to right) t ¼ 1; 2; 3. Top plot represents the observed realization at t ¼ 0 0 Fig. 14 Local measurements of AVaR0:95 ðV D Þ from model (2)–(3) at horizons (left to right) t ¼ 1; 2; 3. Top plot represents the observed realization at t ¼ 0 123 Stochastic Environmental Research and Risk Assessment 0 Fig. 15 Local measurements of AVaR0:95 ðV Au Þ from model (2)–(3) at horizons (left to right) t ¼ 1; 2; 3. Top plot represents the observed realization at t ¼ 0 analyzed and discussed through simulated examples. These refer to the nature of the tails of the marginal distributions, as well as to the range of spatial dependence and local variability of realizations. In the spatiotemporal case, application to predictive risk for different horizons is illustrated based on a wellknown autoregressive model involving blur-generated diffusivity. The methodology proposed can be applied in a wide variety of situations, with minimum requirements except for the ability to perform conditional simulation in the case under study. A number of relevant directions for continuing research are open, among which we particularly stress the following: (1) sensitivity analysis with respect to different sources of perturbation and model specifications; (2) generalized thresholds (functional, stochastic, test-function based, etc.) (3) clustering and dynamics of ‘hotspots’; (4) analytical and methodological developments for secondand higher-order indicators, and (5) extension and application of state distortion and space deformation approaches in relation to model complexities. Acknowledgements The authors are particularly grateful to the Associate Editor and the Reviewers for their detailed comments and constructive suggestions, which have led to significant improvements in the final version of the manuscript. This work has been partially 123 supported by Grants MTM2012-32666 and MTM2015-70840-P of Spanish MINECO/FEDER, EU. References Adler RJ (1981) The geometry of random fields. Wiley, Chichester Adler RJ (2008) Some new random field tools for spatial analysis. Stoch Environ Res Risk Assess 22:809–822 Adler RJ, Taylor JE (2007) Random fields and geometry. Springer, New York Adler RJ, Samorodnitsky G, Taylor JE (2010) Excursion sets of three classes of stable random fields. Adv Appl Probab 42:293–318 Adler RJ, Samorodnitsky G, Taylor JE (2013) High-level excursion set geometry for non-Gaussian infinitely divisible random fields. Ann Probab 41:134–169 Angulo JM, Madrid AE (2010) Structural analysis of spatio-temporal threshold exceedances. Environmetrics 21:415–438 Angulo JM, Madrid AE (2014) A deformation/blurring-based spatiotemporal model. Stoch Environ Res Risk Assess 28:1061–1073 Artzner P, Delbaen F, Eber JM, Heath D (1999) Coherent measures of risk. Math Fin 9:203–228 Azäis JM, Wschebor M (2009) Level sets and extrema of random processes and fields. Wiley, Chichester Bernardi M, Durante F, Jaworski P, Petrella L, Salvadori G (2018) Conditional risk based on multivariate hazard scenarios. Stoch Environ Res Risk Assess 32(1):203–211 Brown PE, Karesen KF, Roberts GO, Tonellato S (2000) Blurgenerated non-separable space–time models. J R Stat Soc Ser B 62:847–860 Stochastic Environmental Research and Risk Assessment Chiu SN, Stoyan D, Kendall WS, Mecke J (2013) Stochastic geometry and its applications. Wiley, Chichester Christakos G (1992) Random field models in earth sciences. Academic Press, San Diego Christakos G (2000) Modern spatio-temporal geostatistics. Oxford University Press, New York Christakos G, Hristopulos DT (1996) Stochastic indicators for waste site characterization. Water Resour Res 32:2563–2578 Christakos G, Hristopulos DT (1997) Stochastic indicator analysis of contaminated sites. J Appl Probab 34:988–1008 Craigmile PF, Cressie N, Santner TJ, Rao Y (2005) A loss function approach to identifying environmental exceedances. Extremes 8:143–159 Filipović D, Kupper M (2008) Optimal capital and risk transfers for group divesification. Math Fin 18:55–76 Föllmer H, Schied A (2002) Convex measures of risk and trading constraints. Fin Stoch 6:429–447 Föllmer H, Schied A (2016) Stochastic finance. An introduction in discrete time. Walter de Gruyter GmbH & Co. KG, Berlin French JP, Sain SR (2013) Spatio-temporal exceedance locations and confidence regions. Ann Appl Stat 7:1421–1449 Frittelli M, Gianin E (2002) Putting order in risk measures. J Bank Fin 26:1473–1486 Gneiting T, Schlather M (2004) Stochastic models that separate fractal dimension and the Hurst effect. SIAM Rev 46:269–282 Haier A, Molchanov I, Schmutz M (2016) Intragroup transfers, intragroup divesification and their risk assessment. Ann Fin 12:363–392 Kleinow J, Moreira F, Strobl S, Vähämaa S (2017) Measuring systemic risk: a comparison of alternative market-based approaches. Fin Res Lett 21:40–46 Klüppelberg C, Straub D, Welpe IM (eds) (2014) Risk—a multidisciplinary introduction. Springer, Berlin Lahiri S, Kaiser MS, Cressie N, Hsu NJ (1999) Prediction of spatial cumulative distribution functions using subsampling. J Am Stat Assoc 94:86–97 Leonenko N, Olenko A (2014) Sojourn measures of Student and Fisher–Snedecor random fields. Bernoulli 20:1454–1483 Li QQ, Li YP, Huang GH, Wang CX (2018) Risk aversion based interval stochastic programming approach for agricultural water management under uncertainty. Stoch Environ Res Risk Assess 32(3):715–732 Madrid AE, Angulo JM, Mateu J (2012) Spatial threshold exceedance analysis through marked point processes. Environmetrics 23:108–118 Madrid AE, Angulo JM, Mateu J (2016) Point pattern analysis of spatial deformation and blurring effects on exceedances. J Agric Biol Environ Stat 21:512–530 Piterbarg VI (1996) Asymptotic methods in the theory of Gaussian processes and fields. American Mathematical Society, Providence Souza SRSD, Silva TC, Tabak BM, Guerra SM (2016) Evaluating systemic risk using bank default probabilities in financial networks. J Econ Dyn Control 66:54–75 Vanmarcke E (2010) Random fields. Analysis and synthesis. World Scientific, Singapore Wright DL, Stern HS, Cressie N (2003) Loss functions for estimation of extrema with an application to disease mapping. Can J Stat 31:251–266 Yaglom AM (1987a) Correlation theory of stationary and related random functions I—basic results. Springer, New York Yaglom AM (1987b) Correlation theory of stationary and related random functions II—supplementary notes and references. Springer, New York Yang Y, Christakos G (2015) Spatio-temporal characterization of ambient PM2.5 concentrations in Shandong province (China). Environ Sci Technol 49:13431–13438 Zhang J, Craigmile PF, Cressie N (2008) Loss function approaches to predict a spatial quantile and its exceedance region. Technometrics 50:216–227 123