Stochastic derivative estimation for max-stable random fields

TL;DR

This paper develops and validates stochastic derivative estimation methods for max-stable random fields, crucial in spatial extremes, with applications demonstrated in risk and dependence measures for insurance.

Contribution

It provides tractable conditions for applying likelihood ratio and infinitesimal perturbation analysis methods to max-stable fields like Brown–Resnick and Smith, which was previously complex.

Findings

Both LRM and IPA perform well in simulations.

Conditions ensure the validity of derivative estimations.

Application to insurance risk measures demonstrates practical utility.

Abstract

We consider expected performances based on max-stable random fields and we are interested in their derivatives with respect to the spatial dependence parameters of those fields. Max-stable fields, such as the Brown--Resnick and Smith fields, are very popular in spatial extremes. We focus on the two most popular unbiased stochastic derivative estimation approaches: the likelihood ratio method (LRM) and the infinitesimal perturbation analysis (IPA). LRM requires the multivariate density of the max-stable field to be explicit, and IPA necessitates the computation of the derivative with respect to the parameters for each simulated value. We propose convenient and tractable conditions ensuring the validity of LRM and IPA in the cases of the Brown--Resnick and Smith field, respectively. Obtaining such conditions is intricate owing to the very structure of max-stable fields. Then we focus on risk and dependence measures, which constitute one of the several frameworks where our theoretical results can be useful. We perform a simulation study which shows that both LRM and IPA perform well in various configurations, and provide a real case study that is valuable for the insurance industry.

Tables

Table 1. Table 1: The left panel gives the values of ∂ R ​ ( 𝜽 ) / ∂ κ | 𝜽 = θ 0 evaluated-at 𝑅 𝜽 𝜅 𝜽 subscript 𝜃 0 \partial R(\mathbf{\boldsymbol{\theta}})/\partial\kappa|_{\boldsymbol{\theta}=\mathbf{\theta}_{0}} and ∂ R ​ ( 𝜽 ) / ∂ ψ | 𝜽 = 𝜽 0 evaluated-at 𝑅 𝜽 𝜓 𝜽 subscript 𝜽 0 \partial R(\boldsymbol{\theta})/\partial\psi|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}} in the various configurations. The middle one displays the previous values normalized by R ​ ( 𝜽 0 ) 𝑅 subscript 𝜽 0 R(\boldsymbol{\theta}_{0}) . The right one gives R ​ ( 𝜽 0 ) 𝑅 subscript 𝜽 0 R(\mathbf{\boldsymbol{\theta}}_{0}) . All values are the theoretical ones.

	$\kappa$		$\psi$		$\kappa$		$\psi$		$R\left({𝜽}_0\right)$
	2	3	2	3	2	3	2	3	2	3
${(1,1)}^{\prime }$	0.048	0.046	0.131	0.126	0.061	0.058	0.167	0.158	0.784	0.797
${(3,2)}^{\prime }$	0.074	0.074	-0.044	-0.044	0.122	0.117	-0.072	-0.070	0.610	0.626
${(9,9)}^{\prime }$	0.087	0.089	-0.439	-0.452	0.306	0.302	-1.552	-1.529	0.283	0.296

Table 2. Table 2: Same as in Table 1 , but for ∂ R ​ ( Σ ) / ∂ σ 11 | Σ = Σ 0 evaluated-at 𝑅 Σ subscript 𝜎 11 Σ subscript Σ 0 \partial R(\Sigma)/\partial\sigma_{11}|_{\Sigma=\Sigma_{0}} , ∂ R ​ ( Σ ) / ∂ σ 12 | Σ = Σ 0 evaluated-at 𝑅 Σ subscript 𝜎 12 Σ subscript Σ 0 \partial R(\Sigma)/\partial\sigma_{12}|_{\Sigma=\Sigma_{0}} and ∂ R ​ ( Σ ) / ∂ σ 22 | Σ = Σ 0 evaluated-at 𝑅 Σ subscript 𝜎 22 Σ subscript Σ 0 \partial R(\Sigma)/\partial\sigma_{22}|_{\Sigma=\Sigma_{0}} .

	${\sigma }_{11}$		${\sigma }_{12}$		${\sigma }_{22}$		${\sigma }_{11}$		${\sigma }_{12}$		${\sigma }_{22}$		$R\left({\Sigma }_0\right)$
	2	3	2	3	2	3	2	3	2	3	2	3	2	3
${(1,1)}^{\prime }$	0.174	0.170	0.06	0.058	0.020	0.020	0.242	0.232	0.083	0.080	0.029	0.027	0.717	0.732
${(3,2)}^{\prime }$	0.233	0.243	0.05	0.053	0.011	0.011	1.669	1.655	0.362	0.359	0.078	0.078	0.139	0.147

Table 3. Table 3: Values of the CLIC and fitted parameters of the Brown–Resnick and Smith max-stable random fields. The values inside the brackets are the standard errors.

Brown–Resnick	CLIC		$\kappa$	$\psi$	$\eta$	$\tau$	$\xi$
	6020169		$3.05\left(0.98\right)$	$0.86\left(0.06\right)$	$26.11\left(0.41\right)$	$2.90\left(0.24\right)$	$-0.11\left(0.03\right)$
Smith	CLIC	${\sigma }_{11}$	${\sigma }_{12}$	${\sigma }_{22}$	$\eta$	$\tau$	$\xi$
	6084169	$0.88\left(0.17\right)$	$0.07\left(0.03\right)$	$2.43\left(0.47\right)$	$26.12\left(0.38\right)$	$2.92\left(0.22\right)$	$-0.10\left(0.01\right)$

Table 4. Table 4: The left panel gives the values of ∂ R ​ ( 𝜽 ) / ∂ κ | 𝜽 = 𝜽 0 evaluated-at 𝑅 𝜽 𝜅 𝜽 subscript 𝜽 0 \partial R(\boldsymbol{\theta})/\partial\kappa|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}} and ∂ R ​ ( 𝜽 ) / ∂ ψ | 𝜽 = 𝜽 0 evaluated-at 𝑅 𝜽 𝜓 𝜽 subscript 𝜽 0 \partial R(\boldsymbol{\theta})/\partial\psi|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}} for the different combinations of sites. The middle one displays the previous values normalized by R ​ ( 𝜽 0 ) 𝑅 subscript 𝜽 0 R(\boldsymbol{\theta}_{0}) . The right one gives R ​ ( 𝜽 0 ) 𝑅 subscript 𝜽 0 R(\boldsymbol{\theta}_{0}) . The sensitivities and normalized sensitivities have been computed using the LRM method with S = 10 6 𝑆 superscript 10 6 S=10^{6} .

${𝐱}_2$	$\kappa$	$\psi$	$\kappa$	$\psi$	$R\left({𝜽}_0\right)$
${(1,1)}^{\prime }$	0.039	0.106	0.046	0.126	0.840
${(3,2)}^{\prime }$	0.068	-0.041	0.100	-0.059	0.685
${(9,9)}^{\prime }$	0.096	-0.486	0.278	-1.409	0.345

Equations

$$\underset{i=1}{\overset{n}{\bigvee }}\frac{X_i-b_n}{a_n}\stackrel{d}{=}X,n=1,2,\dots ,$$

$$\underset{i=1}{\overset{n}{\bigvee }}{n}^{-1}{Y}_i\stackrel{d}{=}Y,n=1,2,\dots ,$$

$$X(\mathbf{x})\overset{d}{=}\left\{\begin{array}[]{ll}\left(\eta(\mathbf{x})-\tau(\mathbf{x})/\xi(\mathbf{x})\right)+\tau(\mathbf{x})Y(\mathbf{x})^{\xi(\mathbf{x})}/\xi(\mathbf{x}),&\quad\xi(\mathbf{x})\neq 0\\ \eta(\mathbf{x})+\tau(\mathbf{x})\log(Y(\mathbf{x})),&\quad\xi(\mathbf{x})=0\end{array}\right..$$

$$Y\stackrel{d}{=}\underset{i=1}{\overset{\infty }{\bigvee }}{U}_i{Z}_i,$$

$$F_{\mathbf{\theta }}\left(\mathbf{y}\right)=\exp \left(-V_{\theta }(\mathbf{y})\right),\mathbf{y}\in (0,\infty {)}^M,$$

$$V_{\mathbf{\theta }}(\mathbf{y})=\mathbb{E}\left[\underset{i=1}{\overset{M}{\bigvee }}\frac{Z_{\mathbf{\theta }}({\mathbf{x}}_i)}{y_i}\right],\mathbf{y}\in (0,\infty {)}^M.$$

$$f_{\mathbf{\theta }}\left(\mathbf{y}\right)=\exp \left(-V_{\theta }(\mathbf{y})\right)\sum _{\pi \in \Pi }{\left(-1\right)}^{|\pi |}\prod _{B\in \pi }\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y}),$$

$$R(\mathbf{\theta })=\mathbb{E}\left[H\left({\mathbf{Y}}_{\mathbf{\theta }}\right)\right],$$

$$\mathbb{E}\left[G\left({\mathbf{X}}_{\mathbf{\theta }}\right)\right],$$

$$R\left(\mathbf{\theta }\right)=\mathbb{E}\left[H\left({\mathbf{Y}}_{\mathbf{\theta }}\right)\right]={\int }_{(0,\infty {)}^M}H\left(\mathbf{y}\right)f_{\mathbf{\theta }}\left(\mathbf{y}\right)\text{d}\mathbf{y}\text{.}$$

$${\frac{\partial R\left(\mathbf{\theta }\right)}{\partial \theta }|}_{\mathbf{\theta }={\mathbf{\theta }}_0}={\int }_{(0,\infty {)}^M}H\left(\mathbf{y}\right)\frac{\partial f_{\mathbf{\theta }}\left(\mathbf{y}\right)}{\partial \mathbf{\theta }}\text{d}\mathbf{y}=\mathbb{E}\left[H\left({\mathbf{Y}}_{{\mathbf{\theta }}_0}\right){\frac{\partial \log f_{\mathbf{\theta }}\left({\mathbf{Y}}_{\mathbf{\theta }}\right)}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }={\theta }_0}\right],$$

$$\frac{\partial \log f_{\mathbf{\theta }}\left(\mathbf{y}\right)}{\partial \mathbf{\theta }}=-\frac{\partial V_{\mathbf{\theta }}(\mathbf{y})}{\partial \mathbf{\theta }}+\frac{\exp \left(-V_{\mathbf{\theta }}(\mathbf{y})\right)}{f_{\mathbf{\theta }}\left(\mathbf{y}\right)}\sum _{\pi \in \Pi }{\left(-1\right)}^{|\pi |}\sum _{B\in \pi }\frac{\partial }{\partial \mathbf{\theta }}\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y})\prod _{B^{\prime }\in \pi ,B^{\prime }\ne B}\frac{{\partial }^{|B^{\prime }|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y}).$$

$${\gamma }_{\mathbf{\theta }}({\mathbf{x}}_1,{\mathbf{x}}_2)={\left(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert /\kappa \right)}^{\psi },{\mathbf{x}}_1,{\mathbf{x}}_2\in {\mathbb{R}}^d,$$

$$V_{\mathbf{\theta }}(\mathbf{y})=\sum _{i=1}^M{y}_i^{-1}{\varphi }_i\left(\mathbf{y},\mathbf{\theta }\right),$$

$${\varphi }_i\left(\mathbf{y},\mathbf{\theta }\right)={\Phi }_{M-1}\left({\lambda }_{\mathbf{\theta }}({\mathbf{x}}_i,{\mathbf{x}}_{-i})+\log \left({\mathbf{y}}_{-i}/y_i\right);{\Omega }_{\mathbf{\theta }}^{(i)}\right).$$

$$B_{{\mathcal{V}}_{{\mathbf{\theta }}_0}}=\underset{i=1,\dots ,M}{\inf }\underset{\mathbf{\theta }\in {\mathcal{V}}_{{\mathbf{\theta }}_0}}{\inf }\underset{\mathbf{y}\in (0,\infty {)}^M}{\inf }\left\{{\varphi }_i\left(\mathbf{y},\mathbf{\theta }\right)-{\varphi }_i\left(\mathbf{y},{\mathbf{\theta }}_0\right)\right\}.$$

$$\mathbb{E}\left[|H\left({\mathbf{Y}}_{{\mathbf{\theta }}_0}\right)|\left(1+\frac{1}{M}\sum _{i=1}^M{Y}_{{\mathbf{\theta }}_0,i}^{-1}\right)\left(\sum _{i=1}^M{Y}_{{\mathbf{\theta }}_0,i}^{-\alpha }\right)\exp \left(-B_{{\mathcal{V}}_{{\mathbf{\theta }}_0}}\sum _{i=1}^M{Y}_{{\mathbf{\theta }}_0,i}^{-1}\right){\Vert {\mathbf{Y}}_{{\mathbf{\theta }}_0}\Vert }^{\alpha }\right]<\infty ,$$

$$\underset{k=1,\dots ,L}{\sup }\underset{1\le i,j\le m}{\sup }\underset{\mathbf{\theta }\in {\mathcal{V}}_{{\mathbf{\theta }}_0}}{\sup }\left|\frac{\partial {\lambda }_{\mathbf{\theta }}({\mathbf{x}}_i,{\mathbf{x}}_j)}{\partial {\mathbf{\theta }}_k}\right|<\infty .$$

$${\frac{\partial R\left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }={\mathbf{\theta }}_0}=\mathbb{E}\left[H\left({\mathbf{Y}}_{{\mathbf{\theta }}_0}\right){\frac{\partial \log f_{\mathbf{\theta }}\left({\mathbf{Y}}_{\mathbf{\theta }}\right)}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }={\mathbf{\theta }}_0}\right].$$

$${\hat{f}}_{\mathbf{\theta }}\left(\mathbf{y}\right)=\exp \left(-V_{\theta }(\mathbf{y})\right)\frac{1}{N}\sum _{i=1}^N\prod _{B\in {\pi }_i}\left(-\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\theta }(\mathbf{y})\right),$$

$$g_{\mathbf{\theta }}(\pi |\mathbf{y})=\frac{\prod _{B\in \pi }\left(-\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y})\right)}{{\sum }_{\pi \in \Pi }\prod _{B\in \pi }\left(-\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y})\right)}\propto \prod _{B\in \pi }\left(-\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y})\right).$$

$$\frac{\partial \log {\hat{f}}_{\mathbf{\theta }}\left(\mathbf{y}\right)}{\partial \mathbf{\theta }}=-\frac{\partial V_{\mathbf{\theta }}(\mathbf{y})}{\partial \mathbf{\theta }}+\frac{\exp \left(-V_{\mathbf{\theta }}(\mathbf{y})\right)}{{\hat{f}}_{\mathbf{\theta }}\left(\mathbf{y}\right)}\frac{1}{N}\sum _{i=1}^N{\left(-1\right)}^{|{\pi }_i|}\sum _{B\in {\pi }_i}\frac{\partial }{\partial \mathbf{\theta }}\frac{{\partial }^{|B|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y})\prod _{B^{\prime }\in {\pi }_i,B^{\prime }\ne B}\frac{{\partial }^{|B^{\prime }|}}{\partial {\mathbf{y}}_B}{V}_{\mathbf{\theta }}(\mathbf{y}).$$

$${\frac{\partial R\left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }={\mathbf{\theta }}_0}=\mathbb{E}\left[{\frac{\partial H\left({\mathbf{Y}}_{\mathbf{\theta }}\right)}{\partial \mathbf{\theta }}|}_{\mathbf{\theta }={\mathbf{\theta }}_0}\right],$$

$$\left|H\left({\mathbf{Y}}_{{\mathbf{\theta }}_1}\right)-H\left({\mathbf{Y}}_{{\mathbf{\theta }}_2}\right)\right|\le \Vert {\mathbf{\theta }}_1-{\mathbf{\theta }}_2\Vert B_{{\mathbf{\theta }}_0}$$

$$\underset{\mathbf{\theta }\in {\mathcal{V}}_{{\mathbf{\theta }}_0}}{\sup }\Vert \frac{\partial H\left({\mathbf{Y}}_{\mathbf{\theta }}\right)}{\partial \mathbf{\theta }}\Vert \le B_{{\mathbf{\theta }}_0},$$

$$Y_{\Sigma }(\mathbf{x})=\underset{i=1}{\overset{\infty }{\bigvee }}{U}_i{\phi }_M(\mathbf{x}-{\mathbf{C}}_i,\Sigma ),\mathbf{x}\in {\mathbb{R}}^d,$$

$${\frac{\partial R\left(\Sigma \right)}{\partial \Sigma }|}_{\Sigma ={\Sigma }_0}={\left({\frac{\partial R\left(\Sigma \right)}{\partial {\sigma }_{ij}}|}_{\Sigma ={\Sigma }_0}\right)}_{ij}.$$

$$\frac{\partial H\left({\mathbf{Y}}_{\Sigma }\right)}{\partial \Sigma }=\sum _{j=1}^M\frac{\partial H\left({\mathbf{Y}}_{\Sigma }\right)}{\partial y_j}\frac{\partial Y_{\Sigma }({\mathbf{x}}_j)}{\partial \Sigma }.$$

$$\underset{\Sigma \in {\mathcal{V}}_{{\Sigma }_0}}{\sup }{\Vert \frac{\partial \log Y_{\Sigma }(\mathbf{x})}{\partial \Sigma }\Vert }^q\le C_{{\Sigma }_0}(\mathbf{x},q)\text{mbox}and\mathbb{E}\left[C_{{\Sigma }_0}(\mathbf{x},q)\right]<\infty .$$

$$\underset{j=1,...,M}{\sup }\mathbb{E}\left[\underset{\Sigma \in {\mathcal{V}}_{{\Sigma }_0}}{\sup }{\left|Y_{\Sigma }({\mathbf{x}}_j)\frac{\partial H\left({\mathbf{Y}}_{\Sigma }\right)}{\partial y_j}\right|}^p\right]<\infty ,$$

Full text

Key words: Infinitesimal Perturbation Analysis; Likelihood Ratio Method; Max-stable random fields; Monte Carlo computation, Risk assessment.

1 Introduction

Sensitivity analysis (SA) is essential in many fields and constitutes an active research area. It consists in the quantitative assessment of how changes in a specific model parameter impact the so-called expected performance. The expected performance is the expectation of a function of an underlying stochastic model, and SA aims at computing its derivatives with respect to the model parameters. If the expected performance is excessively sensitive to some critical parameters, it warns the decision maker not to be overly confident in the obtained value, especially if there is a huge uncertainty on these parameters. It also informs him/her that energy should be invested to find more reliable estimates of these parameters (e.g., by developing new inference methods). By the delta method, the knowledge of the derivative of the expected performance with respect to a parameter of interest allows one to translate confidence bounds for this parameter into confidence bounds for the expected performance. Finally, if the aim is to maximize the expected performance with respect to some parameter, knowing its derivative is often required as most optimization algorithms are based on gradient methods.

To the best of our knowledge, SA was never considered for expected performances based on extreme-value models, although this is of high necessity. Indeed, such models are by essence generally fitted using a small amount of data, resulting in a large uncertainty in the parameter estimates. For a practical introduction about extreme-value theory, see, e.g., Coles, (2001). Extreme-value models are valuable in many applications. Max-stable random fields (e.g., de Haan,, 1984; de Haan and Ferreira,, 2006), which constitute an extension of univariate extreme-value theory to the infinite-dimensional (e.g., spatial) setting, are, for instance, well suited for the modelling of the pointwise maxima of variables having a spatial extent such as environmental ones. Such a modelling is useful for quantifying the impact of natural disasters, which is essential for insurance companies and civil authorities, especially in a context of climate change. However, owing to the complex structure of max-stable fields, their multivariate density is in general not available for more than two or three sites, making estimation highly non-trivial. Commonly used estimators are not asymptotically efficient (in the sense of the Cramér-Rao bound), which, combined with the uncertainty stemming from the scarcity of data mentioned above, can lead to estimates that are far from the true value. It is therefore essential to assess the sensitivity of any expected performance based on max-stable fields with respect to the parameters. The purpose of this paper is precisely to provide tractable conditions enabling such an assessment and to illustrate this in a case where the expected performance is a risk or dependence measure.

We focus on expected performances that do not have any closed-form expression, entailing that assessment of their derivative cannot be achieved using differentiation or finite differences based on the analytical expression, but requires simulation-based approaches. There are essentially three classic stochastic derivative estimation methods. The first one is based on finite-differences. Corresponding estimators involve a bias-variance trade-off and require simulating at multiple parameter values, making this methodology less powerful than the two others. The second approach, referred to as likelihood ratio method (LRM), is based on the derivatives of the density function associated with the simulation model. The third one, called infinitesimal perturbation analysis (IPA), relies on computing the derivative with respect to the parameter of each simulated value (and then averaging them); IPA is thus also referred to as a sample path differentiation method. Contrary to IPA where the parameter is considered as purely structural, in LRM it is viewed as a parameter of the probability measure. LRM requires the density function to be explicit and differentiable and, on the other hand, IPA requires differentiability of the sample paths with respect to the parameter, which is a quite strong condition. It is argued in Glasserman, (2003, Section 7.4) that LRM estimators often have a larger variance than IPA estimators, explaining why IPA is generally considered as the best derivative estimator. A huge literature is dedicated to LRM and IPA; excellent general references are in particular Asmussen and Glynn, (2010), Chapter VII, and Glasserman, (2003), Chapter 7. Both approaches are widely applied to many fields, for instance to finance where many risk hedging strategies involve computing sensitivities of option prices to the underlying assets’ prices and other parameters (the so-called Greeks). Broadie and Glasserman, (1996) and Chen and Fu, (2001) use IPA for option pricing and mortgage-backed securities, respectively. Motivated by financial applications, Glasserman and Liu, (2010) investigate LRM in the case where the relevant densities are only known through their characteristic functions or Laplace transforms.

To our knowledge, SA, and thus LRM and IPA, have not yet been considered for expected performances based on extreme-value models, a fortiori max-stable fields. Among both methods, none can be used for all max-stable fields systematically. LRM is often ruled out when the random performance involves values of the field at more than two or three sites since the multivariate density is not explicitly available. However, there are notable exceptions such as the Brown–Resnick random field (Brown and Resnick,, 1977; Kabluchko et al.,, 2009) or the extremal t random field (Opitz,, 2013). The main technical challenge to apply LRM is to show that interchange between derivative (with respect to the parameter) and integration over the space of possible values of the field at the sites considered is feasible. On the other hand, IPA involves showing that the paths of the random field are differentiable with respect to the spatial dependence parameters, which is often arduous as max-stable fields arise as the pointwise maxima over an infinite number of latent random fields. Nevertheless, when the paths of these latent fields are sufficiently smooth, such as for the Smith random field (Smith,, 1990), it is possible. The Smith field belongs to the class of Brown–Resnick fields, but, being a border case, its properties in terms of availability of a closed-form multivariate density and smoothness of the paths are very different than for most Brown–Resnick fields.

Our main theoretical contribution consists in proposing convenient conditions ensuring the validity of LRM and IPA for estimating the derivatives of expected performances based on the Brown–Resnick field and the Smith field, respectively. The Brown–Resnick field is very flexible and is one of the most appropriate models for environmental data among currently known max-stable models (e.g., Davison et al.,, 2012). We focus on the derivatives with respect to the spatial dependence parameters, since showing the validity of LRM or IPA for estimating the derivatives with respect to the marginal parameters (those of the generalized extreme-value distribution) is easy (under mild assumptions). Our conditions are as tractable as possible; e.g., in the case of IPA for the Smith field, the condition involves the derivative of the expected performance with respect to the field values but not to the spatial dependence parameters. Obtaining such practical conditions is difficult owing to the complex structure of max-stable fields.

Our second contribution pertains to the context of risk assessment in insurance or finance. Our results are insightful for studying the sensitivity of risk or dependence measures based, e.g., on insured losses triggered by extreme events having a spatial extent (typically such as weather events), or the sensitivity of prices of event-linked securities such as catastrophe bonds. Indeed, many risk or dependence measures as well as prices can be written as expected performances. After providing a general discussion, we focus on a specific dependence measure for insured losses due to extreme wind speeds. This measure has analytical derivatives for both the Brown–Resnick and the Smith fields, and hence allows us to compare the LRM and IPA estimates with the true values. We implement LRM and IPA (adapting when necessary existing simulation algorithms) and perform a thorough numerical study which shows that both estimation methods are very accurate. Finally, we propose a concrete application which is valuable for actuarial practice. Using reanalysis wind speed data, we study the sensitivities of the aforementioned dependence measure in an area centred over the Ruhr region in Germany. This application highlights the importance of assessing the sensitivity of risk or dependence measures in real practice.

The remainder of the paper is organized as follows. Section 2 first recalls useful results about max-stable fields and presents the concept of expected performance. Then it introduces LRM and IPA with their associated examples of max-stable fields (Brown–Resnick and Smith, respectively) and states our conditions guaranteeing the validity of both methods in these cases. Section 3 is devoted to risk assessment applications. After a general discussion, we introduce the aforementioned specific dependence measure, we present the simulation study and we finally expose our real case study. Section 4 briefly summarizes our contribution and presents some perspectives. Throughout the paper, we shall use the following notations. Let ′ denote transposition. For $\mathbf{x}=(x_{1},\dots,x_{d})^{\prime}\in\mathbb{R}^{d}$, $\left\|\mathbf{x}\right\|^{2}=\mathbf{x}^{\prime}\mathbf{x}=\sum_{i=1}^{d}x_{i}^{2}$, for a matrix $A\in\mathbb{R}^{d\times d}$, $\left\|A\right\|=\sup\{\left\|A\mathbf{x}\right\|:\mathbf{x}\in\mathbb{R}^{d}$ such that $\left\|\mathbf{x}\right\|=1\}$ and, for a positive definite symmetric matrix $\Sigma\in\mathbb{R}^{d\times d}$, $\left\|\mathbf{x}\right\|_{\Sigma^{-1}}^{2}=\mathbf{x}^{\prime}\Sigma^{-1}\mathbf{x}$. Furthermore, $\bigvee$ denotes the supremum when the latter is taken over a countable set. Finally, $\overset{d}{=}$ denotes equality in distribution. In the case of random fields, by distribution we mean the set of all finite-dimensional distributions. All proofs can be found in the appendix.

2 Stochastic derivative estimators

2.1 Max-stable random fields and associated expected performance

A random field $X$ on $\mathbb{R}^{d}$ with non-degenerate marginals is called max-stable if there are continuous functions $a_{n}(\mathbf{\cdot})>0$ and $b_{n}(\mathbf{\cdot})$ on $\mathbb{R}^{d}$ such that if $X_{1},\ldots,X_{n}$ are independent copies of $X$ then

[TABLE]

(pointwise maxima). A random field $Y$ on $\mathbb{R}^{d}$ with standard Fréchet margins (i.e., $\mathbb{P}(Y(\mathbf{x})\leq y)=\exp(-1/y)$, $y>0$, $\mathbf{x}\in\mathbb{R}^{d}$) is max-stable if

[TABLE]

where $Y_{1},\ldots,Y_{n}$ are independent copies of $Y$; such a field is said to be simple max-stable. It is well-known that there exist continuous functions $\eta(\mathbf{\cdot})$, $\tau(\mathbf{\cdot})>0$ and $\xi(\mathbf{\cdot})$ on $\mathbb{R}^{d}$, called the location, scale and shape functions such that

[TABLE]

Any simple max-stable random field $Y$ on $\mathbb{R}^{d}$ can be written (e.g., de Haan,, 1984) as

[TABLE]

where the $(U_{i})_{i\geq 1}$ are the points of a Poisson point process on $(0,\infty)$ with intensity function $u^{-2}\mathrm{d}u$ and the $Z_{i},i\geq 1$, are independent copies of a non-negative random field $Z$ such that, for all $\mathbf{x}\in\mathbb{R}^{d}$, $\mathbb{E}[Z(\mathbf{x})]=1$. Conversely, any field of the form (2) is simple max-stable. We assume that the distribution of $Z$ depends on a spatial dependence parameter $\boldsymbol{\theta}\in\mathbb{R}^{L}$, and we now denote $Z$ by $Z_{\boldsymbol{\theta}}$ (as well as $Y$ by $Y_{\boldsymbol{\theta}}$). The dependence of $Z_{\boldsymbol{\theta}}$ with respect to $\boldsymbol{\theta}$ will be explicitly described below.

We consider $M$ sites $\mathbf{x}_{1},\ldots,\mathbf{x}_{M}\in\mathbb{R}^{d}$, and let $\mathbf{Y}_{\boldsymbol{\theta}}=\left(Y_{\boldsymbol{\theta}}(\mathbf{x}_{1}),\ldots,Y_{\boldsymbol{\theta}}(\mathbf{x}_{M})\right)^{\prime}$. It is known that the distribution function of $\mathbf{Y}_{\boldsymbol{\theta}}$ is

[TABLE]

where $V_{\boldsymbol{\theta}}$ is the so-called exponent measure function given by

[TABLE]

If the $M$-th order partial derivatives of $V_{\boldsymbol{\theta}}$ exist, then the density function of $\mathbf{Y}_{\mathbf{\theta}}$ may be derived by the Faà di Bruno’s formula for multivariate functions. Let $\mathcal{I}=\{1,\ldots,M\}$. We denote by $\Pi$ the set of all partitions of $\mathcal{I}$ and, for a partition $\pi\in\Pi$, $B\in\pi$ means that $B$ is one of the blocks of the partition $\pi$. Moreover, for any set $B\subset\mathcal{I}$ and $\mathbf{y}\in(0,\infty)^{M}$, we let $\mathbf{y}_{B}=\left(y_{j}\right)_{j\in B}$. The density function of $\mathbf{Y}_{\boldsymbol{\theta}}$ is then given by

[TABLE]

where $|\pi|$ denotes the number of blocks of the partition $\pi$, $|B|$ the cardinality of the set $B$ and $\partial^{|B|}/\partial\mathbf{y}_{B}$ the partial derivative with respect to $\mathbf{y}_{B}$.

The random performance (sometimes also called model output) is defined by $H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)$, where $H$ is a function from $\mathbb{R}^{M}$ to $\mathbb{R}$, and the expected performance is its expectation, i.e.,

[TABLE]

provided that $\mathbb{E}\left[|H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)|\right]<\infty$. Below we give conditions ensuring the applicability of LRM and IPA to estimate $\partial R\left(\boldsymbol{\theta}\right)/\partial\boldsymbol{\theta}$ at a specific $\boldsymbol{\theta}_{0}\in\mathbb{R}^{L}$, for $\mathbf{Y}_{\boldsymbol{\theta}}$ built from the Brown–Resnick and the Smith fields, respectively. The quantity in (4) is pretty general and the results developed in this section can thus be applied for various purposes. One of them, particularly important in finance and insurance, is to estimate the sensitivities of risk and dependence measures; this is done in Section 3.

When fitted to data, a max-stable field is not simple but written as in (1) with location, scale and shape functions not uniformly equal to unity. Hence, denoting by $X_{\boldsymbol{\theta}}$ such a field and letting $\mathbf{X}_{\boldsymbol{\theta}}=\left(X_{\boldsymbol{\theta}}(\mathbf{x}_{1}),\ldots,X_{\boldsymbol{\theta}}(\mathbf{x}_{M})\right)^{\prime}$, the expected performance should be written

[TABLE]

where $G$ is a function from $\mathbb{R}^{M}$ to $\mathbb{R}$ such that $\mathbb{E}\left[|G\left(\mathbf{X}_{\boldsymbol{\theta}}\right)|\right]<\infty$. However, we see by (1) that the representation (4) encompasses the case of (5) by letting $H$ depend on the location, scale and shape parameters of $X_{\boldsymbol{\theta}}$ at the $M$ sites. As technical issues arise for the derivatives with respect to the spatial dependence parameters only, we do not consider differentiation with respect to the marginal parameters, (although this is of practical interest), and focus on (4) for convenience.

2.2 Likelihood ratio method

2.2.1 General methodology

We first introduce LRM to the context of expected performances based on max-stable fields. Let $\mathbf{Y}_{\boldsymbol{\theta}}$ be a simple max-stable random vector, $H$ be a function from $\mathbb{R}^{M}$ to $\mathbb{R}$ such that $\mathbb{E}\left[|H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)|\right]<\infty$, $R\left(\boldsymbol{\theta}\right)=\mathbb{E}[H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)]$, and $\boldsymbol{\theta}_{0}$ be a possible value of the parameter $\boldsymbol{\theta}$. LRM requires $\mathbf{Y}_{\boldsymbol{\theta}}$ to have a density function $f_{\boldsymbol{\theta}}$ that can be differentiated with respect to $\boldsymbol{\theta}$. In that case, the expected performance satisfies

[TABLE]

We assume that there exists some non-random neighbourhood of $\boldsymbol{\theta}_{0}$, $\mathcal{V}_{\boldsymbol{\theta}_{0}}$, such that

• •

for each $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$, $\mathbb{E}\left[|H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)|\right]<\infty$,

• •

for almost all $\mathbf{y}\in(0,\infty)^{M}$, $\partial f_{\boldsymbol{\theta}}\left(\mathbf{y}\right)/\partial\boldsymbol{\theta}$ exists for all $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$,

• •

there is an integrable function $\Psi:\mathbb{R}^{M}\rightarrow\mathbb{R}$ such that $|H\left(\mathbf{y}\right)|\sup_{j=1,...,M}\left|\partial f_{\boldsymbol{\theta}}\left(\mathbf{y}\right)/\partial\theta_{j}\right|\leq\Psi\left(\mathbf{y}\right)$ for all $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$ and almost every $\mathbf{y}\in(0,\infty)^{M}$.

Then, by the dominated convergence theorem, differentiation and integration can be interchanged, giving

[TABLE]

where, provided that, for all $B\subset\mathcal{I}$, $\partial\left(\partial^{|B|}V_{\boldsymbol{\theta}}(\mathbf{y})/\partial\mathbf{y}_{B}\right)/\partial\boldsymbol{\theta}$ exists for all $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$,

[TABLE]

Then LRM consists in computing $\partial R\left(\boldsymbol{\theta}\right)/\partial\boldsymbol{\theta}|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}}$ by estimating the expectation in the right-hand side of (6) by Monte Carlo. Note that the above assumptions are quite usual in the theory of maximum likelihood estimation where $\boldsymbol{\theta}\mapsto\partial\log f_{\boldsymbol{\theta}}\left(\mathbf{y}\right)/\partial\boldsymbol{\theta}$ is called the score function.

2.2.2 The case of the Brown–Resnick random field

We now focus on the case of the Brown–Resnick field. Let $W_{\boldsymbol{\theta}}$ be a centred Gaussian random field on $\mathbb{R}^{d}$ with stationary increments and with semivariogram $\gamma_{\boldsymbol{\theta}}$, and let us define $Z_{\boldsymbol{\theta}}(\mathbf{x})=\exp\left(W_{\boldsymbol{\theta}}(\mathbf{x})-\mathrm{Var}(W_{\boldsymbol{\theta}}(\mathbf{x}))/2\right)$, $\mathbf{x}\in\mathbb{R}^{d}$, where $\mathrm{Var}$ denotes the variance. Then the field $Y_{\boldsymbol{\theta}}$ defined by (2) with that $Z_{\boldsymbol{\theta}}$ is referred to as the Brown–Resnick random field associated with the semivariogram $\gamma_{\boldsymbol{\theta}}$ (Brown and Resnick,, 1977; Kabluchko et al.,, 2009). It is stationary333Throughout the paper, stationarity refers to strict stationarity. and its distribution only depends on the variogram (Kabluchko et al.,, 2009, Theorem 2 and Proposition 11, respectively). The case where $W_{\boldsymbol{\theta}}$ is a fractional Brownian motion leads to the commonly used semivariogram

[TABLE]

where $\kappa>0$ and $\psi\in(0,2)$ are the range and the smoothness parameters, respectively, and $\boldsymbol{\theta}=\left(\psi,\kappa\right)^{\prime}$. Generally, $\psi=2$ is allowed in (7), but we exclude that value here as it makes the multivariate density unavailable under a closed-form expression for $M>d+1$, and thus LRM inapplicable in this case.

For any semivariogram $\gamma_{\boldsymbol{\theta}}$, let $\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})=\left(\gamma_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})/2\right)^{1/2}$, $\mathbf{x}_{i},\mathbf{x}_{j}\in\mathbb{R}^{d}$, $i,j=1,\ldots,M$. Let also $\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{-i})=\left(\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})\right)^{\prime}_{j\neq i}$ , $\log\left(\mathbf{y}_{-i}/y_{i}\right)=\left(\log\left(y_{j}/y_{i}\right)\right)_{j\neq i}^{\prime}$, and $\Omega_{\boldsymbol{\theta}}^{(i)}$ be the matrix with $(j,m)$-th entry $2[\lambda_{\boldsymbol{\theta}}^{2}(\mathbf{x}_{i},\mathbf{x}_{j})+\lambda_{\boldsymbol{\theta}}^{2}(\mathbf{x}_{i},\mathbf{x}_{m})-\lambda_{\mathbf{\theta}}^{2}(\mathbf{x}_{j},\mathbf{x}_{m})]$, $j,m\neq i$. We denote, for $p\in\mathbb{N}_{\ast}$, by $\Phi_{p}\left(\cdot;\Omega\right)$ and $\varphi_{p}\left(\cdot;\Omega\right)$ the $p$-dimensional Gaussian distribution and density functions with covariance matrix $\Omega$, respectively. Provided that the matrices $\Omega_{\boldsymbol{\theta}}^{(i)}$, $i=1,\ldots,M$, are positive definite, then the exponent measure function of the Brown–Resnick random field is given by (e.g., Huser and Davison,, 2013)

[TABLE]

where

[TABLE]

Combined with (3), (8) shows the existence of a closed-form expression for the density $f_{\boldsymbol{\theta}}$.

In the following theorem, we provide convenient conditions ensuring the applicability of LRM for estimating $\partial R\left(\boldsymbol{\theta}\right)/\partial\boldsymbol{\theta}|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}}$.

Theorem 1**.**

For a non-random neighbourhood of $\boldsymbol{\theta}_{0}$, $\mathcal{V}_{\boldsymbol{\theta}_{0}}$, let

[TABLE]

Assume that there exist a non-random neighbourhood $\mathcal{V}_{\mathbf{\theta}_{0}}$ and a constant $\alpha>0$ such that

[TABLE]

where $Y_{\boldsymbol{\theta}_{0},i}=Y_{\boldsymbol{\theta}_{0}}\left(\mathbf{x}_{i}\right)$. Assume moreover that

[TABLE]

Then

[TABLE]

These technical but tractable conditions are pretty natural ones to guarantee the validity of the interchange between differentiation and integration.

Realizations of the term inside the expectation in (12) can be obtained since the Brown–Resnick field can be simulated at a finite number of sites and the score can be computed. The simulation methods for the Brown–Resnick field are either exact (e.g., Dombry et al.,, 2016; Oesting et al.,, 2018) or approximate (Schlather,, 2002, Theorem 4).

Nonetheless, the computation of the density (and thus of the score) may be challenging as the number of terms in the sum equals the $M$-th Bell number, which grows super-exponentially in the dimension $M$. To circumvent this issue, a solution is to approximate (3) by Monte Carlo simulation; i.e., for $N\geq 1$,

[TABLE]

where the partitions $\pi_{1},\ldots,\pi_{N}$ form an ergodic sequence with stationary distribution given by

[TABLE]

Dombry et al., (2013) design a Gibbs sampler to generate approximate simulations $\pi_{1},\ldots,\pi_{N}$ without explicitly computing the constant factor in the denominator of (13). We refer to that paper for more details about the practical implementation. Theoretically, the ergodicity of the resulting Markov chain implies that the precision of the approximation is arbitrarily high as $N\rightarrow\infty$. In practice, the number of iterations of the Gibbs sampler, $N$, is typically much smaller than the cardinality of $\Pi$, but the approximation is reasonably good even for moderate values of $N$ because generally only a few partitions $\pi\in\Pi$ are compatible with the data; see, e.g., Huser et al., (2019). For the Brown–Resnick field associated with the semivariogram (7), they conclude that the Gibbs sampler converges quickly and that about $10\times M$ iterations are enough for the algorithm to converge for a large number of parameter configurations.

Using the Monte Carlo based idea, the approximation of the score function is then given by

[TABLE]

Analytical expressions of $\partial(\partial^{|B|}V_{\boldsymbol{\theta}}(\mathbf{y})/\partial\mathbf{y}_{B})/\partial\boldsymbol{\theta}$ are known for the Brown–Resnick random field (e.g., Dombry et al.,, 2017, Section B.4).

2.3 Infinitesimal perturbation analysis

2.3.1 General methodology

We first introduce IPA to the context of expected performances based on max-stable fields. Let $\mathbf{Y}_{\boldsymbol{\theta}}$ be a simple max-stable random vector, $H$ be a function from $\mathbb{R}^{M}$ to $\mathbb{R}$ such that $\mathbb{E}\left[|H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)|\right]<\infty$, $R\left(\boldsymbol{\theta}\right)=\mathbb{E}[H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)]$, and $\boldsymbol{\theta}_{0}$ be a possible value of the parameter $\boldsymbol{\theta}$. IPA requires the random performance $H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)$ to be differentiable with respect to $\boldsymbol{\theta}$. However, the density $f_{\boldsymbol{\theta}}$ does not need to be explicit, making IPA a potential solution when LRM is invalid. When both LRM and IPA can be applied, IPA may be preferable, although the optimal choice depends on the specific form of the random performance. For stochastic processes and random fields, IPA is also called pathwise derivative estimation because it uses differentiation of sample path functionals.

The essence of IPA is to assess the derivative of interest using

[TABLE]

provided that the derivative and the expectation can be interchanged. The derivative can then be computed through estimation of the right-hand side of (14) by Monte Carlo. A sufficient condition for the interchange to hold in the general case is given, e.g., in Asmussen and Glynn, (2010), Chapter VII, Proposition 2.3.

Proposition 1**.**

Assume that $\boldsymbol{\theta}\mapsto H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)$ is an almost surely (a.s.) differentiable function at $\boldsymbol{\theta}_{0}$ and that a.s. $\boldsymbol{\theta}\mapsto H\left(\mathbf{Y}_{\boldsymbol{\theta}}\right)$ satisfies the Lipschitz condition

[TABLE]

for $\boldsymbol{\theta}_{1}$, $\boldsymbol{\theta}_{2}$ in a non-random neighbourhood of $\boldsymbol{\theta}_{0}$, where $\mathbb{E}\left[B_{\mathbf{\theta}_{0}}\right]<\infty$. Then (14) holds.

If $g:I\rightarrow\mathbb{R}$ is differentiable on an open set $I\subset\mathbb{R}^{d}$, and satisfies $\left\|\partial g(\mathbf{x})/\partial\mathbf{x}\right\|\leq K$ for all $\mathbf{x}$ in $I$, then $g$ is Lipschitz continuous with Lipschitz constant at most $K$ over $I$. Therefore, we immediately deduce that, if there exists a random variable $B_{\boldsymbol{\theta}_{0}}$ satisfying $\mathbb{E}\left[B_{\boldsymbol{\theta}_{0}}\right]<\infty$ and such that a.s.

[TABLE]

where $\mathcal{V}_{\boldsymbol{\theta}_{0}}$ is a non-random neighbourhood of $\boldsymbol{\theta}_{0}$, then (14) holds.

2.3.2 The case of the Smith random field

We illustrate IPA in the case where the random performance is based on the Smith random field (Smith,, 1990), for which a closed-form expression for the density is only known when $M$, the number of stations in $\mathbb{R}^{d}$, is smaller or equal than $d+1$ (Genton et al.,, 2011). Letting $(U_{i},\mathbf{C}_{i})_{i\geq 1}$ be the points of a Poisson point process on $(0,\infty)\times\mathbb{R}^{d}$ with intensity function $u^{-2}\mathrm{d}u\times\mathrm{d}\mathbf{c}$, the Smith random field with covariance matrix $\Sigma=\left(\sigma_{ij}\right)_{ij}$ is defined by

[TABLE]

i.e., by taking $Z_{i}\left(\mathbf{x}\right)=\varphi_{M}(\mathbf{x}-\mathbf{C}_{i},\Sigma)$ in (2). It is stationary and simple max-stable, and corresponds to the Brown–Resnick field associated with the semivariogram $\gamma(\mathbf{x})=\mathbf{x}^{\prime}\Sigma^{-1}\mathbf{x}/2$, $\mathbf{x}\in\mathbb{R}^{d}$ (e.g., Huser and Davison,, 2013). Such a semivariogram leads to the impossibility of characterizing the density for $M>d+1$, as mentioned above. As was the case of $\boldsymbol{\theta}$ for the Brown–Resnick random field, $\Sigma$ completely characterizes the dependence structure of the Smith field. For ease, the vector $\boldsymbol{\theta}$ is replaced by the positive definite matrix $\Sigma$ in the following. The derivative of an expected performance $R\left(\Sigma\right)$ written as in (4) with respect to $\Sigma$ at some positive definite matrix $\Sigma_{0}$ is defined by (e.g., Dwyer,, 1967)

[TABLE]

Note that results concerning differentiation with respect to a scalar or a vector also hold in the case of differentiation with respect to a matrix.

Assume now that $\mathbf{y}\mapsto H\left(\mathbf{y}\right)$ is differentiable. The differentiability of the function $\Sigma\mapsto\mathbf{Y}_{\Sigma}$ in a neighbourhood of $\Sigma_{0}$ will be shown in the proof of Theorem 3 (Appendix A.2). Then the chain rule gives

[TABLE]

We shall prove (Theorem 4 in Appendix A.2) that there exists some non-random neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, such that, for any $q>1$, there exists a random variable $C_{\Sigma_{0}}(\mathbf{x},q)$ satisfying a.s.

[TABLE]

This technical outcome will allow us to derive our main result.

Theorem 2**.**

Assume that $\mathbf{y}\mapsto H\left(\mathbf{y}\right)$ is a differentiable function and that there exists $p>1$ such that

[TABLE]

where $\mathcal{V}_{\Sigma_{0}}$ is a non-random neighbourhood of $\Sigma_{0}$. Then

[TABLE]

This non-trivial theorem provides a sufficient condition to use IPA to compute the derivatives of $R\left(\Sigma\right)$. This condition is much more tractable and easier to check than that in Proposition 1. The simplification stems from the fact that we take care in Theorems 3 and 4 of the intricate term of (17), $\partial Y_{\Sigma}(\mathbf{x}_{j})/\partial\Sigma$, which involves the sample path properties with respect to differentiation of the Smith random field. Theorems 3 and 4 are delicate to establish precisely due to the inherent structure of max-stable fields.

In practice we have to simulate realizations of the random matrix inside the expectation in (19), which can be done by simulating the Smith field and using (17). The expression of $\partial Y_{\Sigma}(\mathbf{x}_{j})/\partial\Sigma$ appearing in (17) is given in (29) (Appendix A.2). As the Brown–Resnick field, the Smith field can be simulated exactly (e.g., Oesting et al.,, 2018; Dombry et al.,, 2016) or approximately (Schlather,, 2002, Theorem 4); the latter approach is very accurate.

3 Application to risk assessment

We now focus on one framework (among several others) where the results of Section 2 are useful, which is the context of risk and dependence measures. After detailing the link with (4), we consider a dependence measure which is particularly suited to insurance of damage triggered by extreme wind speeds. We show that the conditions of Theorems 1 and 2 hold in that case and confirm through a simulation study that both LRM and IPA perform very well. Finally, we consider concrete data and show that the sensitivity of our dependence measure can be very high, highlighting the practical importance of studying the sensitivity of functions of max-stable fields in the context of risk assessment.

3.1 Risk measures based on max-stable fields

Here we show that the quantity defined in (4) encompasses many risk and dependence measures, and our focus is mainly on actuarial applications. A univariate risk measure is a mapping from a set of random variables to the real numbers. A dependence measure summarises the strength of dependence between several elements of such a set of random variables. In finance, these random variables often represent portfolio returns, and in an insurance context, they might be the claims associated with insurance policies. When claims are triggered by environmental events, a possible model for the insured cost field is (Koch,, 2017, Section 2.3)

[TABLE]

where $E$ is the (deterministic) insured exposure (i.e., insured value) field, $D_{\mathbf{x}}$ the damage function at site $\mathbf{x}$ and $X$ the random field of the environmental variable generating risk. The application of the damage function $D_{\mathbf{x}}$ to $X(\mathbf{x})$ yields the insured cost ratio (i.e., the insured cost divided by the insured value) at site $\mathbf{x}$, which, multiplied by the insured exposure, gives the corresponding insured cost.

Among the risk measures that can be written as in (4), one can find many sophisticated examples, e.g., in insurance/reinsurance pricing or regulation. For $j=1,\dots,M$, let $C_{j}$ denote the claim of an insurance company at $\mathbf{x}_{j}$ and assume that $C_{j}$ can be written as a function of the max-stable field (e.g., as in the right-hand side of (20)). Premium loadings that are proportional to specific moments of the sum of the $C_{j}$ at two or more sites constitute excellent examples in insurance pricing. In reinsurance, the premium is sometimes based on order statistics of the claims, as in the case of the “excédent du coût moyen relatif” (ECOMOR) or largest claims reinsurance (LCR) treaties. Let us consider, e.g., $M\geq 3$ sites and assume that each of those is associated with an insurance policy whose corresponding claim is $C_{j}$ as above. We consider the ordered values of those claims, $C^{(1:M)}\geq C^{(2:M)}\geq\cdots\geq C^{(M:M)}$. For instance, the risk measure $\mathbb{E}[(C^{(1:M)}-C^{(3:M)})+(C^{(2:M)}-C^{(3:M)})]$ would be involved in the pricing of an ECOMOR reinsurance treaty having the third largest claim as priority. These quantities can be written under the form (4) but behave in a non-linear way and do not have any analytical expression. A valuable example of dependence measure, which will be considered until the end, is presented in the next section.

3.1.1 A specific dependence measure for wind speed

We present in this section a dependence measure for the costs due to high wind speeds that can be written as in (4). We model the cost by (20) in the case $d=2$, where the field of wind speed extremes, $X_{\boldsymbol{\theta}}$, is assumed to be max-stable. Moreover, for any site $\mathbf{x}\in\mathbb{R}^{2}$, we choose as damage function $D_{\mathbf{x}}(x)=(x/u)^{\beta(\mathbf{x})}$, $x\leq u$, where $u>0$ and $\beta(\mathbf{x})\in\mathbb{N}_{\ast}$, which is utterly appropriate in the case of wind. Indeed, as the force exerted by the wind and the corresponding rate of work are proportional to the second and third powers of wind speed, respectively, the total cost for a specific structure is expected to increase as the square or the cube of the maximum wind speed. For studies 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 for the cube, see Lamb and Frydendahl, (1991, Chapter 2, p.7), Emanuel, (2005) and Kantha, (2008). However, in the case of insured costs, several authors have recently found power-laws with much higher exponents; e.g., Prahl et al., (2012) obtained exponents spanning from $8$ to $12$ for insured losses on residential buildings in Germany. Indeed, the presence of a deductible in the insurance contract increases the exponent from $2$ or $3$ to a larger value depending on the deductible (e.g., Prahl et al.,, 2015; Koch,, 2019). The level $u$ corresponds to that value of wind speed which triggers an insured cost ratio equal to unity, and is thus assumed to be larger than the wind speeds observed in practice. For a more detailed review of wind damage functions, see, e.g., Koch, (2019).

We consider two sites ($M=2$), $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{2}$. Recall that the link between $X_{\boldsymbol{\theta}}$ and the associated simple max-stable field $Y_{\boldsymbol{\theta}}$ is given by (1), and let $\eta_{i}=\eta(\mathbf{x}_{i}),\tau_{i}=\tau(\mathbf{x}_{i}),\xi_{i}=\xi(\mathbf{x}_{i})$ and $\beta_{i}=\beta(\mathbf{x}_{i})$, $i=1,2$, such that $\beta_{i}\xi_{i}<1/2$. The shape parameter for wind speed maxima is often slightly negative, implying that the event $X_{\boldsymbol{\theta}}(\mathbf{x}_{i})<0$ occurs with probably not exactly [math] (although extremely close owing to the values of the location and scale parameters). Anyway, this is not a problem as $\beta_{i}\in\mathbb{N}_{\ast}$. Moreover, we take $Y_{\boldsymbol{\theta}}$ to be the Brown–Resnick or the Smith random field and we consider as dependence measure the correlation between the costs due to extreme winds at the two sites, i.e.,

[TABLE]

where $\boldsymbol{\theta}$ must be replaced by $\Sigma$ when the Smith field is considered. The condition $\beta_{i}\xi_{i}<1/2$, $i=1,2$, ensures the existence of the correlation in (21). Correlation is used a lot by practitioners in the finance/insurance industry and the measure $R(\boldsymbol{\theta})$ is of practical interest for any insurance/reinsurance company handling the risk of damage caused by extreme wind speeds; among others, it provides insight about potential spatial diversification. For details, see Koch, (2019) where (21) is thoroughly studied.

We now explain why the measure in (21) is of the form (4) with $M=2$. It is well-known that, if $Y$ is a random variable following the standard Fréchet distribution, then $\mathbb{E}[Y^{\tilde{\beta}}]=\Gamma(1-\tilde{\beta})$ for any $\tilde{\beta}<1$, where $\Gamma$ denotes the gamma function. We now also assume that $\xi_{i}\neq 0$ for $i=1,2$. Consequently, using (1) and the binomial theorem, we obtain, for $i=1,2$,

[TABLE]

where

[TABLE]

Moreover, Corollary 1 in Koch, (2019) gives

[TABLE]

where

[TABLE]

with

[TABLE]

Therefore, (22) and (23) yield

[TABLE]

with

[TABLE]

where $x_{i}\left(y_{i}\right)=\left(\eta_{i}-\tau_{i}/\xi_{i}\right)+\tau_{i}y_{i}^{\xi_{i}}/\xi_{i}$ since $\xi_{i}\neq 0$, $i=1,2$.

On top of being useful for actuarial applications, the measure (21) has a closed-form derivative with respect to the dependence parameters of the max-stable field, allowing us to compare the values of $\partial R(\boldsymbol{\theta})/\partial\boldsymbol{\theta}|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}}$ obtained using LRM or IPA with their true values. The analytical formulas of these derivatives are given in Appendix B.1.

3.1.2 Validity of the assumptions

We now show that the assumptions of Theorems 1 and 2 are valid in the case of (21), for $Y$ being the Brown–Resnick and the Smith random field, respectively. Consequently, LRM and IPA can be used to assess the derivative of (21) with respect to $\boldsymbol{\theta}$ and $\Sigma$, respectively.

Proposition 2**.**

Let $H$ be defined as in (25).

i) Assume that $Y_{\boldsymbol{\theta}}$ is the Brown–Resnick random field with semivariogram given by (7). There exist a non-random neighbourhood $\mathcal{V}_{\boldsymbol{\theta}_{0}}$ of $\boldsymbol{\theta}_{0}$, $\mathcal{V}_{\boldsymbol{\theta}_{0}}$, and a constant $\alpha>0$ such that (10) and (11) hold true.

ii) Assume that $Y_{\Sigma}$ is the Smith random field. Then there exist $p>1$ and a non-random neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, such that (18) holds true.

3.2 Simulation study

We numerically assess the accuracy of the LRM and IPA for computing the sensitivities of the dependence measure $R\left(\boldsymbol{\theta}\right)$ in (21), where $Y$ is the Brown–Resnick field associated with the semi-variogram (7), and the Smith field, respectively. The number of simulations used to approach the expectations in (12) and (19) is denoted by $S$. We display the boxplots of the relative errors of $100$ estimates in different configurations. We take $\beta_{1}=\beta_{2}=\beta$ with $\beta=2,3$, $S=10^{4},10^{5}$, and we look at different combinations of sites $\mathbf{x}_{1},\mathbf{x}_{2}$. We recall that the relative errors can be calculated since the sensitivities of $R\left(\boldsymbol{\theta}\right)$ can be obtained analytically with an integral form. The integrals in (53) and (54) were computed using adaptive quadrature with a relative tolerance of $10^{-7}$ to allow an accurate approximation.

3.2.1 LRM for the Brown–Resnick field

We examine three combinations of sites: we set $\mathbf{x}_{1}=(0,0)^{\prime}$ and $\mathbf{x}_{2}=(1,1)^{\prime},(3,2)^{\prime},(9,9)^{\prime}$. These have been chosen in order to cover a wide range of sensitivities and relative sensitivities. Moreover, we take $\boldsymbol{\theta}_{0}=\left(\psi_{0},\kappa_{0}\right)^{\prime}=(3.05,0.86)^{\prime}$, $\eta_{1}=\eta_{2}=26.11$, $\tau_{1}=\tau_{2}=2.90$ and $\xi_{1}=\xi_{2}=-0.11$, which are the estimates obtained on the data used in the application below (Section 3.3). We simulated the Brown–Resnick field using the rmaxstab function of the SpatialExtremes R package by Ribatet et al., (2018).

Table 1 shows that, for $\kappa$ and $\psi$, the theoretical values of $R\left(\boldsymbol{\theta}_{0}\right)$, its sensitivity and its relative sensitivity do not evolve much when increasing $\beta$ from $2$ to $3$; the absolute values of the relative sensitivities slightly decrease whereas $R\left(\boldsymbol{\theta}_{0}\right)$ weakly increases. The absolute values of the sensitivities and of their relative counterparts are the highest for $\mathbf{x}_{2}=(9,9)^{\prime}$; in that case, the relative sensitivities are very high (about $30\%$ for $\kappa$ and more than $150\%$ for $\psi$), probably owing to the fairly low value of $R\left(\boldsymbol{\theta}_{0}\right)$. Such large values typically explain the practical importance of properly assessing sensitivities, e.g., in an insurance context.

As $M=2$, the density of the Brown–Resnick field is explicitly known and thus the score function has a closed-form expression (available upon request) with a small number of terms; it is therefore needless to apply the MCMC methodology through Gibbs sampler described above.

It is important to remark that $\partial\log f_{\boldsymbol{\theta}}(y_{1},y_{2})/\partial\psi$ and $\partial\log f_{\boldsymbol{\theta}}(y_{1},y_{2})/\partial\kappa$ are proportional (Appendix B.2), which entails that the LRM estimates of $\partial R(\boldsymbol{\theta})/\partial\psi$ and $\partial R(\boldsymbol{\theta})/\partial\kappa$ are also proportional. Moreover, by (54) and (55), the true values of those derivatives are proportional with the same factor, implying that the relative errors of the estimates of $\partial R(\boldsymbol{\theta})/\partial\psi$ and $\partial R(\boldsymbol{\theta})/\partial\kappa$ are equal. Consequently, we only display the results for $\kappa$.

Figure 1 shows that the LRM estimator is unbiased, as expected. Moreover its variability substantially decreases when increasing $S$ from $10^{4}$ to $10^{5}$ and becomes very small (the relative error of most estimates is less than $5\%$). This variability is the lowest for $\mathbf{x}_{2}=(9,9)^{\prime}$, perhaps due to the fact that this configuration corresponds by far to the largest sensitivities and relative sensitivities. To some extent, we may expect the accuracy of the estimation to be an increasing function of the absolute value of the relative sensitivity. Nevertheless, this seems to be more complex: the variability is lower for $\mathbf{x}_{2}=(3,2)^{\prime}$ than for $\mathbf{x}_{2}=(1,1)^{\prime}$ and, for all site combinations, for $\beta=3$ than for $\beta=2$. Overall, the LRM estimator is very satisfying in all configurations, whether the sensitivity (or relative sensitivity) is low or large.

3.2.2 IPA for the Smith field

We consider two combinations of sites: we set $\mathbf{x}_{1}=(0,0)^{\prime}$ and $\mathbf{x}_{2}=(1,1)^{\prime},(3,2)^{\prime}$. Additionally, we take

[TABLE]

$\eta_{1}=\eta_{2}=26.12$, $\tau_{1}=\tau_{2}=2.92$ and $\xi_{1}=\xi_{2}=-0.10$, which are the estimates obtained in Section 3.3. Unlike in Section 3.2.1, we do not consider the case $\mathbf{x}_{2}=(9,9)^{\prime}$ as $R(\Sigma_{0})$ is approximately equal to [math] in that configuration. Note that we do not display the results for $\sigma_{21}$ as they are exactly the same as those for $\sigma_{12}$, consistently with the theory. The analytical computation of the terms $\partial Y_{\Sigma}(\mathbf{x}_{j})/\partial\Sigma$, $j=1,2$, which is necessary to implement IPA (see (17)), requires the coordinates of the centers of the “storms” (see Smith,, 1990, for the interpretation of the Smith field in terms of storms) realizing the maxima at the sites $\mathbf{x}_{j}$ (see (29) in Section A.2). To the best of our knowledge, these quantities cannot be obtained from the simulation algorithms available on the Web (e.g., in R packages like SpatialExtremes by Ribatet et al., (2018) or in the code by Dombry et al., (2016) available on the Biometrika website). To overcome this impediment and for other technical reasons, we programmed the simulation algorithm of the Smith random field ourselves by adapting the approach of Schlather, (2002, Theorem 4). Regarding the quantity $r$ appearing in that approach, we took the value $r=15$ in order to ensure an accurate simulation. The corresponding code will be available.

Table 2 shows that, when increasing $\beta$ from $2$ to $3$, the sensitivities and relative sensitivities slightly decrease whereas $R(\Sigma_{0})$ weakly increases. The relative sensitivities are the highest for $\mathbf{x}_{2}=(3,2)^{\prime}$; they take very high values for $\sigma_{11}$ and $\sigma_{12}$ (about $36\%$ for $\sigma_{12}$ and more than $160\%$ for $\sigma_{11}$), probably due to the low value of $R(\Sigma_{0})$.

Figures 2–4 indicate that the IPA estimator is unbiased (at least for $S$ large enough) and, as expected, the variability decreases when increasing $S$ from $10^{4}$ to $10^{5}$. In most cases, it reaches a low level (with a relative error of most estimates less than $5\%$). The variability is slightly lower in the case of $\mathbf{x}_{2}=(3,2)^{\prime}$, which can be explained by the fact that the relative sensitivity is much larger in that case. It is however the converse for $\sigma_{22}$, coefficient for which the increase of relative sensitivity compared to the case $\mathbf{x}_{2}=(1,1)^{\prime}$ is lower than for the other coefficients. Moreover, the variability is systematically higher for $\sigma_{22}$ than for the other coefficients, and, whatever the coefficient and the combination of sites, it tends to be lower for $\beta=3$ than for $\beta=2$. We intuitively expect the method to perform the best for high relative sensitivities, but we see that it seems more complex. Overall, IPA performs very well on this example, whether the sensitivity (or relative sensitivity) is low or high.

3.3 Application

This section is devoted to the computation of sensitivities of the measure $R(\boldsymbol{\theta})$ (see (21)) in a real case study that is valuable for the insurance/reinsurance industry. We consider concrete wind speed data on a region where insurance/reinsurance plays a significant role and for which power-law type damage functions for losses stemming from extreme wind speeds are clearly documented in the literature.

3.3.1 Data

We consider publicly available data from the European Centre for Medium Range Weather Forecasting (ECMWF), more precisely a subset of the ERA5 (ECMWF Reanalysis 5th Generation) data set. More specifically the data we study consist in hourly $10$ m wind gust time series from 1 January 1979 at 07:00 to 1 June 2019 at 00:00. The region we consider is a rectangle from $6^{\circ}$ to $9.75^{\circ}$ longitude and $49.75^{\circ}$ to $52.25^{\circ}$ latitude, and the resolution is $0.25^{\circ}$ longitude and $0.25^{\circ}$ latitude; see Figure 5. Thus, the rectangle contains 176 grid points and is basically centred over the Ruhr region in Germany. This area exhibits a high total residential insured value per unit of surface, and, according to Prahl et al., (2012, Figure 2), a damage function of power-law type with an exponent around 8 seems appropriate for insured losses (due to high wind speeds) on residential buildings in that region. At each grid point, we compute the $41$ seasonal (from October to March) maxima, leading to the dataset we use to fit different max-stable models for $Y_{\boldsymbol{\theta}}$. For the first season, the maximum is computed over January–March. Considering the period October–March allows us to remove seasonal non-stationarity in the wind speed time series and to focus mainly on winter storms (and not on highly localized extreme winds occurring during intense summer thunderstorms).

3.3.2 Results

The max-stable models we consider are both the Brown–Resnick random field with semi-variogram specified in (7) and the Smith field. Regarding the location, scale and shape parameters, it is reasonable to consider them as constant over the region (shown in a more detailed analysis of the same data that will be soon available in a subsequent paper). Modelling these parameters using trend surfaces (instead of considering specific parameters at each grid point) is common as it reduces the estimation time and allows prediction at sites where no observations are available.

Both models are fitted using maximum pairwise likelihood estimation (e.g., Padoan et al.,, 2010; Davison et al.,, 2012) implemented in the fitmaxstab function of SpatialExtremes (Ribatet et al.,, 2018). The marginal and dependence parameters are estimated in a single step. Then we perform model selection by minimization of the composite likelihood information criterion (CLIC); see Varin and Vidoni, (2005).

Table 3 shows that, according to the CLIC, the Brown–Resnick random field is more compatible with our data than the Smith random field is. We also fitted the Schlather random field (Schlather,, 2002) with various correlation functions (Whittle-Matérn, powered exponential and Cauchy) and the Brown–Resnick random field is more appropriate according to the CLIC. This is in agreement, e.g., with the results obtained by Davison et al., (2012) in the case of rainfall.

Figure 6 shows that the theoretical extremal coefficient function (e.g., Schlather and Tawn,, 2003) of the fitted Brown–Resnick model agrees reasonably well with the empirical extremal coefficients, though it is slightly above their binned estimates. This low underestimation of the spatial dependence may come from the choice of very parsimonious trend surfaces for the location, scale and shape parameters. The comparison done in Figure 6 constitutes a classical graphical goodness-of-fit diagnostic for max-stable models and it shows here that the proposed model fits the data sufficiently well for the purpose of this application.

The final dependence measure is

[TABLE]

where $X_{\boldsymbol{\theta}}$ is a general Brown–Resnick field having as marginal parameters those in Table 3, and we aim at providing its derivative at $\boldsymbol{\theta}_{0}=(3.05,0.86)^{\prime}$. As the Brown–Resnick random field appears as a suitable model for wind speed maxima and the power-law with $8$ as exponent is an appropriate damage function, this dependence measure is well suited in practice for (re)insurance companies on the region considered. The values of the sensitivities obtained by LRM and their relative counterparts are given in Table 4 for the same combinations of sites as in Section 3.2.1. The combinations $\mathbf{x}_{1}=(0,0)^{\prime},\mathbf{x}_{2}=(1,1)^{\prime}$ and $\mathbf{x}_{1}=(0,0)^{\prime},\mathbf{x}_{2}=(3,2)^{\prime}$ correspond to pairs of sites within our region of interest if we put the origin $(0,0)^{\prime}$ on the left part and in the lower-left corner of the rectangle, respectively; by stationarity of the field, the origin can be chosen arbitrarily. The site $\mathbf{x}_{2}=(9,9)^{\prime}$ is far to the East of our region even we take as origin the lower-left corner of the rectangle, but we chose it to emphasize that sensitivities can be very high for sites that are highly distant. The relative sensitivities are slightly lower in absolute value than in the cases $\beta=2$ and $\beta=3$ (see Table 1). However they can be quite substantial (e.g., for $\psi$ in the case $\mathbf{x}_{2}=(1,1)^{\prime}$ and $\kappa$ in the case $\mathbf{x}_{2}=(3,2)^{\prime}$) or very high (for both $\kappa$ and $\psi$ in the case $\mathbf{x}_{2}=(9,9)^{\prime}$), which highlights that looking at sensitivities is strongly recommended in concrete risk assessment studies.

4 Discussion

Max-stable random fields are particularly suitable to model extreme spatial (e.g., environmental) events and, if applicable, LRM as well as IPA are powerful techniques to estimate the derivatives of an expected performance with respect to the model parameters. In this paper, we introduce these two methods to extreme-value theory by applying them to expected performances based on the Brown–Resnick and the Smith max-stable random fields. In a first part, we give convenient and tractable conditions on these fields that ensure the validity of LRM and IPA. Obtaining such tractable conditions is non-trivial due to the complex structure of max-stable fields. Then, we focus on one of the several frameworks where our theoretical results can be valuable, more precisely the context of risk assessment. We apply LRM and IPA to a specific example of dependence measure which is suited to assess the dependence between damage due to extreme wind speeds at two stations, and is thus worthwhile for actuarial practice. We show through a simulation study that both methods perform very well in various configurations. We finally present a concrete application involving reanalysis wind speed data; it emphasizes the importance of considering sensitivities in a risk assessment context when using max-stable fields.

It would be interesting to see how the LRM and IPA estimators behave in other configurations. In some situations (e.g., where the derivatives are very close to [math]), appropriate improvements of the classical Monte Carlo method (such as importance sampling) might be considered to reduce the variability of the estimators. Future appealing work would also consist in extending our conditions to other max-stable random fields and possibly even proposing other methods for max-stable fields (such as, e.g., the one by Peng et al.,, 2018) for which LRM and IPA are not applicable.

Acknowledgements

Both authors gratefully acknowledge Anthony C. Davison and anonymous referees for insightful comments. Erwan Koch would like to thank the Swiss National Science Foundation (project number $200021\_178824$) for financial support.

Appendix A Proofs of the main results

A.1 Proof of Theorem 1

Proof.

From Section 4.4 in Dombry et al., (2017), we know that $\mathbf{Y}_{\boldsymbol{\theta}}$ has an Hüsler-Reiss distribution with a strictly conditionally negative definite matrix given by $\Lambda_{\boldsymbol{\theta}}=\left(\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})\right)_{1\leq i,j\leq m}$. By Lemma B.4. in Dombry et al., (2017), we deduce that, for any neighbourhood $\mathcal{V}_{\boldsymbol{\theta}_{0}}$ of $\boldsymbol{\theta}_{0}$ and any $\alpha>0$, there exists a constant $A_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}>0$ such that, for all $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$, $\mathbf{y}\in(0,\infty)^{M}$ with $\left\|\mathbf{y}\right\|=1$ and $B\subset\mathcal{I}$,

[TABLE]

where

[TABLE]

Equation (14) in Dombry et al., (2017) then gives, for all $\mathbf{y}\in(0,\infty)^{M}$,

[TABLE]

Now let us consider the ratio

[TABLE]

By (9), we have

[TABLE]

where $B_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}$ is given in (9). From Equation (28) in Asadi et al., (2015), we know that for all $B\subset\mathcal{I}$ and $\mathbf{y}\in(0,\infty)^{M}$

[TABLE]

with

[TABLE]

Let $\boldsymbol{\Theta}^{\left(\mathbf{y}\right)}=\mathbf{y}/\left\|\mathbf{y}\right\|$. It follows that, for $\pi\in\Pi$,

[TABLE]

with

[TABLE]

Let

[TABLE]

We can note that $C\left(\pi,\boldsymbol{\Theta}^{\left(\mathbf{y}\right)},\boldsymbol{\theta}\right)$ are positive and uniformly bounded for any partition $\pi$, any $\boldsymbol{\Theta}^{\left(\mathbf{y}\right)}\in\mathcal{S}_{+}=\{\mathbf{y}\in(0,\infty)^{M}:\left\|\mathbf{y}\right\|=1\}$, and any $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$. Then

[TABLE]

The numerator and the denominator are polynomials of order $M$ in $\left\|\mathbf{y}\right\|^{-1}$. If $|\pi|=1$, then $\pi$ $=\{\{1,...,M\}\}$ and

[TABLE]

If $|\pi|=M$, then $\pi$ $=\{\{1\},...,\{M\}\}$ and

[TABLE]

The ratio of both polynomials is bounded over any compact set included in $(0,\infty)$ and has limits as $\left\|\mathbf{y}\right\|$ tends to [math] or $\infty$. Moreover these limits are also bounded for any $\boldsymbol{\Theta}^{\left(\mathbf{y}\right)}\in\mathcal{S}_{+}$. We can deduce that there exists a constant $C_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}$ such that

[TABLE]

It follows from (27) that

[TABLE]

Let us choose

[TABLE]

to get

[TABLE]

for all $\boldsymbol{\theta}\in\mathcal{V}_{\boldsymbol{\theta}_{0}}$ and almost every $\mathbf{y}\in(0,\infty)^{M}$. The result follows by the dominated convergence theorem.

∎

A.2 Proof of Theorem 2

Recall that $(U_{i},\mathbf{C}_{i})_{i\geq 1}$ are the points of a Poisson point process on $(0,\infty)\times\mathbb{R}^{d}$ with intensity function $u^{-2}\mathrm{d}u\times\mathrm{d}\mathbf{c}$. For $i\geq 1$, let $\varphi_{i,\Sigma}$ be the function from $\mathbb{R}^{d}$ to $\mathbb{R}$ defined by $\varphi_{i,\Sigma}\left(\mathbf{x}\right)=U_{i}\varphi_{M}(\mathbf{x}-\mathbf{C}_{i},\Sigma)$.

Let us begin by noting that for each $j=1,\ldots,M$, the supremum $Y_{\Sigma}(\mathbf{x}_{j})=\bigvee_{i=1}^{\infty}U_{i}\varphi_{M}(\mathbf{x}_{j}-\mathbf{C}_{i},\Sigma)$ is a.s. attained by a unique function $\varphi_{i,\Sigma}$ at $\mathbf{x}_{j}$.

Proposition 3**.**

Let $\mathbf{x}_{1},\dots,\mathbf{x}_{M}\in\mathbb{R}^{d}$, and define for $\mathbf{x}\in\mathbb{R}^{d}$,

[TABLE]

Then a.s., we have $|\mathcal{I}_{\mathbf{x}_{j}}|=1$ for all $j\geq 1$, where we recall that $|.|$ stands for the cardinality of a set.

Proof.

Note that for every $\mathbf{x}\in\mathbb{R}^{d}$, $\left\{\varphi_{i,\Sigma}\left(\mathbf{x}\right),i=1,2,\ldots\right\}$ is a Poisson random measure on $(0,\infty)$, which has no atoms. Therefore, the maxima will have no ties with probability one.

∎

Remark 1**.**

The proof is also a direct consequence of Proposition 2.5 in Dombry and Eyi-Minko, (2013).

By Proposition 3, we can define, for $j=1,\ldots,M$, the a.s. unique indexes $i_{\mathbf{x}_{j},\Sigma}$ satisfying

[TABLE]

We now fix $\omega\in\Omega$ and consider one realization of the point process $(U_{i},\mathbf{C}_{i})_{i\geq 1}$, denoted by $(U_{i}\left(\omega\right),\mathbf{C}_{i}\left(\omega\right))_{i\geq 1}$. For ease of exposition, we do not mention $\omega$ further below. Let us consider the set

[TABLE]

i.e., the set of sites $\mathbf{x}\in\mathbb{R}^{d}$ for which the maximum $\bigvee_{i=1}^{\infty}U_{i}\varphi_{M}(\mathbf{x}-\mathbf{C}_{i},\Sigma)$ is attained by at least two distinct functions $\varphi_{j,\Sigma}$ and $\varphi_{k,\Sigma}$.

Remark 2**.**

If $\Sigma=Id_{d}$ (identity matrix with dimension $d\times d$), the set $\mathcal{B}$ is the set of boundaries of the cells of a Poisson Laguerre Tessellation (e.g., Dombry and Kabluchko,, 2018).

A key result is that $\mathcal{B}$ has a null Lebesgue measure.

Proposition 4**.**

The set $\mathcal{B}$ has zero Lebesgue measure with probability one.

Proof.

A point $\mathbf{x}\in\mathcal{B}$ is characterized as follows:

[TABLE]

Let $\mathbf{h}\in\mathbb{R}^{d}$ such that $\mathbf{x}+\mathbf{h}$ is in a neighbourhood of $\mathbf{x}$ but still belongs to $\mathcal{B}$, i.e., such that

[TABLE]

Using (28), we obtain that the previous equality is equivalent to

[TABLE]

Therefore, $\mathbf{x}+\mathbf{h}\in\mathcal{B}$ in a neighbourhood of $\mathbf{x}$ implies that $\mathbf{h}$ is orthogonal to the vector $(\mathbf{C}_{j}-\mathbf{C}_{k})$ for the inner product induced by $\Sigma^{-1}$. Thus, only one direction is suitable for $\mathbf{h}$. $\mathcal{B}$ is an union of segments in $\mathbb{R}^{d}$. Moreover there is no ball around $\mathbf{x}$ belonging to $\mathcal{B}$, which implies that the interior of $\mathcal{B}$ is empty.

Let

[TABLE]

It is easily seen that

[TABLE]

and hence that $\lim_{j\rightarrow\infty}\nu(\mathcal{B}_{j,k})=0$, where $\nu$ denotes the Lebesgue measure in $\mathbb{R}^{d}$. Therefore,

[TABLE]

∎

Let $\mathbf{x}\in\mathbb{R}^{d}/\mathcal{B}$. We are now interested in the existence of a neighbourhood of $\Sigma_{0}$ over which $i_{\mathbf{x},\Sigma}$ is constant and thus $\Sigma\mapsto Z_{\Sigma}(\mathbf{x})$ becomes differentiable with respect to $\Sigma$.

Theorem 3**.**

Let $\mathbf{x}\in\mathbb{R}^{d}/\mathcal{B}$. There exists a neighbourhood of $\Sigma_{0}$, $\mathcal{W}_{\Sigma_{0}}$, such that $i_{\mathbf{x},\Sigma}$ is constant over this neighbourhood. Moreover, the function $\Sigma\mapsto\log Y_{\Sigma}(\mathbf{x})$ is differentiable over $\mathcal{W}_{\Sigma_{0}}$ and

[TABLE]

Proof.

It follows from Proposition 3 that

[TABLE]

and, for $\mathbf{x}\in\mathbb{R}^{d}/\mathcal{B}$, we have, for all $j\neq i_{\mathbf{x},\Sigma}$,

[TABLE]

or equivalently

[TABLE]

i.e.,

[TABLE]

Let $\zeta>0$ and define

[TABLE]

such that $\mathcal{I}_{1}\cup\mathcal{I}_{2}\cup\mathcal{I}_{3}=\{j\geq 1,j\neq i_{\mathbf{x},\Sigma}\}$. Moreover it is easy to see, using the definition of $i_{\mathbf{x},\Sigma}$ and the form of the intensity function of the point process $(U_{i},\mathbf{C}_{i})_{i\geq 1}$, that $|\mathcal{I}_{1}|$ is finite, $|\mathcal{I}_{2}|$ is finite, but $|\mathcal{I}_{3}|$ is infinite.

Let $\Theta$ be a positive-definite matrix of size $d\times d$. We have, for any $j\geq 1$,

[TABLE]

In addition,

[TABLE]

Since $\left\|\cdot\right\|$ and $\left\|\mathbf{\cdot}\right\|_{\Sigma^{-1}}$ are equivalent norms, there exists a positive constant $D$ such that

[TABLE]

It follows that, for any $j\geq 1$, there exists a function $\mathbf{x}\mapsto a_{j}\left(\mathbf{x},\Sigma,\Theta\right)$ such that

[TABLE]

and

[TABLE]

Using (31), we obtain

[TABLE]

Using (32), we see that there exists $\kappa>0$ such that, for $\left\|\Sigma^{-1}-\Theta^{-1}\right\|<\kappa$, we have, for all $j\geq 1$ such that $j\neq i_{\mathbf{x},\Sigma}$, that

[TABLE]

and, for all $j\in\mathcal{I}_{3}$,

[TABLE]

Hence, using (33), we obtain, for all $j\in\mathcal{I}_{3}$,

[TABLE]

Now, using (30), the continuity of $\Sigma^{-1}\mapsto\left\|\mathbf{\cdot}\right\|_{\Sigma^{-1}}$ and the fact that $|\mathcal{I}_{1}|$ and $|\mathcal{I}_{2}|$ are finite, there exists $\kappa^{\prime}>0$ such that, for $\left\|\Sigma^{-1}-\Theta^{-1}\right\|<\kappa^{\prime}$, we have, for all $j\in$ $\mathcal{I}_{1}\cup\mathcal{I}_{2}$, that

[TABLE]

Combining (34) and (35), we obtain that, for all positive definite matrix $\Theta$ satisfying $\left\|\Sigma^{-1}-\Theta^{-1}\right\|<\min\{\kappa,\kappa^{\prime}\}$, we have, for all $j\geq 1$ such that $j\neq i_{\mathbf{x},\Sigma}$,

[TABLE]

Accordingly, $i_{\mathbf{x},\Theta}=i_{\mathbf{x},\Sigma}$. Hence we can choose the neighbourhood of $\Sigma_{0}$

[TABLE]

to define the derivative of $\Sigma\mapsto\log Y_{\Sigma}(\mathbf{x})$ at $\Sigma_{0}$.

We now compute the corresponding derivative. Using Proposition 3, we have

[TABLE]

and hence

[TABLE]

Formula (11.7) in Dwyer, (1967) gives, for any symmetric matrix $\Sigma$, that

[TABLE]

Moreover, since $i_{\mathbf{x},\Sigma}$ is constant over $\mathcal{W}_{\Sigma_{0}}$, Equation (11.8) in Dwyer, (1967) provides, for any symmetric matrix $\Sigma$,

[TABLE]

Combining (36) and (37), we finally obtain

[TABLE]

∎

We now prove that the derivative of $\Sigma\mapsto\log Y_{\Sigma}(\mathbf{x})$ can be uniformly bounded by an integrable random variable over a neighbourhood of $\Sigma_{0}$.

Theorem 4**.**

There exists a non-random neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, such that, for any $q>1$, there exists a random variable $C_{\Sigma_{0}}(\mathbf{x},q)$ satisfying a.s.

[TABLE]

and $\mathbb{E}\left[C_{\Sigma_{0}}(\mathbf{x},q)\right]<\infty$.

Proof.

First recall that a.s.

[TABLE]

which gives a.s.

[TABLE]

Consequently, using the well-known fact that, for all $a,b\in\mathbb{R}$ and $q\geq 1$, $|a-b|^{q}\leq 2^{q-1}(|a|^{q}+|b|^{q})$, we obtain

[TABLE]

Let $A=\Sigma^{-1}(\mathbf{x}-\mathbf{C}_{i_{\mathbf{x},\Sigma}})$, such that $\Sigma^{-1}(\mathbf{x}-\mathbf{C}_{i_{\mathbf{x},\Sigma}})(\mathbf{x}-\mathbf{C}_{i_{\mathbf{x},\Sigma}})^{{}^{\prime}}\Sigma^{-1}=AA^{{}^{\prime}}$. Note that $AA^{{}^{\prime}}$ is a non-negative symmetric matrix of rank $1$ and thus it only has one positive eigenvalue given by

[TABLE]

(since we have that $(AA^{{}^{\prime}})A=A(A^{{}^{\prime}}A)=(A^{{}^{\prime}}A)A$). It follows that

[TABLE]

For a positive definite symmetric matrix $\Theta$, we denote by $\lambda_{\max}\left(\Theta\right)$ and $\lambda_{\min}\left(\Theta\right)$ respectively the maximum and the minimum of its positive eigenvalues. We have that

[TABLE]

and

[TABLE]

which yield

[TABLE]

Combining (38), (39) and (40), we have, for any $q>1$, that a.s.

[TABLE]

where

[TABLE]

In order to control $C_{\Sigma_{0}}(\mathbf{x},q)$, it is sufficient to control $\|\mathbf{x}-\mathbf{C}_{i_{\mathbf{x},\Sigma}}\|_{\Sigma^{-1}}^{2q}$. As the center of the “storm” realizing the maximum at point $\mathbf{x}$, $\mathbf{C}_{i_{\mathbf{x},\Sigma}}$ is characterized by

[TABLE]

and therefore

[TABLE]

Additionally, we have, for all $i\geq 1$,

[TABLE]

which yields, by equivalence of the norms $\|.\|$ and $\|.\|_{\Sigma^{-1}}$, that

[TABLE]

for some positive constant $D_{0}$. Hence let us now choose $\mathcal{V}_{\Sigma_{0}}$ such that

[TABLE]

Now, for two real-valued random variables $U$ and $V$, $U+V\geq\lambda$ implies that $U\geq\lambda/2$ or $V\geq\lambda/2$, giving $\mathbb{P}(U+V\geq\lambda)\leq\mathbb{P}\left(U\geq\lambda/2\right)+\mathbb{P}\left(V\geq\lambda/2\right).$ Thus, using (43) and (44), we obtain

[TABLE]

We first deal with the first term of the right-hand side of (46). Since the $(U_{i},\mathbf{C}_{i})_{i\geq 1}$ are the points of a Poisson process on $(0,\infty)\times\mathbb{R}^{d}$ with intensity function $u^{-2}\mathrm{d}u\times\mathrm{d}\mathbf{c}$, we have

[TABLE]

where

[TABLE]

Making the change of variable $v=\lambda/2+2\log(u)\mbox{, yielding }u=\exp\left(v/2-\lambda/4\right)$ and $\mathrm{d}u=\exp\left(v/2-\lambda/4\right)\mathrm{d}v/2$, we obtain

[TABLE]

The latter integral is finite and it follows from (47) that

[TABLE]

Let us now deal with the second term of the right-hand side of (46). Observe that

[TABLE]

The mapping $u\rightarrow u^{-1}$ applied to the points of the Poisson point process yields a new Poisson point process with intensity function $\mathrm{d}u$. Furthermore, such a Poisson point process on $(0,\infty)$ is homogeneous and can be represented as the sum of independent standard exponential random variables. Thus, we can write $\min_{i\geq 1}U_{i}^{-1}{\overset{d}{=}}\mathrm{Exp}(1)$. Moreover, $\inf_{\Sigma\in\mathcal{V}_{\Sigma_{0}}}U_{i_{\mathbf{x},\Sigma}}^{-1}\geq\min_{i\geq 1}U_{i}^{-1}$, implying

[TABLE]

which gives for large $\lambda$

[TABLE]

Now, for any $q>1$, we have

[TABLE]

We carry out the change of variable $\lambda=u^{1/q}$, which gives $u=\lambda^{q}$ and hence $\mathrm{d}u=q\lambda^{q-1}\mathrm{d}\lambda$. Accordingly, we obtain

[TABLE]

and therefore, using (46), (49) and (50),

[TABLE]

Finally, let us recall that

[TABLE]

Consequently, using (45), we finally deduce that any neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, satisfying

[TABLE]

leads to

[TABLE]

∎

We finally provide the proof of Theorem 2.

Proof.

Let us choose $\mathcal{V}_{\Sigma_{0}}$ as in the proof of Theorem 4. It follows from (17) that

[TABLE]

which gives

[TABLE]

Hence we choose

[TABLE]

We have

[TABLE]

Let $q>1$ such that $p^{-1}+q^{-1}=1$. By Hölder inequality, we have

[TABLE]

By Theorem 4 and the assumptions, the result follows. ∎

A.3 Proof of Proposition 2

A.3.1 LRM for the Brown–Resnick random field

Let us first recall that $M=2$ and that $H$ is given by (25) with $\beta_{i}\in\mathbb{N}_{\ast}$, $i=1,2$. Let us now note that $\left\|\mathbf{Y}_{\boldsymbol{\theta}_{0}}\right\|^{\alpha}\leq 2^{-\alpha/2}\left(Y_{\boldsymbol{\theta}_{0},1}\vee Y_{\boldsymbol{\theta}_{0},2}\right)^{\alpha}\leq 2^{-\alpha/2}\left(Y_{\boldsymbol{\theta}_{0},1}^{\alpha}+Y_{\boldsymbol{\theta}_{0},2}^{\alpha}\right)$. Therefore there exists a positive constant $C$ such that

[TABLE]

where, for $i=1,2$, $\Delta_{i}$ is a set of constants whose the maximal value is $\xi_{i}\beta_{i}+\alpha$. Moreover,

[TABLE]

Since $Y_{\boldsymbol{\theta}_{0},i}$ has a standard Fréchet distribution, the expectation $\mathbb{E}\left[\left(Y_{\boldsymbol{\theta}_{0},i}^{2\delta_{i}}\exp\left(-2B_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}Y_{\boldsymbol{\theta}_{0},i}^{-1}\right)\right)\right]$ is finite, if $2\left(\xi_{i}\beta_{i}+\alpha\right)<1$ and $-2B_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}<1$. But it is possible to choose $\alpha$ and $\mathcal{V}_{\boldsymbol{\theta}_{0}}$ to satisfy such constraints since $\xi_{i}\beta_{i}<1/2$ and $B_{\mathcal{V}_{\boldsymbol{\theta}_{0}}}$ is defined by (9).

Finally note that the derivatives $\partial\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{1},\mathbf{x}_{2})/\partial\psi$ and $\partial\lambda_{\boldsymbol{\theta}}(\mathbf{x}_{1},\mathbf{x}_{2})/\partial\kappa$ are easily obtained and therefore it is straightforward to conclude that the condition in (11) also holds with such a choice of neighbourhood $\mathcal{V}_{\boldsymbol{\theta}_{0}}$.

A.3.2 IPA for the Smith random field

From (25), there exist positive constants $C_{i}$, $i=1,2$, such that

[TABLE]

Therefore, we have for $p>1$,

[TABLE]

Next result will allow us to prove that Condition (18) holds.

Proposition 5**.**

There exists a non-random neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, such that, for any $\beta<1$ and $\mathbf{x\in}\mathbb{R}^{d}$,

[TABLE]

Proof.

We have

[TABLE]

and thus

[TABLE]

It is well-know that $\bigvee_{i=1}^{\infty}U_{i}$ has a standard Fréchet distribution and therefore $\mathbb{E[}(\bigvee_{i=1}^{\infty}U_{i})^{\beta}]<\infty$ since $\beta<1$. We deduce that it suffices to choose $\mathcal{V}_{\Sigma_{0}}$ such that $\sup_{\Sigma\in\mathcal{V}_{\Sigma_{0}}}\det\left(\Sigma\right)^{-1/2}<\infty$ and $\sup_{\Sigma\in\mathcal{V}_{\Sigma_{0}}}\det\left(\Sigma\right)^{1/2}<\infty$. This is possible since any symmetric invertible matrix $A$ admits a neighbourhood $\mathcal{V}$ of matrices such that, for all $\tilde{A}\in\mathcal{V}$, $\det(A)-\varepsilon<\det(\tilde{A})<\det(A)+\varepsilon$, where $0<\varepsilon<\det(A)$. ∎

By Proposition 5, $\left(\ref{Eq_Bound_sup_derivative}\right)$ and the Hölder inequality, there exists a non-random neighbourhood of $\Sigma_{0}$, $\mathcal{V}_{\Sigma_{0}}$, such that for some $p>1$ satisfying $2p\xi_{1}\beta_{1}<1$ and $2p\xi_{2}\beta_{2}<1$, we have

[TABLE]

Appendix B Analytical formulas

B.1 Expressions for the derivatives of $R(\boldsymbol{\theta})$

when $H$ is given by (25)

For $\tilde{\beta_{1}},\tilde{\beta_{2}}<1/2$, we introduce the function $g_{\tilde{\beta_{1}},\tilde{\beta_{2}}}$ defined by

[TABLE]

where, for $t,h>0$,

[TABLE]

with $\Phi$ and $\varphi$ denoting the standard univariate Gaussian distribution and density functions, respectively.

Let $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{2}$ be our sites of interest and $Y_{\boldsymbol{\theta}}$ be a Brown–Resnick random field with dependence parameters $\psi$ and $\kappa$ gathered in $\boldsymbol{\theta}$, i.e., $\boldsymbol{\theta}=\left(\psi,\kappa\right)^{\prime}$. Let $\boldsymbol{\theta}_{0}$ be the parameter at which we wish to compute the derivative of $R(\boldsymbol{\theta})$. As before, let $\eta_{i},\tau_{i},\xi_{i}$ be the location, scale and shape parameters of $X_{\boldsymbol{\theta}}(\mathbf{x}_{i})$, $i=1,2$, and $\beta_{i}=\beta(\mathbf{x}_{i})\in\mathbb{N}_{\ast}$ such that $\beta_{i}\xi_{i}<1/2$. Theorem 2 in Koch, (2019) yields

[TABLE]

where

[TABLE]

This yields

[TABLE]

We obtain from (52) that

[TABLE]

where $h_{0}=\sqrt{2\gamma_{\boldsymbol{\theta}_{0}}(\mathbf{x}_{2}-\mathbf{x}_{1})}$. The term

[TABLE]

has a closed (although complicated) expression (available upon request). The corresponding integral can be computed using numerical methods such as adaptive quadrature. Finally, we obtain the true values of $\partial R(\boldsymbol{\theta})/\partial\boldsymbol{\theta}|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{0}}$ by plugging the appropiate values of (55) in (54).

Now, if $Y_{\Sigma}$ is the Smith random field with covariance matrix $\Sigma$ and $\Sigma_{0}$ is a symmetric positive-definite matrix at which we want to compute the derivative of $R\left(\boldsymbol{\theta}\right)$, exactly the same formulas and procedure can be applied upon replacement of $\boldsymbol{\theta}$ with $\Sigma$, $\boldsymbol{\theta}_{0}$ with $\Sigma_{0}$, and $\gamma_{\boldsymbol{\theta}}(\mathbf{x}_{2}-\mathbf{x}_{1})$ with $(\mathbf{x}_{2}-\mathbf{x}_{1})^{\prime}\Sigma^{-1}(\mathbf{x}_{2}-\mathbf{x}_{1})/2$.

B.2 Proportionality of the components of the score function for the

Brown–Resnick random field

The bivariate density of the (simple) Brown–Resnick random field (at $\mathbf{x}_{1}$ and $\mathbf{x}_{2}\in\mathbb{R}^{2}$) satisfies, for $y_{1},y_{2}>0$,

[TABLE]

where

[TABLE]

This easily yields

[TABLE]

and it follows from (6) that the LRM estimates of $\partial R(\boldsymbol{\theta})/\partial\psi$ and $\partial R(\boldsymbol{\theta})/\partial\kappa$ are proportional by the factor given in the right-hand side of (57).