ExtremalDep: Modelling extremal dependence in high-dimensional extremes

Let ${\bm{X}}_{1},{\bm{X}}_{2},\ldots$ be independent and identically distributed (i.i.d.) random vectors, where ${\bm{X}}_{i}=(X_{i,1},\ldots,X_{i,d})\in{\mathbb{R}}^{d}$ with $d=2,3,\ldots$, whose joint distribution $F$ is in the domain of attraction of a multivariate Generalised Extreme Value (GEV) distribution $G$ (de2006, chap. 6), shortly denoted as $F\in\mathcal{D}(G)$. This implies that for $n=1,2,\ldots$, there are norming sequences ${\bm{a}}_{n}>{\bf 0}=(0,\ldots,0)$ and ${\bm{b}}_{n}\in{\mathbb{R}}^{d}$ such that

$$\lim_{n\to\infty}F^{n}({\bm{a}}_{n}{\bm{x}}+{\bm{b}}_{n})=G_{\bm{\gamma}}({\bm{x}}\mid D),\quad{\bm{x}}\in{\mathbb{R}}^{d},$$

(2.1)

where ${\bm{\gamma}}=(\gamma_{1},\ldots,\gamma_{d})^{\top}\in{\mathbb{R}}^{d}$. The limiting distribution can be written as

$$G_{{\bm{\gamma}}}({\bm{x}}\mid D)={C_{\scalebox{0.65}{EV}}}\left(G_{\gamma_{1}}(x_{1}),\ldots,G_{\gamma_{d}}(x_{d})\mid D\right),\qquad{\bm{x}}\in{\mathbb{R}}^{d},$$

(2.2)

where ${C_{\scalebox{0.65}{EV}}}$ is an extreme-value copula (beirlant2004, chap. 8) and for all $j\in\{1,\ldots,d\}$ the marginal distributions are of the form

$$G_{\gamma_{j}}(x_{j})=\exp\left(-\left(1+\gamma_{j}x_{j}\right)^{-\frac{1}{\gamma_{j}}}\right),\quad 1+\gamma_{j}x_{j}>0,$$

(2.3)

i.e., univariate GEV distributions with extreme value index $\gamma_{j}\in{\mathbb{R}}$ describing the tail heaviness of the distribution (e.g., de2006, chap. 1). In particular, the extreme-value copula takes the form

$${C_{\scalebox{0.65}{EV}}}({\bm{u}}\mid D)=\exp\left(-L\left((-\ln u_{1}),\ldots,(-\ln u_{d})\right)\right),\quad{\bm{u}}\in(0,1]^{d},$$

(2.4)

where $L:[0,\infty)^{d}\mapsto[0,\infty)$ is a homogeneous function of order $1$ named stable-tail dependence function (beirlant2004, chap. 8). This function is fully characterised by its restriction on ${\mathcal{R}}:=\{{\bm{t}}\in[0,1]^{d-1}:\|{\bm{t}}\|_{1}\leqslant 1\}$ which is known as the Pickands dependence function (e.g., falk2010, Ch. 4) and defined as

$$A({\bm{t}})=L(1-t_{1}-\dots-t_{d-1},{\bm{t}})=d\int_{\mathcal{S}}\max\{(1-t_{1}-\cdots-t_{d-1})w_{1},\ldots,t_{d-1}w_{d}\}\mathrm{d}H({\bm{w}}),$$

(2.5)

where $H$ is a probability measure on the $d$-dimensional unit simplex $\mathcal{S}:=\{{\bm{w}}\geqslant{\bf 0}:\|{\bm{w}}\|_{1}=1\}$, named angular measure (e.g., falk2010, Ch. 4). A valid angular measure is a probability measure that satisfies the mean constraint (C1), while a valid Pickands dependence function is a function that must satisfies the convexity and boundary conditions (C2)-(C3). Specifically,

1. (C1)

   $\int_{\mathcal{S}}w_{j}H(\mathrm{d}{\bm{w}})=1/d,\;\forall\;j\in\{1,\ldots,d\}$.

2. (C2)

   $A(a{\bm{t}}_{1}+(1-a){\bm{t}}_{2})\leqslant aA({\bm{t}}_{1})+(1-a)A({\bm{t}}_{2}),\,a\in[0,1],\,\forall\,{\bm{t}}_{1},{\bm{t}}_{2}\in{\mathcal{R}}$,

3. (C3)

   $1/d\leqslant\max\left(t_{1},\ldots,t_{d-1},1-t_{1}-\cdots-t_{d-1}\right)\leqslant A({\bm{t}})\leqslant 1,\,\forall\,{\bm{t}}\in{\mathcal{R}}$.

Conditions (C2)-(C3) are necessary and sufficient to define a valid Pickands dependence function only in the bivariate case, (for details see beirlant2004, chap. 8). Since $A$ and $H$ are related through the one-to-one map in (2.5), the dependence parameter $D$ of the copula ${C_{\scalebox{0.65}{EV}}}(\cdot\mid D)$ is commonly meant to be either the angular measure or the Pickands dependence function. Both are possible means to describe the so-called extremal dependence.

Given a sample $({\bm{X}}_{1},\ldots,{\bm{X}}_{N})$ with joint distribution $F\in\mathcal{D}(G)$, for each margin $j\in\{1,\ldots,d\}$, $k$ maxima are computed over blocks of size $n$, such that $N=n\cdot k$, where the $i$th maximum given by $M_{i,j}^{(n)}=\max(X_{(i-1)n+1,j},\ldots,X_{in,j})$, $i\in\{1,\ldots,k\}$. The vector of componentwise maxima is then defined as ${\bm{M}}_{n,i}=(M_{i,j}^{(n)},1\leqslant j\leqslant d)$. Setting ${\bm{\mu}}={\bm{b}}_{n}$, ${\bm{\sigma}}={\bm{a}}_{n}$, ${\bm{y}}={\bm{a}}_{n}{\bm{x}}+{\bm{b}}_{n}$, by (2.1), the joint distribution of ${\bm{M}}_{n,i}$ can be approximated, for large enough $n$, as $F^{n}({\bm{y}})\approx G_{{\bm{\theta}}}({\bm{y}}\mid D),$ where ${\bm{\theta}}=({\bm{\mu}},{\bm{\sigma}},{\bm{\gamma}})^{\top}$ and each one-dimensional marginal distribution regarding the maximum $M_{i,j}^{(n)}$ can be approximated as

$$F^{n}(y_{j})\approx G_{\theta_{j}}(y_{j})=\exp\left(-\left(1+\gamma_{j}\frac{y_{j}-\mu_{j}}{\sigma_{j}}\right)^{-\frac{1}{\gamma_{j}}}\right),\quad j\in\{1,\ldots,d\}.$$

(2.6)

The dependence among the maxima $\left(M_{i,1}^{(n)},\ldots,M_{i,d}^{(n)}\right)$ can be approximated by the extreme-value copula ${C_{\scalebox{0.65}{EV}}}(\cdot\mid D)$, where if the parameter $D$ is the angular measure $H$, it describes the dependence among maxima as follows: the more mass $H$ places around the centre of the simplex $(1/d,\ldots,1/d)$ the more dependent the maxima are. Conversely, the maxima becomes less dependent as $H$ concentrates its mass close to the vertices of the simplex. Alternatively, $D$ can be the Pickands dependence function with the lower and upper bounds in the inequality in (C3) respectively representing the cases of complete dependence and independence. A useful summary of the dependence among the maxima is given by the so-called extremal coefficient (e.g., beirlant2004, Ch. 8), defined as

$$\eta=dA(1/d,\ldots,1/d),$$

(2.7)

which represents the (fractional) number of independent block maxima. The extremal coefficient satisfies $1\leqslant\eta\leqslant d$, with the lower and upper bounds describing the case of complete dependence and independence.

Beyond these results, the domain of attraction provides also a useful basis to derive very small joint tail probabilities corresponding to high-dimensional extremes relative to the original distribution $F$ or to derive extreme quantile regions corresponding to very small joint tail probabilities. These results are summarized as follows. The limit in (2.1) holds if and only if $n(1-F({\bm{a}}_{n}{\bm{x}}+{\bm{b}}_{n}))\to-\log G({\bm{x}})$ as $n\to\infty$, for all ${\bm{x}}\in{\mathbb{R}}^{d}$, with the convention that $-\log(0)=\infty$. Together with the weak convergence of the copula of the vector of componentwise maxima, namely $\{C({\bm{u}}^{1/n})\}^{n}\to{C_{\scalebox{0.65}{EV}}}({\bm{u}})$ as $n\to\infty$, for all ${\bm{u}}\in(0,1]^{d}$, this result implies that for all ${\bm{x}}\in[0,\infty)^{d}$, as $n\to\infty$,

$$n\,{\mathbb{P}}\left(\bigcup_{j=1}^{d}F_{j}(X_{j})>1-\frac{x_{j}}{n}\right)\to L({\bm{x}})=d\int_{\mathcal{S}}\max_{j=1,\ldots,d}{(x_{j}}w_{j})\mathrm{d}H({\bm{w}}).$$

Recall that a copula $C$ defined on $[0,1]^{d}$ and associated to a joint distribution function $F$, is a suitable function that allows the definition $F({\bm{x}})=C(F_{1}(x_{1}),\ldots,F_{d}(x_{d}))$ (see e.g., joe2014dependence, for details). For simplicity we assume that the margins of $F$ are continuous. The result in the above display implies the following important approximation. Let $p_{j}=p_{n,j}=x_{j}/n$ for $j=1,\ldots,d$ and $Q_{j}(p_{j})$ be the ($1-p$)-quantile of $F_{j}$, then by homogeneity of the stable tail dependence function, the probability that at least one among the $d$ variables $X_{j}$ exceeds a high quantile of its own distribution, can be approximated for large $n$ as

$${\mathbb{P}}(X_{1}>Q_{1}(p_{1})\,\text{or}\cdots\text{or}\,X_{d}>Q_{d}(p_{d}))\approx L(p_{1},\ldots,p_{d}).$$

(2.8)

By the inclusion and exclusion principle, the weak convergence result of the extreme value copula also allows to obtain

$$n\,{\mathbb{P}}\left(\bigcap_{j=1}^{d}F_{j}(X_{j})>1-\frac{x_{j}}{n}\right)\to R({\bm{x}})=d\int_{\mathcal{S}}\min_{j=1,\ldots,d}{(x_{j}}w_{j})\mathrm{d}H({\bm{w}}),\quad n\to\infty,$$

where $R$ is known as the tail copula function, which is also a homogeneous function of order $1$, and therefore the probability that all the $d$ variables $(X_{1},\ldots,X_{d})$ exceed simultaneously high quantiles of their own distributions can be approximated, for a sufficiently large $n$, as

$${\mathbb{P}}(X_{1}>Q_{1}(p_{1})\,\text{and}\cdots\text{and}\,X_{d}>Q_{d}(p_{d}))\approx R(p_{1},\ldots,p_{d}).$$

(2.9)

Finally, a related question is, how to determine an extreme region given a small joint probability $p$ of falling in it? The approach discussed by cai2011, einmahl2013 and he2017 is summarized as follows with a focus on the bivariate case explored by einmahl2013. Let $f$ be the density on ${\mathbb{R}}_{+}^{2}$ of a bivariate distribution $F$ satisfying $F\in\mathcal{D}(G)$, with$\gamma_{1},\gamma_{2}>0$. Define quantile regions by the level sets of $f$ such that

$$\mathcal{Q}=\{{\bm{x}}\in{\mathbb{R}}^{2}_{+}:f({\bm{x}})\leqslant\alpha\},$$

(2.10)

with $\alpha>0$. Then, for a small probability $p>0$ we derive the region $\mathcal{Q}$ such that ${\mathbb{P}}(\mathcal{Q})=p$ and $f$ everywhere on $\mathcal{Q}$ is smaller than $f$ everywhere on $Q^{\complement}$ and therefore the latter is the smallest region satisfying ${\mathbb{P}}(\mathcal{Q}^{\complement})=1-p$. Accordingly, an extreme quantile region is defined by a level set $\mathcal{Q}_{N}$ (with $\alpha=\alpha_{N}$ in (2.10)), which is such that ${\mathbb{P}}(\mathcal{Q}_{N})=p$ where $p=p_{N}\to 0$ and the expected number of points falling in it is $Np\to c>0$ as $N\to\infty$. Note that when $c<1$, the region is so extreme that almost no observations are expected to fall in it.

The domain of attraction condition implies the existence of a measure $\xi$, named exponent measure, such that

$$N{\mathbb{P}}(\{Q_{1}(1/N)x_{1}^{\gamma_{1}},Q_{2}(1/N)x_{2}^{\gamma_{2}}:{\bm{x}}\in B\})\to\xi(B),\quad N\to\infty$$

and satisfying $\xi(\partial B)=0$, for all Borel set $B\subset[0,\infty]^{2}$ that are bounded away from the origin. Assume that $f$ is bounded away from zero on $(0,M]^{2}$ and, outside this region, $f$ is decreasing in each coordinate, for some $M>0$. Assume also that the exponent measure admits a density function $g$ such that

$$\xi(\{{\bm{v}}:v_{1}>x_{1}\,\text{ or }\,v_{2}>x_{2}\})=\int\int_{\{v_{1}>x_{1}\text{ or }v_{2}>x_{2}\}}g({\bm{v}})\mathrm{d}{\bm{v}},$$

and

$$NQ_{1}(1/N)Q_{2}(1/N)f(Q_{1}(1/N)x^{\gamma_{1}},Q_{2}(1/N)x^{\gamma_{2}})\to\frac{x_{1}^{1-\gamma_{1}}x_{2}^{1-\gamma_{2}}}{\gamma_{1}\gamma_{2}}g({\bm{x}})=:q({\bm{x}}),\quad N\to\infty.$$

In particular, we have that $g$ is related to the density $h$ of the angular measure $H$ by $h(w)=2^{-1}g(w,1-w)$, where $w=x_{1}/r$ and $r=x_{1}+x_{2}$. The basic set $S$ is defined as $S=\{{\bm{x}}:q({\bm{x}})\leqslant 1\}=\{{\bm{x}}:r\geqslant q_{\star}^{-1}(w),w\in[0,1]\}$, where

$$q_{\star}(w)=\left(2\gamma_{1}^{-1}\gamma_{2}^{-1}w^{1-\gamma_{1}}(1-w)^{1-\gamma_{2}}h(w)\right)^{-\frac{1}{1+\gamma_{1}+\gamma_{2}}},$$

and, according to the exponent measure, its size is equal to

$$\xi(S)=2\int_{[0,1]}q_{\star}(w)h(w)\mathrm{d}w.$$

The set $S$ is then suitably transformed in order to obtain for $\mathcal{Q}_{N}$ the approximation

$$\widetilde{\mathcal{Q}}_{N}=\left\{\left(\mu_{1}+\sigma_{1}\frac{\left(\frac{k\,\xi(S)x_{1}}{Np}\right)^{\gamma_{1}}-1}{\gamma_{1}},\mu_{2}+\sigma_{2}\frac{\left(\frac{k\,\xi(S)x_{2}}{Np}\right)^{\gamma_{2}}-1}{\gamma_{2}}\right):{\bm{x}}\in S\right\},$$

(2.11)

where $k=k_{N}$ with $k\to\infty$ as $N\to\infty$ and $k=o(N)$. Furthermore, $\mu_{j}=Q_{j}(k/N)$ and $\sigma_{j}=a_{j}(N/k)$ with $j=1,2$, where $a_{j}(\cdot)$ is a suitable positive scale function (de2006, Ch. 1-2). $\widetilde{\mathcal{Q}}_{N}$ is a good approximation of $\mathcal{Q}_{N}$ in the sense that ${\mathbb{P}}(\widetilde{\mathcal{Q}}_{N})\approx p$ and ${\mathbb{P}}(\mathcal{Q}_{N}\Delta\widetilde{\mathcal{Q}}_{N})\to 0$ as $N\to\infty$ (see einmahl2013), where $A\Delta B=A\setminus B\cup B\setminus A$.

From a statistical perspective, the multivariate GEV model defines a complex semiparametric family of distributions of the form $G_{{\bm{\theta}}}({\bm{y}}\mid D)={C_{\scalebox{0.65}{EV}}}(G_{\theta_{1}}(y_{1}),\ldots,G_{\theta_{d}}(y_{d})\mid D)$, ${\bm{y}}\in{\mathbb{R}}^{d}$, where ${\bm{\theta}}\in\Theta\subset({\mathbb{R}},{\mathbb{R}}_{+},{\mathbb{R}})^{d}$ is a finite-dimensional vector of marginal parameters and $D$ is an infinite-dimensional dependence parameter of the copula ${C_{\scalebox{0.65}{EV}}}$, named extremal dependence, living in a suitable space, depending whether $D$ is the angular measure or the Pickands dependence function. The estimation of such infinite-dimensional parameter, accounting for the corresponding constraints (see conditions (C1)-(C3) of Section 2.1), is far from being a simple task. To reduce complexity, classes of parametric models were initially proposed to model $D$, see Beranger2015 for a review. In order to provide valid extremal dependence structures, the proposed models satisfy the required constraints. Section 2.3 describes how to fit some of the most popular parametric models for $D$ using ExtremalDep.

Over time, more sophisticated semi- and non-parametric approaches for modelling and estimating the extremal dependence have been proposed, under both parametrisations. Some examples include: projection estimators (e.g., fils2008), polynomials (e.g., Marcon2017a) and splines (e.g., cormier2014), just to name a few. Section 2.4 details how to use ExtremalDep for non-parametric estimation of the extremal dependence, simulation of bivariate extreme values with a flexible semi-parametric dependence structure, estimation of small probabilities to belong to certain extreme sets and estimation extreme quantile regions corresponding to small joint tail probabilities.

A non-exhaustive list of popular parametric models for the extremal dependence includes the Asymmetric-Logistic, Pairwise-Beta, Tilted-Dirichlet, Hüsler-Reiss, Extremal-$t$ and Extremal-Skew-$t$ models (see e.g., Beranger2015; Beranger2017).

The Poisson-Point-Process (PPP) is a simple estimation method for statistical models used to represent the extremal dependence, when the latter is parametrized according to the angular measure and provided that it allows for a density $h(w|{\bm{\varphi}})$ on $\mathcal{S}$, where ${\bm{\varphi}}$ is a suitable vector of parameters (see e.g., coles1991; engelke2015). In particular, let ${\bm{x}}_{1},\ldots,{\bm{x}}_{N}$ be i.i.d. observations from $F\in\mathcal{D}(G)$ and ${\bm{y}}_{1},\ldots,{\bm{y}}_{N}$ be their equivalent with unit Fréchet margins. The pseudopolar transformation $T$ is defined as $T({\bm{y}}_{i})=(r,{\bm{w}}_{i})$ where $r={\bm{y}}_{1}+\cdots+{\bm{y}}_{N}$ is the radial part and ${\bm{w}}_{i}={\bm{y}}_{i}/r$ is the angular part of ${\bm{y}}_{i}$, $i=1,\ldots,N$. The domain of attraction condition implies that if $r_{0}$ is a large threshold and $W_{r_{0}}=\{(r,{\bm{w}}):r>r_{0}\}$ is a set of points with a radial component larger than $r_{0}$, then the number of points falling in $W_{r_{0}}$ is described approximatively by a Poisson random variable with rate $1/\psi(W_{r_{0}})$, where $\psi=\xi(T^{-1}(\cdot))$. Thus, conditionally on observing $k$ points $\{(r_{(i)},{\bm{w}}_{(i)}),i=1,\ldots,k\}$ with radial part greater than $r_{0}$, such points are approximately independent and generated from the density function $(r^{-2}\mathrm{d}r\times H(\mathrm{d}{\bm{w}}))/\psi(W_{r_{0}})$. The likelihood function of the exceedances can then be approximated as

	$\displaystyle L({\bm{\varphi}}\mid(r_{(i)},{\bm{w}}_{(i)}),i=1,\ldots,k)$	$\displaystyle\approx\frac{e^{-\psi(W_{r_{0}})}(\psi(W_{r_{0}}))^{k}}{k!}\prod_{i=1}^{k}\frac{r_{(i)}^{-2}h({\bm{w}}_{(i)}\mid{\bm{\varphi}})}{\psi(W_{r_{0}})}$
		$\displaystyle\propto h({\bm{w}}_{(i)}|{\bm{\varphi}}).$

Estimates of ${\bm{\varphi}}$ can be obtained by maximizing the log-likelihood $\ell({\bm{\varphi}})=\sum_{i=1}^{k}\log h({\bm{w}}_{(i)}\mid{\bm{\varphi}})$. Alternatively, one can appeal to a Bayesian approach (sabourin2013), based on the approximate posterior density of the angular density parameters given by

$$\pi({\bm{\varphi}}\mid{\bm{w}}_{(1)},\ldots,{\bm{w}}_{(k)})=\frac{\prod_{i=1}^{k}h({\bm{w}}_{(i)}\mid{\bm{\varphi}})\phi({\bm{\varphi}})}{\int_{\Psi}\prod_{i=1}^{k}h({\bm{w}}_{(i)}\mid{\bm{\varphi}})\phi({\bm{\varphi}})\mathrm{d}{\bm{\varphi}}},$$

where $\phi({\bm{\varphi}})$ is a suitable prior distribution on ${\bm{\varphi}}$ and $\Psi$ is a suitable parameter space.

The routine fExtDep with argument method="PPP" allows to estimate ${\bm{\varphi}}$ by likelihood maximisation. The model class for $h$ is specified with the argument model with the list of available options reported in Table 2.3. Some classes of models are characterised by an angular density $h({\bm{w}}\mid{\bm{\varphi}})$ that places a continuous mass only in the interior of $\mathcal{S}$ (PB, TD, HR), while others also allow to place densities and atoms on subspaces of $\mathcal{S}$. Beranger2017 describe an extended version of the ML method that allows for the estimation of these more complex models up to the three-dimensional case (AL, ET, EST) by exploiting the following approach. In the bivariate case, an observation ${\bm{w}}=(w,1-w)$ falls at an endpoint of $\mathcal{S}$ if $w<c$ or $w>1-c$ for some small $c\in(0,1/2)$ and on the interior of $\mathcal{S}$ otherwise. A trivariate observation ${\bm{w}}=(w_{1},w_{2},1-w_{1}-w_{2})$ falls: at the $j$th corner of $\mathcal{S}$ if ${\bm{w}}\in\mathcal{C}_{j}=\left\{w_{j}>1-c\right\}$, on the edge between the $j$th and $l$th components if ${\bm{w}}\in\mathcal{E}_{j,l}=\left\{w_{j},w_{l}<1-c,w_{j}+w_{l}>1-c,w_{j}>1-2w_{l},w_{l}>1-2w_{j}\right\}$, and in the interior of $\mathcal{S}$ otherwise (i.e., if $w_{j}>c,$ for all $j$) (see Beranger2017, for details). This approach is implemented when specifying a positive value for the argument c of fExtDep.

We now illustrate the use of fExtDep on the analysis of high levels of air pollution recorded in Leeds, UK, over the winter period November 1st to February 27/28th, between 1994 and 1998 (see heffernan2004). The object pollution contains the daily maximum of the five air pollutants: particulate matter (PM10), nitrogen oxide (NO), nitrogen dioxide (NO2), ozone (03), and sulfur dioxide (SO2), in unit Fréchet scale (for a description on the used transformation see Beranger2015). Furthermore, the datasets PNS, NSN and PNN contains the angular components corresponding to the 100 largest radial components of the triplets of pollutants (PM10, NO, SO2), (NO2, SO2, NO) and (PM10, NO, NO2), respectively. We run the following simple commands

R> data(pollution)R> f.et05 <- fExtDep(method="PPP", data=PNS, model = "ET",+                    par.start = c(0.5, 0.5, 0.5, 3), trace=2, c=0.05)

initial  value -118.269640iter  10 value -146.907539iter  10 value -146.907539final  value -146.907539converged

R> f.et <- fExtDep(method="PPP", data=PNS, model = "ET",+                  par.start = c(0.5, 0.5, 0.5, 3), trace=2 )

initial  value -161.078869iter  10 value -225.654851iter  20 value -234.524086iter  30 value -235.285461iter  40 value -235.360084final  value -235.382293converged The first call to the fExtDep routine fits the Extremal-$t$ angular density (model="ET"), to the pollution data (data=PSN), allowing for point masses at the corners by specifying the argument c, while the second call considers a density defined only in the interior of $\mathcal{S}$. The argument par.start gives starting values for the parameters and trace=2 allows to monitor optimization progress. The optimization method is set by the argument optim.method ("BFGS" by default) on which basis the routine optim from the stats library is used. fExtDep returns a list with the estimates of the model parameters (par), their standard errors (SE), the log-likelihood maximum (LL) and the Takeuchi Information Criterion (TIC). fExtDep is first ran for the extremal-$t$ model using different values of $c$ ($0,0.05$ and $0.1$) to take into account point masses at the corners of $\mathcal{S}$. The TICs indicate that only assuming mass in the interior of $\mathcal{S}$ provides the best results (lowest TIC). A performance comparison between the PB, AL, TD, HR, ET and EST models defined on the interior only, is performed with TIC scores reported in Table 1 and as a result the Extremal skew-$t$ model fits the data best (lowest TIC).

Bayesian inference for the parameter ${\bm{\varphi}}$ is implemented through a model averaging approach (sabourin2013) and available in the routine fExtDep when method="BayesianPPP". Table 2.3 reports the implemented prior densities for each angular density model. In particular, $\phi(\cdot,a,b)$ stands for a one-dimensional Gaussian density with mean $a$ and variance $b$, $t(\cdot)$ is a transformation applied to the parameter (to map it to the real line). For each component of ${\bm{\varphi}}$, the proposal density is also Gaussian but of the form $\phi(t(\varphi_{j}^{\star});t(\varphi_{j}^{(i)}),\tt{MCpar})$, where $\varphi_{j}^{\star}$ and $\varphi_{j}^{(i)}$ are the proposed value and the value at the $i$th iteration of the algorithm of the $j$ element of ${\bm{\varphi}}$, respectively, and MCpar is the variance term that needs to be specified by the user. The next display illustrates how to perform Bayesian inference for the Hüssler-Reiss model (model="HR") considering Nsim=5e+4 iterations, a burn-in period of Nbin=3e+4 iterations, hyper-parameters Hpar=Hpar.hr, MCpar=0.35 and using the argument seed for reproducibility . Refer to the help page for a description of the list returned as output from fExtDep.

R> Hpar.hr <- list(mean.lambda=0, sd.lambda=3)R> PNS.hr <- fExtDep(method="BayesianPPP", data=PNS, model="HR", Nsim=5e+4,+                    Nbin=3e+4, Hpar=Hpar.pb, MCpar=0.35, seed=14342)

|=====================================================================| 100%

R> labs <- c(expression(PM[10]),expression(NO),expression(SO[2]))R> plot.ExtDep(object="angular", model="HR", par=PNS.hr$emp.mean,+              data=PNS,  cex.lab=2, labels=labs)