Correlation of powers of Hüsler–Reiss vectors and Brown–Resnick fields, and application to insured wind losses

Hüsler–Reiss vectors and Brown–Resnick fields are popular models in multivariate and spatial extreme-value theory, respectively, and are widely used in applications. We provide analytical formulas for the correlation between powers of the components of the bivariate Hüsler–Reiss vector, extend these to the case of the Brown–Resnick field, and thoroughly study the properties of the resulting dependence measure. The use of correlation is justified by spatial risk theory, while power transforms are insightful when taking correlation as dependence measure, and are moreover very suited damage functions for weather events such as wind extremes or floods. This makes our theoretical results worthwhile for, e.g., actuarial applications. We finally perform a case study involving insured losses from extreme wind speeds in Germany, and obtain valuable conclusions for the insurance industry.


Introduction
Extreme-value theory offers many statistical techniques and models useful in various fields such as finance, insurance and environmental sciences.Max-stable random vectors (e.g., de Haan and Resnick, 1977) naturally arise when extending univariate extremevalue theory to the multidimensional setting, and several parametric multivariate maxstable distributions, such as the Hüsler-Reiss model (Hüsler and Reiss, 1989), have been proposed.Max-stable random fields (e.g., de Haan, 1984;de Haan and Ferreira, 2006;Davison et al., 2012) constitute an infinite-dimensional generalization and are particularly suitable to model the temporal maxima of a given variable at all points in space since they represent the only possible non-degenerate limiting field of pointwise maxima taken over suitably rescaled independent copies of a field (e.g., de Haan, 1984).One famous example is the Brown-Resnick field (Brown and Resnick, 1977;Kabluchko et al., 2009) which, owing to its flexibility, is generally a good model for spatial extremes of environmental variables.Finite-dimensional distributions of the Brown-Resnick field are Hüsler-Reiss distributions so there is a natural and close link between Hüsler-Reiss vectors and Brown-Resnick fields.
Our main theoretical contributions are explicit formulas for the correlation between powers of the components of bivariate Hüsler-Reiss random vectors, analytical expressions of the spatial correlation function of powers of Brown-Resnick fields, and a careful study of its properties; some results are rather technical to obtain.Studying the correlation function of a field is prominent as it naturally appears when computing the variance of the spatial integral of that field (e.g., Koch, 2019b).If the field models an insured cost, its spatial integral models the total insured loss over the integration region, and its variance is thus of interest for any insurance company.The correlation function also explicitly shows up in the standard deviation of the central limit theorem of the field, and is thus key for the behaviour of the spatial integral when the size of the integration region becomes large (e.g., Koch, 2019b).Moreover, despite its drawbacks, correlation is commonly used in the finance/insurance industry, making its study useful from a practical viewpoint.Finally, the criticism that it does not properly capture extremal dependence is somewhat irrelevant here as we consider the correlation between random variables which already model extreme events.
It is often insightful to consider the correlation between various powers of two random variables rather than focusing only on the correlation between these variables.First, applying simple non-linear transformations such as the absolute value or powers before taking the correlation sometimes allows one to detect and characterize a strong dependence that would not have been spotted using otherwise; this partially alleviates the defect that correlation only captures linear dependence.In finance, it is common to look at the autocorrelation of powers of the absolute values of asset returns.Returns generally do not exhibit any significant autocorrelation (e.g., Cont, 2001) whereas their squares or other power values (see, e.g., Ding et al., 1993, who consider powers ranging from 0.125 to 5) show a significantly positive serial correlation.Second, taking powers may be useful for estimation.Let X 1 , X 2 be random variables whose joint distribution function depends on a parameter θ.If an explicit formula is available for the correlation between X β 1 and X β 2 , where β belongs to an appropriate space, then one can think of estimating θ by equating that expression with its empirical counterpart, and searching for the value of β leading to the optimal corresponding estimator.Such an approach may be notably useful for max-stable random fields, for which estimation is arduous.
Considering powers of random variables is also valuable when these variables are used to model the impact of natural disasters such as, e.g., windstorms or floods.According to physics, the total cost arising from damaging wind to a specific structure should increase as the square (e.g., Simiu and Scanlan, 1996) or the cube (e.g., Lamb and Frydendahl, 1991;Emanuel, 2005;Powell and Reinhold, 2007) of the maximum wind speed.Moreover, several studies exploring insured costs have found that power-laws with much higher exponents are appropriate (e.g., Prahl et al., 2012).A commonly used damage function for flood is D(z) = z/(z + 1), where z > 0 is the water level measured in meters (e.g., Hinkel et al., 2014;Prahl et al., 2016) and so the destruction percentage approximately follows a power-law with exponent unity for levels much below one meter.Thus, as max-stable vectors and fields are suited to model componentwise and pointwise maxima, studying their powers is worthwhile for assessing costs from extreme wind or flood events.
In the second part of the paper, we use our theoretical results to study the spatial dependence of insured losses from extreme wind speed for residential buildings over a large part of Germany.We use ERA5 (European Centre for Medium-Range Weather Forecasts Reanalysis 5th Generation) wind speed reanalysis data on 1979-2020 to derive seasonal pointwise maxima, we fit the Brown-Resnick and Smith random fields, and use the appropriate power damage function for the considered region, according to Prahl et al. (2012).The best fitted model leads to a correlation displaying a slow decrease with the distance.We also consider other power values and we find that, for a fixed distance, the correlation between insured costs evolves only slightly with the value of the damage power; this is useful information for insurance companies.
The rest of the paper is organized as follows.Section 2 first briefly reviews Hüsler-Reiss vectors and Brown-Resnick fields, and then details our main theoretical contributions.We present our case study in Section 3, and Section 4 summarizes our main findings and provides some perspectives.All the proofs are gathered in the Appendix.The code and data required to reproduce the results of the case study will be available in a publication on the Zenodo repository.Note that some elements of this article are revised versions of results from Sections 2.2 and 3 and Appendix A of the unpublished work by Koch (2019a).Throughout the paper, designates transposition and N * = N\{0}.

Preliminaries
having standard Fréchet marginals is said to follow the bivariate Hüsler-Reiss distribution (Hüsler and Reiss, 1989) (1) This is a popular and flexible distribution for max-stable random vectors, and the parameter h interpolates between complete dependence (h = 0) and independence (h = ∞).
The i-th component, i = 1, 2, of any bivariate max-stable vector follows the generalized extreme-value (GEV) distribution with location, scale and shape parameters where (Z 1 , Z 2 ) is a max-stable vector with standard Fréchet marginal distributions.
In the following, a max-stable random field with standard Fréchet margins will be called simple.The class of Hüsler-Reiss distributions is tightly linked to the Brown-Resnick random field (Brown and Resnick, 1977;Kabluchko et al., 2009) which is a flexible and widely used max-stable model.It is very suited to model, e.g., extremes of environmental data (e.g., Davison et al., 2012, Section 7.4, in the case of rainfall) as it allows realistic realizations as well as independence when distance goes to infinity.If {W (x)} x∈R d is a centred Gaussian random field with stationary increments and with semivariogram γ W , then the Brown-Resnick random field associated with the semivariogram γ W is defined by where the (U i ) i≥1 are the points of a Poisson point process on (0, ∞) with intensity function u −2 du and the Y i , i ≥ 1, are independent replications of where Var denotes the variance.It is a stationary1 and simple max-stable field whose distribution only depends on the semivariogram (Kabluchko et al., 2009, Theorem 2 and Proposition 11, respectively).Its finite-dimensional distribution functions are Hüsler-Reiss distributions (Kabluchko et al., 2009, Remark 24) and, in particular, for any A commonly used semivariogram is where κ > 0 and ψ ∈ (0, 2] are the range and the smoothness parameters, respectively, and • denotes the Euclidean norm.The Smith random field with positive definite covariance matrix Σ (Smith, 1990) corresponds to the Brown-Resnick field associated with the semivariogram see, e.g., Huser and Davison (2013).
If {X(x)} x∈R d is max-stable, there exist functions η(•) ∈ R, τ (•) > 0 and ξ(•) ∈ R defined on R d , called the location, scale and shape functions, such that where {Z(x)} x∈R d is simple max-stable.In the following, if {X(x)} x∈R d is defined by (7) with {Z(x)} x∈R d being the Brown-Resnick field associated with the semivariogram γ W , then X will be referred to as the Brown-Resnick field associated with the semivariogram γ W and with GEV functions η(x), τ (x) and ξ(x).If, for all x ∈ R d , η(x) = η, τ (x) = τ and ξ(x) = ξ, then X will be termed the Brown-Resnick field associated with the semivariogram γ W and with GEV parameters η, τ and ξ.

Theoretical contributions
Several dependence measures for max-stable vectors and fields have been introduced in the literature: the extremal coefficient (e.g., Schlather and Tawn, 2003), the F-madogram (Cooley et al., 2006) and the λ-madogram (Naveau et al., 2009), among others.Here we propose a new spatial dependence measure which is the correlation of powers of maxstable vectors/fields and not of max-stable vectors/fields themselves.As explained in Section 1, taking power transforms when using correlation is standard practice when dealing with financial time series.For X being defined by ( 2) with (Z 1 , Z 2 ) following the Hüsler-Reiss distribution (1), we study Corr(X β 1 1 , X β2 2 ), where β i ∈ N * such that β i ξ i < 1/2, and Corr denotes the correlation.This allows obtaining the expression of Corr(X β(x 1 ) (x 1 ), X β(x 2 ) (x 2 )), x 1 , x 2 ∈ R 2 , where X is the Brown-Resnick field associated with any semivariogram and with GEV functions η(x), τ (x), ξ(x), and β(x) is a function taking values in N * such that β(x)ξ(x) < 1/2 for all x ∈ R 2 .If those GEV functions and β(x) are not spatially constant, the field {X β(x) (x)} x∈R 2 is not second-order stationary and its correlation function does not only depend on the lag vector.Taking constant GEV and power functions as in the case study is however reasonable when the region considered is fairly homogeneous (in terms, e.g., of altitude, weather influences and distance to a coastline) or not too large.Moreover, every non-stationary random field can be approximated by piecewise stationary fields; see Koch (2019b) and references therein.Therefore, our main focus will be on where X is the Brown-Resnick field with GEV parameters η, τ , ξ, and β ∈ N * such that βξ < 1/2; in this setting, X β is second-order stationary.
Rescaled powers of max-stable random fields constitute appropriate models for the field of insured costs from high wind speeds (see Section 3.1 for details) and so (8) can be viewed as the correlation function of insured wind costs, thus being useful for actuarial practice.The formulas we derive in this section also make possible the estimation of the parameters of Hüsler-Reiss distributions and Brown-Resnick random fields by equalizing the theoretical correlation and the empirical one computed on the dataset (method of moments); this may be investigated in a subsequent work.Appendix B, which deals with simple Brown-Resnick fields, can be useful in this respect.
Before presenting the main results, we recall the importance of correlation for risk assessment in a spatial context, which justifies studying the correlation despite the existence of dependence measures specifically designed for max-stable fields.Anyway, powers of max-stable fields are not necessarily max-stable themselves, making these measures not directly usable.
Denote by C the set of all real-valued and measurable 2 random fields on R 2 having almost surely (a.s.) locally integrable sample paths.Furthermore, let A denote the set of all compact subsets of R 2 with a strictly positive Lebesgue measure and A c be the set of all convex elements of A. For any A ∈ A c , let b A denote its barycenter and λA be the area obtained by applying to A a homothety with center b A and ratio λ > 0.
Let C ∈ C model the insured cost per surface unit triggered by events belonging to a specific class (e.g., European windstorms) during a given period of time.The total insured loss on a given region A ∈ A can thus be modelled by and Theorem 4 in Koch (2019b) yields Hence the correlation is explicitly involved in the variance of the total insured loss, which is a key quantity for an insurance company.Moreover, assuming that C belongs to C, has a constant expectation and satisfies the CLT (see Koch et al., 2019, Section 3.1) (which holds for C = X β if X is the Brown-Resnick field associated with the semivariogram (5) and with GEV parameters η, τ and ξ such that βξ < 1/2), is the standard deviation of the normal distribution appearing in the CLT of C and is thus (Koch, 2019b, Theorems 2 and 5) essential for the asymptotic distribution of L(λA, C) and the asymptotic properties of spatial risk measures induced by the field C and associated with value-at-risk and expected shortfall.The analysis of ( 8) is thereby insightful for the risk assessment of wind damage; the formulas derived in this paper are used in an ongoing study.
As (2) specifies a transformation of simple max-stable random vectors, we first deal with such vectors.In the next theorem, we take a random vector Z = (Z 1 , Z 2 ) following the Hüsler-Reiss distribution (1).If β ∈ R and Z is a standard Fréchet random variable, it is easily shown that Z β has a finite second moment if and only if β < 1/2, which imposes, in order for the covariance Cov(Z β 1 1 , Z β 2 2 ) to exist, that β 1 , β 2 < 1/2.This covariance and other expressions throughout this section involve, for β 1 , β 2 < 1/2, where Γ denotes the gamma function, and, for θ, h > 0, with Φ and φ denoting the standard Gaussian distribution and density functions, respectively.
In order to obtain the next result, we take advantage of the radius/angle decomposition of multivariate extreme-value distributions.
Theorem 1.Let Z = (Z 1 , Z 2 ) follow the Hüsler-Reiss distribution (1) with parameter h.Then, for all Remark 1. Theorem 1, which is a cornerstone of this section, stems from unpublished work in Section 4.5.1 of the PhD thesis by Koch (2014).
We adapt Theorem 1 to the more realistic setting where the margins are general GEV distributions with non-zero shape parameters.The support of such margins possibly includes strictly negative values, and we thus consider powers which are strictly positive integers.
Let Z be a simple max-stable vector with continuous exponent function and let Using similar arguments, we get lim ξ→0 Var(X This result obviously applies if Z follows the Hüsler-Reiss distribution (1).
Next proposition, which is an immediate corollary of Theorem 2, provides all the necessary ingredients for the computation of our dependence measure D X,β in (8).
Proposition 2. Under the same assumptions as in Theorem 2 but with and, for i = 1, 2, The following theorem, which is a direct consequence of (4) and Proposition 2, gives the expression of D X,β .
Note that the case ξ = 0 is easily recovered as explained above.
Remark 2. The combination of (4) and Theorem 2 yields the following more general result than Theorem 3. Let {X(x)} x∈R 2 be the Brown-Resnick field associated with the semivariogram γ W and with GEV functions η(x) ∈ R, τ (x) > 0, ξ(x) = 0, and let β(x) be a function taking values in N * such that β(x)ξ(x) < 1/2 for any x ∈ R 2 .Then, The analytical formulas in Theorems 1, 2, 3, Proposition 2, and Remark 2 allow a more accurate and much faster computation of the respective quantities than using Monte Carlo methods as the involved integrals can be computed fast and with high precision using, e.g., adaptive quadrature.The Smith field being a member of the class of Brown-Resnick fields, Theorem 3 and Remark 2 also apply for X being the Smith field with any covariance matrix.
The influence of the marginal parameters and of the power β merits some theoretical comments.Let Z = (Z 1 , Z 2 ) and X = (X 1 , X 2 ) be as in Theorem 2 and suppose that X 1 and X 2 are a.s.strictly positive (i.e., ξ 1 , ξ 2 > 0 and η 1 − τ 1 /ξ 1 , η 2 − τ 2 /ξ 2 > 0).For η ∈ R, τ > 0 and ξ = 0, the transformation z → η − τ /ξ + τ z ξ /ξ, z > 0, is strictly increasing and the same applies for x → x β , x > 0, with β ∈ N, and z → z β * , z > 0, with 0 < β * < 1/2.Thus, owing to the invariance of the copula of a distribution under strictly increasing transformations of the margins, the copula of (X β 1 1 , X β 2 2 ) is the same whatever the values of β i ∈ N * , and is the same as the copula of (Z However, the correlation between two random variables does not only depend on their copula but also on their margins, and is typically not invariant under non-linear transformations.We do not have equality between Corr(X β 1 1 , X β 2 2 ) and Corr(Z 2 ) in general, as can also be seen directly from the formulas, and this also holds in the particular case where Z, X and β 1 , β 2 are as in Proposition 2 and 8) is not invariant with respect to the marginal parameters η, τ , ξ and the power β.Taking the appropriate values of those quantities is necessary when using D X,β for concrete risk assessment problems, and studying its sensitivity with respect to β is also of interest.The conclusions of this paragraph regarding the correlations are a fortiori true if X 1 , X 2 are not a.s.strictly positive; in that case, even the mentioned equalities of copulas do not hold in general.
We now investigate the behaviour of the function g β,η,τ,ξ defined in (15) in order to derive useful conclusions about D X,β and because we need it in an ongoing work about spatial risk measures.The proof of next proposition is appealing as it first involves showing a result (Proposition 6 in Appendix A.4.1) about the correlation order, which is a classical concept of dependence comparison in actuarial risk theory (e.g., Denuit et al., 2005, Section 6.2).
Proposition 5.For all η ∈ R, τ > 0, ξ = 0 and β ∈ N * such that βξ < 1/2, the function g β,η,τ,ξ defined in (15) satisfies By Theorem 3, D X,β (x 1 , x 2 ) depends on x 1 and x 2 through γ W (x 2 − x 1 ) only.As a variogram is a non-negative conditionally negative definite function, it follows from Berg et al. (1984, Chapter 4, Section 3, Proposition 3.3) defines a metric.For many common models of isotropic semivariogram γ W , γ W (x 2 − x 1 ) is a strictly increasing function of x 2 − x 1 , which implies by ( 17) and Proposition 3 that D X,β (x 1 , x 2 ) is a strictly decreasing function of x 2 − x 1 ; such a decrease of the correlation with the distance seems natural.Moreover ( 17) and ( 19) give that lim x 2 −x 1 →0 D X,β (x 1 , x 2 ) = 1, and ( 17) and (20) imply, provided The faster the increase of γ W to infinity, the faster the convergence of D X,β (x 1 , x 2 ) to 0. These results are consistent with our expectations.For a function

Case study
We focus on insured losses from wind extremes for residential buildings over a large part of Germany, more precisely over the rectangle from 5.75 • to 12 • longitude and 49 • to 52 • latitude (see Figure 1).We apply the results developed in Section 2.2 for assessing the spatial dependence of those losses.For the insured cost field, we use the model introduced in Koch (2017, Section 2.3), that is where E is the strictly positive and deterministic field of insured value per surface unit, D : R → [0, 1] is the damage function, and X is the model for the random field of the environmental variable generating risk.Applying the damage function D to X allows getting at each site the insured cost ratio, which, multiplied by the insured value, gives the corresponding insured cost.We assume the risk to be generated by wind speed maxima and we model the latter with a Brown-Resnick and a Smith max-stable model.Section 3.1 outlines and thoroughly justifies the power damage function D that we will use.In Section 3.2 we describe the wind speed data and perform model estimation, selection and validation.Finally, we apply in Section 3.3 the results of Section 2.2 using the derived insured cost model.

Power damage function
We consider the damage function where β ∈ N * and c 1 > 0. The quantity c 1 corresponds to the wind speed above which the insured cost ratio equals unity and can be assumed to be much larger than possible wind speed values in Germany, especially as the distribution of wind speed maxima is bounded in this application (see below).In practice it is hence equivalent to take and this is our choice in the following.
Power functions are perfectly suited to the case of wind.The total cost for a specific structure should increase as the square or the cube of the maximum wind speed since wind loads and dissipation rate of wind kinetic energy are proportional to the second and third powers of wind speed, respectively.For arguments supporting the use of the square, see, e.g., Simiu and Scanlan (1996, Equations (4.7.1), (8.1.1)and (8.1.8)and the interpretation following Equation (4.1.20)).Regarding the cube, see, among others, Lamb and Frydendahl (1991, Chapter 2, p.7) where the cube of the wind speed appears in the severity index, and Emanuel (2005).In his discussion of the paper by Powell and Reinhold (2007), Kantha (2008) states that wind damage for a given structure must be proportional to the rate of work done (and not the force exerted) by the wind and therefore strongly argues in favour of the cube.In addition to this debate about whether the square or cube is more appropriate for total costs, several studies in the last two decades have found power-laws with much higher exponents when insured costs are considered.For instance, Prahl et al. (2012) find powers ranging from 8 to 12 for insured losses on residential buildings in Germany (local damage functions).Prahl et al. (2015) argue that, if the total cost follows a cubic law but the insurance contract is triggered only when that cost exceeds a strictly positive threshold (e.g., in the presence of a deductible), then the resulting cost for the insurance company is of power-law type but with a higher exponent.We have validated this statement using simulations and observed that the resulting exponent depends on the threshold (not shown).
Several authors (e.g., Klawa and Ulbrich, 2003;Pinto et al., 2007;Donat et al., 2011) use, even in the case of insured losses, a cubic relationship that they justify with the physical arguments given above.However, they apply the third power to the difference between the wind speed value and a high percentile of the wind distribution and not to the effective wind speed; as shown by Prahl et al. (2015, Appendix A3), this is equivalent to applying a much higher power to the effective wind speed.Note that exponential damage functions are sometimes also encountered in the literature (e.g., Huang et al., 2001;Prettenthaler et al., 2012); we do not consider such functions here.
According to Prahl et al. (2012) who use ( 22) as well, a spatially-constant exponent of 10 seems appropriate in our region for insured losses on residential buildings; see their As will be seen, the normalization does not play any role in our application.

Wind data
We consider hourly maxima of the 3 s wind gust at 10 m height (as defined by the World Meteorological Organization) from 1 January 1979 08:00 central European time (CET) to 1 January 2020 at 00:00 CET.This is publicly available data from the European Centre for Medium-Range Weather Forecasts (ECMWF); more precisely we use the "10 m wind gust since previous post-processing" variable in the ERA5 (ECMWF Reanalysis 5th Generation) dataset.The covered region is a rectangle from 5.75 • to 12 • longitude and 49 • to 52 • latitude and the resolution is 0.25 • latitude and 0.25 • longitude, leading to 338 grid points.We randomly choose 226 of them to fit the models and use the remaining 112 for model validation; see Figure 1.This area encompasses the Ruhr region in Germany and is associated with high residential insured values per surface unit.We derive at each grid point the 42 seasonal (from October to March) maxima and fit the models to the resulting pointwise maxima.For the first and last season, the maxima are computed over January-March and October-December, respectively.Focusing on October-March allows us to get rid of seasonal non-stationarity in the wind speed time series and to mainly account for winter storms rather than intense summer thunderstorms.

Model
We consider both the Brown-Resnick field with semivariogram (5) and the Smith field.
As mentioned above, max-stable models are very natural ones for pointwise maxima, and the Brown-Resnick field generally shows good performance on environmental data.We model the location, scale and shape parameters as constant across the region, which is reasonable here (this can be explained by the homogeneity in terms of altitude and weather influences).Using trend surfaces for these parameters rather than fitting them separately at each grid point is standard practice as it reduces parameter uncertainty, allows a joint estimation of all marginal and dependence parameters in a reasonable amount of time and enables prediction at sites where no observations are available.Allowing anisotropy in the semivariogram of the Brown-Resnick model would be pertinent but would not modify our main conclusions.Isotropy already leads to a very satisfying model and makes our dependence measure (8) isotropic in the original space, which facilitates our discussions in Section 3.3.Both models are fitted using maximum pairwise likelihood (e.g., Padoan et al., 2010) implemented in the fitmaxstab function of the SpatialExtremes R package (Ribatet, 2020); marginal and dependence parameters were jointly estimated using the Nelder-Mead algorithm with a relative convergence tolerance of 1.49 × 10 −8 .We then perform model selection by minimization of the composite likelihood information criterion (CLIC); see Varin and Vidoni (2005).According to that criterion, the Brown-Resnick field is the most compatible with the data; see Figure 2 shows that the theoretical pairwise extremal coefficient function of the fitted Brown-Resnick model agrees reasonably well with the empirical pairwise extremal coefficients for the validation grid points.It is slightly above their binned estimates when those are computed using the empirical distribution functions.This small underestimation of the spatial dependence likely comes from the choice of parsimonious trend surfaces for the location, scale and shape parameters, and disappears when we compute the empirical extremal coefficients using the marginal parameters' estimates.Overall Figure 2 indicates that the proposed model fits the extremal dependence structure of the data fairly well.Figure 3 is complementary as it assesses both the marginal and dependence components; it shows that the distributions of several summary statistics in various dimensions are very similar for our model and the data.Finally Figure 4 suggests, for two seasons with different ranges of values, that realizations from our model have similar patterns as observed pointwise maxima, although being slightly rougher.The combination of these goodness-of-fit assessments shows that the proposed model is well-suited to our data, and so that this case study is useful in practice.
Having shown that our model performs well, we fit it to the data corresponding to all grid points in order to get as accurate parameters' estimates as possible; see Table 2. Our estimates are in line with general findings on wind speed extremes.Many studies point out that the shape parameter ξ is usually slightly negative, entailing that the distribution of wind speed maxima has a finite right endpoint.E.g., Ceppi et al. (2008) obtain a ξ ranging from −0.2 to 0 by fitting a generalized Pareto distribution (GPD) to in situ observations over Switzerland.Similarly, Della-Marta et al. ( 2007) fit a GPD to ERA-40 (ECMWF Reanalysis originally intended as a 40-year reanalysis) windstorms data over Europe and find negative values, between −0.1 and −0.3 on most of land areas; see their Figure 4.15.Typical values for the location and scale parameters η and τ for yearly maxima over Europe are about 25 m s −1 and 3 m s −1 , respectively; e.g., considering annual maxima at 35 weather stations in the Netherlands, Ribatet (2013) obtains trend surfaces whose intercepts are about 27 m s −1 for η and 3.25 m s −1 for τ .Finally, a value of the smoothness parameter ψ between 0.2 and 1 seems reasonable; e.g, Ribatet (2013) obtains 0.24 on the Netherlands data and, on similar ones, Einmahl et al. (2016) find 0.40.We obtain a higher value perhaps because reanalysis data tend to be smoother than in situ observations.2: Parameters' estimates (standard errors inside brackets) when using all grid points for the fit.

Results
Using (21), ( 22) and the facts that E(x) > 0 for any x ∈ R 2 and c β 1 > 0, we get Therefore, our dependence measure (8) naturally appears in concrete assessments of the spatial risk associated with extreme wind speed.In this section, we thoroughly study the evolution of D X,β (x 1 , x 2 ) with respect to x 2 − x 1 , where X is the Brown-Resnick model fitted to the data in Section 3.2.2,i.e., with semivariogram (5) and parameters in Table 2, and where β has the proper value on our region, i.e., 10.The integral in I β 1 ,β 2 (see ( 10)) has no closed form and therefore a numerical approximation is required.For this purpose, we use adaptive quadrature with a relative accuracy of 10 −13 .Figure 5 shows a decrease of D X,β from 1 to 0 as the Euclidean distance increases, in agreement with our theoretical results of Section 2.2.The decrease is quite slow owing to fairly large range κ and rather low smoothness ψ.For two grid points 5 • and 10 • away, D X,10 is still as high as 0.65 and 0.48, respectively.The latter conclusion is however hypothetical as the largest distance between two grid points in our region is about 6.93 • ; fitting our model on a wider region would be possible, but the assumption of spatially-constant GEV parameters and power might be less suitable.This slow decrease points out the need for an insurer to cover a wider region than the one considered here in order to benefit from sufficient spatial diversification.As already mentioned, various values (basically from 2 to 12) of damage powers have been proposed in the literature and the appropriate one may depend on the insurance contract.Moreover, as explained in Section 1, taking powers (such as the square) of the variables of interest is worthwhile when using correlation as dependence measure.For example, if the true power is 6, it may also be valuable to study Corr([X 6 (x 1 )] 2 , [X 6 (x 2 )] 2 ) = Corr(X 12 (x 1 ), X 12 (x 2 )).For these reasons, investigating how D X,β (x 1 , x 2 ) varies with β for a given max-stable model X and various values of x 1 − x 2 is useful.Figure 6 shows  that whatever the model considered (including the one fitted to our data) and for any given Euclidean distance, D X,β is only faintly sensitive to the value of β; more precisely, it very slightly increases in a concave way with β.On top of being potentially insightful for the understanding of max-stable fields, this finding is valuable for actuarial practice as it shows that making a small error on the evaluation of β is not very impactful as far as correlation is concerned.Nonetheless this does not imply that the computations should be done with β = 1 regardless of the true power value.First, although evolving little with β, our dependence measure is not constant with β and so using the right value is recommended for accuracy.Second, β strongly affects Var(X β (x)) for any x ∈ R 2 , and thus for instance the covariance function and the variance in (9).5) and parameters in Table 2.
Although the smoothness parameter ψ has been estimated on the data, we also consider various values since ψ heavily affects the rate of decrease of D X,β as the distance between the two sites increase, and thereby the rate of spatial diversification for an insurance company.This allows us to figure out the impact of the use of rougher or smoother data, of estimation error, and of model misspecification.We take ψ = 0.5, 0.81, 1.5, 2; the value 0.81 is the one we obtained on our data, ψ = 2 corresponds to the Smith field with Σ = I 2 (see ( 6)), ψ = 1.5 is intermediate between these two settings, and ψ = 0.5 corresponds to a quite rough field.In accordance with the discussion at the end of Section 2.2, Figure 6 shows that D X,β decreases from 1 to 0 as the Euclidean distance increases, and this at a higher rate for larger values of ψ.The decrease is faster for the Smith field than for all Brown-Resnick fields having ψ < 2, and if the true value of ψ is close to 0.5 or even 0.81, using the Smith model leads to a serious underestimation of the dependence between insured costs.The minimum Euclidean distance required for D X,10 to be lower than 0.1 equals 43.60 • for ψ = 0.81, instead of around 9.54 • for ψ = 2 (not shown).
The results outlined in the two previous paragraphs remain qualitatively unchanged with other values of η, τ , ξ, and choosing a specific value for κ does not induce any loss of generality in our study; should κ be different, the appropriate plots would be the same as in Figure 6 with the values on the x-axis multiplied by the ratio between the true value and the one chosen here.
Finally we briefly study the extension of ( 8) where the marginal parameters and the power are site-specific.We consider two sites x 1 , x 2 that are 3 • away, but our findings hold more generally.We successively investigate the effects of a spatially-varying power, location, scale and shape; more precisely we evaluate (18) where X is the Brown-Resnick model with semivariogram (5) • with parameters in Table 2 and β(x 1 ), β(x 2 ) ∈ {1, . . ., 12}.
The ranges for the GEV parameters have been chosen to be approximately centred on the estimates obtained on the data.Figure 7 shows that, for a fixed β(x 1 ), the correlation increases with β(x 2 ) on [1, β(x 1 )] and then decreases.The highest correlation is thus obtained for β(x 2 ) = β(x 1 ) = β, and, as already seen, slightly increases in a concave way when β increases.Also, the higher the difference between β(x 1 ) and β(x 2 ), the lower the correlation.Similar conclusions hold for the scale and shape parameters, although the variations of the correlation are smaller for the chosen range of values.For τ (x 1 ) = τ (x 2 ) = τ , the increase with respect to τ is concave, whereas for ξ(x 1 ) = ξ(x 2 ) = ξ, the increase with respect to ξ is linear.The findings for the location are similar to those for the scale and shape although, for η(x 1 ) = η(x 2 ) = η, the correlation slowly decreases in a concave way as η increases.

Conclusion
Hüsler-Reiss vectors and Brown-Resnick fields are popular and widely used models for componentwise and pointwise maxima.We provide explicit formulas for the correlation between powers of the components of bivariate Hüsler-Reiss vectors and deduce analytical expressions for the correlation function of powers of Brown-Resnick fields.Although extremal models are considered, studying the correlation function makes sense as the latter is required when we are interested in the variance or the asymptotic distribution of the spatial integral of a field, which is typically the case in spatial risk assessment.The application of a power transform to random variables before taking the correlation allows detection of part of non-linear dependence and is therefore common practice in financial time series analysis.Moreover, the relevance of powers as damage functions for natural disasters is largely documented in the literature.In the second part of the paper, we use our theoretical contributions and reanalysis wind gust data to study the spatial dependence of modelled insured losses from extreme wind speeds for residential buildings in Germany.We find that the dependence decreases slowly with the distance and that our dependence measure is not very sensitive to the power value.
The theoretical results obtained here are used in Koch and Robert (2022) as well as in an ongoing study where spatial risk measures (Koch, 2017(Koch, , 2019b) ) are applied to concrete assessment of the risk of impacts from extreme wind speeds.Other potentially interesting applications of the derived expressions are flood risk assessment and momentbased estimation of the parameters of Hüsler-Reiss vectors or Brown-Resnick fields.A more detailed study, both theoretically and numerically, of the correlation function expressed in Remark 2 (non-stationary case) would be welcome, and deriving analytical formulas of (8) for other classes of max-stable fields such as the extremal t model (Opitz, 2013) as well as r-Pareto fields (e.g., de Fondeville and Davison, 2018) would be useful for applications.
Although our insured loss model is supported by the literature, thoroughly assessing its performance on insured loss data is prominent for practice, and this is done in an ongoing work.Finally, in the case study we take the value of β obtained in other papers as given, and an approach that would involve estimating β as well would consist in fitting powers of rescaled max-stable fields to insured loss data directly.

A Proofs
A.1 For Theorem 1 Proof.First, we show the result for h = 0.In that case, Z 1 = Z 2 a.s.(e.g., Hüsler and Reiss, 1989, Section 2).Hence, since Z 1 and Z 2 follow the standard Fréchet distribution, E[Z β i i ] = Γ(1 − β i ), i = 1, 2, and thus Now, we prove the result for h > 0. We have where l denotes the bivariate density of Z.In order to take advantage of the radius/angle decomposition of multivariate extreme-value distributions, we make the change of variable The corresponding Jacobian matrix is written and its determinant is thus det(J Ψ (u, θ)) = u.Therefore, introducing we have Differentiation of (1) yields (see, e.g., Padoan et al., 2010, Equation (4)), for z 1 , z 2 > 0, where Therefore, for any u, θ > 0, l(u, θu) We denote by F s f the Fréchet distribution with shape and scale parameters 1 and s f > 0, i.e., if X ∼ F s f , P(X ≤ x) = exp(−s f /x), x > 0. Using ( 24) and ( 26) and the fact that the density of X ∼ F s f is l f (x) = s f /x 2 exp (−s f /x), we obtain where µ k (F ) stands for the k-th moment of a random variable having F as distribution.
It is immediate to see that which, combined with (27), yields the result.
Using this together with (2) and the binomial theorem yields (13).

A.3 For Proposition 1
Proof.For i = 1, 2, X i,ξ follows the GEV distribution with parameters η i , τ i and ξ, the density of which we denote by f i .Let us assume that ξ ∈ S β 1 ,β 2 ,ε and ξ > 0. We have for all α > 0 and thus for any subset S of (0, ∞).We deal with the second integral in (28), for which there is a potential problem at ∞.We have where we used the change of variable A well-known inequality states that log(z −ξ ) ≥ 1 − 1/z −ξ for any z ≥ 0, which yields z −ξ (log(z −ξ ) − 1) ≥ −1 and thus g (ξ) ≥ 0. Combined with the fact that 0 ≤ z ≤ 1, this gives for any 0 < ξ ≤ ξ * and therefore, taking α = β i (1 + ε), Combining this result with a similar reasoning for the first integral in (28) and using (29) yields sup ξ∈(0 It follows from Billingsley (1999, p.31) that the (X 1,ξ ) ξ , (X 2,ξ ) ξ and (Y ξ ) ξ are uniformly integrable for ξ around 0 (from the right).Now, it is well-known that X i,ξ d → X i,0 , i = 1, 2, which implies by the continuous mapping theorem that X where V is the exponent function of (Z 1 , Z 2 ) .Thus, by continuity of V , and therefore Consequently, the continuous mapping theorem yields and hence, applied again, Y ξ d → Y 0 .Finally, Theorem 3.5 in Billingsley (1999) The result follows immediately.Similar arguments give the same conclusion for ξ < 0.

A.4 For Proposition 3
In this section, we denote by F X the distribution function of any random variable X and by F X 1 ,X 2 the distribution function of any random vector X = (X 1 , X 2 ) .

A.4.1 Preliminary result
We first need the following result.
The field X β is stationary by stationarity of X and has a finite second moment since βξ < 1/2.Accordingly, X β is second-order stationary.Moreover, X β is samplecontinuous and thus, by the same arguments as in the proof of Proposition 1 in Koch et al. (2019), continuous in quadratic mean.Hence, the covariance function of X β is continuous at the origin.It implies by Theorem 3 that which, combined with (16), yields (19).This easily gives lim h→0 g β,η,τ,ξ (h) = g β,η,τ,ξ (0), which implies that g β,η,τ,ξ is continuous at h = 0.The continuity of g β,η,τ,ξ at any h > 0 comes from the fact that the covariance function of a field which is second-order stationary can be discontinuous only at the origin.
Proof.The field X β is obviously measurable.Furthermore, as X has identical univariate marginal distributions, the function x → E[|X(x) β |] is constant and hence locally integrable.Therefore, Proposition 1 in Koch (2019b) yields that X β has a.s.locally integrable sample paths.

B Case of simple Brown-Resnick fields and β < 1/2
This appendix explains that the results obtained in Sections 2.2 and 3.3 are similar if the Brown-Resnick field considered is simple and the power satisfies β < 1/2.As standard Fréchet margins are rarely encountered in practice, the interest of this section mostly lies in a better understanding of some properties of simple Brown-Resnick fields and in possible applications to inference (using, e.g., the method of moments).

Figure 1 :
Figure 1: The grey and white cells correspond to the 226 and 118 calibration and validation grid points, respectively.

Figure 2 :
Figure 2: Model's performance on the validation grid points.Theoretical pairwise extremal coefficient function from the fitted Brown-Resnick model (red line) and empirical pairwise extremal coefficients (dots).The grey and black dots are pairwise and binned estimates, respectively.The empirical extremal coefficients have been computed using the empirical distribution functions (left) and the obtained GEV parameters (right).

Figure 3 :
Figure 3: Performance of the fitted Brown-Resnick model on the validation grid points.The top row concerns maxima for pairs of validation grid points separated by a low (left), moderate (middle) and long (right) distance.The middle row focuses on minima (left), mean (middle) and maxima (right) for a group of 40 validation grid points chosen randomly.The bottom row concerns minima (left), mean (middle) and maxima (right) for all 118 validation grid points.Overall envelopes at the 95% confidence level are depicted in dark grey.

Figure 4 :
Figure 4: Comparison between observed fields of pointwise maxima and realizations from the fitted Brown-Resnick model.On the left, pointwise maxima over the period October 2005-March 2006 (top) and the period October 2002-March 2003 (bottom).On the right, examples of realizations from the model having values comparable with those in the first column.

Figure 5 :
Figure5: Evolution of D X,10 (x 1 , x 2 ) with respect to x 2 − x 1 for X being the Brown-Resnick field with semivariogram (5) and parameters in Table2.

Figure 6 :
Figure6: Evolution of D X,β (x 2 − x 1 ) with respect to the distance x 2 − x 1 and the power β, where X is the Brown-Resnick field with semivariogram (5) with ψ = 0.5 (top left), 0.81 (top right), 1.5 (bottom left) and 2 (bottom right), and whose other parameters are given in Table2.

Figure 8 :
Figure 8: Evolution of the function I β with respect to the distance h and the power β for β ∈ [−1.6, 0.45].

Table 1 :
CLIC values and parameters' estimates (standard errors inside the brackets) of the Brown-Resnick and Smith models.