Extreme Value Statistics in Semi-Supervised Models

Abstract

We consider extreme value analysis in a semi-supervised setting, where we observe, next to the n data on the target variable, n + m data on one or more covariates. This is called the semi-supervised model with n labeled and m unlabeled data. By exploiting the tail dependence between the target variable and the covariates, we derive estimators for the extreme value index and extreme quantiles of the target variable in this setting and establish their asymptotic behavior. Our estimators substantially improve the univariate estimators, based on only the n target variable data, in terms of asymptotic variances whereas the asymptotic biases remain unchanged. A simulation study confirms the substantially improved behavior of both estimators. Finally the estimation method is applied to rainfall data in France. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

KEYWORDS:

• Asymptotic normality
• Extreme value index
• Semi-supervised inference
• Tail dependence
• Variance reduction

1 Introduction

The semi-supervised model, initially introduced in machine learning, deals with unbalanced datasets, when the labeled data are harder (more expensive or more time consuming) to obtain than the unlabeled data. Consider a dataset with one variable of interest, sometimes referred to as the target variable or outcome variable, and one or more covariates. The difficulty for collecting labeled data stems from collecting the target variable, whereas unlabeled data containing only the covariates, that is with the target variable missing, can be easily collected. Semi-supervised learning focuses on uncovering the (nonlinear) relation between the target variable and the covariates. Estimations and predictions based on such relations and using the additional unlabeled data often show substantially improved performance. For example, for classification analysis see Vapnik (2013) and Zhu and Goldberg (2009); for regression analysis see Wasserman and Lafferty (2008), Azriel et al. (2022), and Chakrabortty and Cai (2018).

Semi-supervised inference aims at estimating parameters or quantities regarding the target variable in the semi-supervised model. Zhang, Brown, and Cai (2019) investigates the general semi-supervised framework and shows how to use the unlabeled data to improve the estimation of the mean of the target variable; for inference on heavy tailed distributions in this framework, see Ahmed and Einmahl (2019).

Extreme value statistics deals with estimation of parameters or quantities related to the tail of a distribution, only making semi-parametric assumptions on this tail. Consequently, most of extreme value methods start with a relatively large number of observations n, but select only $k\ll n$ extreme observations from the full sample for statistical inference. Two techniques are often used in selecting the extreme observations: the peaks-over-threshold (POT) approach which selects the highest k observations, and the block maxima (BM) approach which splits the full sample into k blocks and selects the maxima of each block. Since only k observations are used in estimation, typically consistent estimators have a speed of convergence of $1/\sqrt{k}$. In practice, to obtain accurate estimators for tail parameters/quantities, one needs a sample with a relatively large sample size n to guarantee a sufficient number of extreme observations. In contrast, the semi-supervised model is greatly suitable for statistics of extremes in case data on the target variable are hard to obtain.

The main goal of this article is 2-fold. First, we derive in this semi-supervised setting a new, improved pseudo-maximum likelihood estimator (MLE) for a general extreme value index γ and establish its asymptotic behavior. This extreme value index describes the tail heaviness of a probability distribution. If $\gamma >0$ the distribution is heavy tailed and has an infinite right endpoint, if γ = 0 the distribution is light tailed and may have an infinite or finite endpoint, and if $\gamma <0$ the endpoint is finite, see, for example, Beirlant et al. (2004) or de Haan and Ferreira (2006) for a thorough treatment of extreme value theory and the corresponding statistical inference. Second, we use the adapted estimator of γ, together with an adapted estimator of the so-called scale, to obtain an improved estimator of an extreme quantile in the semi-supervised model. We establish the asymptotic behavior of this adapted extreme quantile estimator. Simulations show that the improvement for the new extreme quantile estimator is actually larger than that for the new estimator of γ.

For ease of explanation of our novel estimator of the extreme value index, let us assume that there is only one covariate. We estimate γ for the variable of interest initially (that is ignoring the covariate) using the pseudo-MLE $\widehat{\gamma }$, see Smith (1987) and Drees, Ferreira, and de Haan (2004). Then, we choose a number g and for the covariate we transform the labeled data empirically (using all the labeled and unlabeled data) such that they obtain an artificial extreme value index g. Using the transformed covariates of the labeled data, we estimate the known g by the pseudo-MLE $\widehat{g}$, say, and use the difference $\widehat{g}-g$ to adapt and substantially improve the initial estimator $\widehat{\gamma }$ for the extreme value index γ of the variable of interest. For this adaptation the tail dependence between the target variable and the covariate is crucial.

The specific contributions of this article are as follows. First, we provide a general result for the relevant, nonstandard tail quantile process (see Lemma 6.5). Based on this tail quantile process result, one could improve other estimators based on the POT approach in the semi-supervised model. Also, we impose no assumptions on the tail of the covariates. When analyzing the tail of the target variable, it is crucial to assume regularity in its tail such as the max-domain of attraction condition in extreme value analysis. However, requiring such conditions for the covariates can be restrictive in applications. Finally, our main results are valid for a broad class of distributions for the target variable: we deal with a general extreme value index γ, which can be positive, zero, or negative. This is particularly important for applications where the sign of γ is not known beforehand. For example, when analyzing extreme weather, various studies find that the extreme value index is around zero for different meteorologic variables: for hourly surge level on the English east coast (Coles and Tawn 1991), for hourly maximum wind speed in Sheffield, UK (Coles and Walshaw 1994), for wave height and still water level on the Dutch coast (de Haan and de Ronde 1998) and for daily rainfall in North Holland, The Netherlands (Buishand, de Haan, and Zhou 2008).

This article is organized as follows. In Section 2, for clearness of the exposition, we first introduce our adapted estimator for the extreme value index in the semi-supervised model with one covariate and we establish its asymptotic normality. Subsequently we introduce the adapted estimator of an extreme quantile and analyze its asymptotic behavior. In Section 3 we consider the general multivariate semi-supervised setting and present and establish asymptotic normality of the adapted estimator of the extreme value index. Section 4 is devoted to a simulation study. The improved performance, in terms of variance, of the adapted estimators compared with the initial estimators is shown. An application to rainfall in France can be found in Section 5 and the detailed proofs are deferred to Section 6 and the supplementary material.

2 Main Results: One Covariate

2.1 Estimation of the Extreme Value Index

Let F be a bivariate distribution with marginals F₁ and F₂. We assume that F₁ is in the max-domain of attraction of an extreme value distribution $G_{\gamma }$, where γ is the extreme value index, our parameter of interest. Let the pairs $(X_1,Y_1),\dots ,(X_n,Y_n)$ be a random sample from F, and let $(Y_{n+1},\dots ,Y_{n+m})$ be a random sample from F₂, independent from the n pairs. This is the semi-supervised model. Assume that the tail copula R of (X₁, Y ₁) exists: (1) $$\begin{array}{c}R(x,y)=\underset{t\downarrow 0}{\mathrm{\text{lim}}}\frac{1}{t}\mathbb{P}\left(1-F_1(X_1)\le \mathit{\text{tx}},1-F_2(Y_1)\le \mathit{\text{ty}}\right),\\ (x,y)\in {[0,\infty ]}^2\setminus \{(\infty ,\infty )\}.\end{array}$$(1)

Denote the order statistics of $X_i,i=1,\dots ,n,$ with $X_{1:n}\le \cdots \le X_{n:n},$ and similarly for the $Y_i,i=1,\dots ,n.$ We estimate $\gamma >-\frac{1}{2}$ with the often used pseudo-MLE $\widehat{\gamma }$ based on $X_{n-k:n},\dots ,X_{n:n}$, for $k\in \{1,\dots ,n-1\}$; see sec. 3.4 in de Haan and Ferreira (2006).

Define for $i=1,\dots ,n$, (2) $${\tilde{Y}}_i=\{\begin{array}{cc}\frac{{\left(1-\left(F_{n+m}(Y_i)-\frac{1}{2(n+m)}\right)\right)}^{-g}-1}{g}& ,\hspace{1em}g\ne 0,\\ -\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(1-\left(F_{n+m}(Y_i)-\frac{1}{2(n+m)}\right)\right)& ,\hspace{1em}g=0,\end{array}\hspace{1em}$$(2) where $F_{n+m}$ is the empirical distribution function based on $Y_l,l=1,\dots ,n+m$, and $g>-\frac{1}{2}$ is a number we may choose that mimics an extreme value index. Let the order statistics of ${\tilde{Y}}_i,i=1,\dots ,n,$ be denoted by ${\tilde{Y}}_{1:n}\le \dots \le {\tilde{Y}}_{n:n},$ and let $\widehat{g}$ be the pseudo-MLE of g based on ${\tilde{Y}}_{n-k:n},\dots ,{\tilde{Y}}_{n:n}$, using the same k as before. Of course, since we choose and hence know g, there is no direct need to estimate it. We will show below, however, that the dependence of the difference $\widehat{g}-g$ and $\widehat{\gamma }$, helps to improve the estimator of γ in the semi-supervised setting.

For the asymptotic theory, we assume that $m=m(n)$ and (3) $$k\to \infty ,\frac{k}{n}\to 0,\sqrt{\frac{n}{n+m}}\to \nu \in (0,1),\hspace{1em}\text{as }n\to \infty .$$(3)

We begin with establishing the joint asymptotic normality of $\widehat{\gamma }$ and $\widehat{g}$, a crucial result for deriving and showing asymptotic normality of our semi-supervised estimator (SSE) of γ. For that purpose we need the usual second order condition on the marginal distribution F₁. Let $U_1=F_1^{-1}(1-1/·)$ be the tail quantile corresponding to F₁. We assume that there exist a positive scale function a, a positive or negative function A, with $\underset{t\to \infty }{\mathrm{\text{lim}}}A(t)=0$, and $\rho \le 0$, such that for $x>0,$ (4) $$\underset{t\to \infty }{\mathrm{\text{lim}}}\frac{\frac{U_1(\mathit{\text{tx}})-U_1(t)}{a(t)}-\frac{x^{\gamma }-1}{\gamma }}{A(t)}=\Psi (x),\hspace{1em}\gamma \in \mathbb{R},$$(4) where $$\Psi (x)=\{\begin{array}{cc}\frac{x^{\gamma +\rho }-1}{\gamma +\rho },& \rho <0,\\ \frac{1}{\gamma }{x}^{\gamma }\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}x,& \mathit{\hspace{1em}\gamma }\ne \rho =0,\\ \frac{1}{2}{\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}^{2}x,& \gamma =\rho =0,\end{array}$$ see de Haan and Ferreira (2006), p. 46.

Proposition 2.1.

Assume $\gamma >-\frac{1}{2}$ and choose $g>-\frac{1}{2}$. Assume that F₂ is continuous, (1, 3), and (4) hold, and $\sqrt{k}A\left(\frac{n}{k}\right)\to \lambda \in \mathbb{R},$ as $n\to \infty$, then with probability tending to 1, there exist unique maximizers of the likelihood functions based on ${\left\{X_i\right\}}_{i=1}^n$ and ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$, denoted as $(\widehat{\gamma },\widehat{g})$, such that $$\left(\mathrm{}\sqrt{k}\left(\mathrm{}\widehat{\gamma }-\gamma \mathrm{}\right),\sqrt{k}\left(\widehat{g}-g\right)\mathrm{}\mathrm{}\right)\stackrel{d}{\to }N\left(\mathrm{}\left[\mathrm{}\frac{\lambda (1+\gamma )}{(1-\rho )(1+\gamma -\rho )},0\mathrm{}\right],\Sigma \mathrm{}\mathrm{}\right)$$where $$\Sigma =\left[\mathrm{}\mathrm{}\begin{array}{cc}{(1+\gamma )}^2& (1-{\nu }^2)(1+\gamma )(1+g)R_g\\ (1-{\nu }^2)(1+\gamma )(1+g)R_g& (1-{\nu }^2){(1+g)}^2\end{array}\mathrm{}\mathrm{}\right],$$ with (5) $$\begin{array}{c}{R}_g=R(1,1)+\frac{g-\gamma }{\gamma +g+1}\\ \hspace{1em}\times \left((2\gamma +1){\int }_0^1\frac{R(s,1)}{s^{1-\gamma }}\mathit{\text{ds}}-(2g+1){\int }_0^1\frac{R(1,t)}{t^{1-g}}\mathit{\text{dt}}\right).\end{array}$$(5)

Based on Proposition 2.1, we derive the SSE of γ. For this derivation only, take $\lambda =0$. Then the approximate bivariate normal distribution of $(\widehat{\gamma },\widehat{g}-g),$ has mean $[\gamma ,0]$ and estimated covariance matrix $$\frac{1}{k}\widehat{\Sigma }=\frac{1}{k}\left[\text{}\begin{array}{cc}{(1+\widehat{\gamma })}^2& (1-\frac{n}{n+m})(1+\widehat{\gamma })(1+g){\widehat{R}}_g\\ (1-\frac{n}{n+m})(1+\widehat{\gamma })(1+g){\widehat{R}}_g& (1-\frac{n}{n+m}){(1+g)}^2\end{array}\text{}\right],$$where ${\widehat{R}}_g$ is the estimator of $R_g,$ obtained by replacing γ with $\widehat{\gamma }$ and the tail copula R, at 3 places, with its natural estimator (6) $$\widehat{R}(x,y)=\frac{1}{k}\sum _{i=1}^n{1}_{[X_i\ge X_{n-\left[\mathit{\text{kx}}\right]+1:n},Y_i\ge Y_{n-\left[\mathit{\text{ky}}\right]+1:n}]},\hspace{1em}x,y\ge 0,$$(6) see, for example, Drees and Huang (1998). Maximizing the thus obtained approximate likelihood function of the single “data point” $(\widehat{\gamma },\widehat{g}-g)$ with respect to the unknown γ we obtain as SSE for γ: (7) $${\widehat{\gamma }}_g=\widehat{\gamma }-\frac{1+\widehat{\gamma }}{1+g}{\widehat{R}}_g(\widehat{g}-g).$$(7)

Now we present the main result of this subsection, the asymptotic normality of the SSE of γ.

Theorem 2.1.

Under the conditions of Proposition 2.1, as $n\to \infty$, (8) $$\begin{array}{c}\sqrt{k}({\widehat{\gamma }}_g-\gamma )\\ \hspace{1em}\stackrel{d}{\to }N\left(\frac{\lambda (1+\gamma )}{(1-\rho )(1+\gamma -\rho )},{(1+\gamma )}^2\left[1-(1-{\nu }^2)R_g^2\right]\right).\end{array}$$(8)

Remark 2.1.

Note that the asymptotic bias of the SSE ${\widehat{\gamma }}_g$ is the same as that of the pseudo-MLE $\widehat{\gamma }$ (in Proposition 2.1). Therefore, when comparing both estimators we can and will focus on the (relative) reduction of the asymptotic variance which is equal to $(1-{\nu }^2)R_g^2$. The value of the crucial $R_g\in [-1,1]$ depends on the known g and the unknown γ and R. Note that R_g can indeed be positive, zero, or negative and R_g can exceed R(1, 1) even when R is symmetric in its arguments. Nevertheless it is appealing to consider $g=\gamma$, reducing R_g to R(1). Since γ is unknown, this would lead to the choice $g=\widehat{\gamma }$. We will not pursue this, however, since simulation results show that this random g leads to much smaller variance reductions than when taking a deterministic g not too far away from γ. Also observe that when g is close to γ and R is symmetric, then $R_g-R(1,1)$ is of order ${(g-\gamma )}^2$, that is, the variance reduction does not change much with the choice of g and is close to R(1). We will see in Section 4 through simulations that the variance reduction is substantial in “standard” semi-supervised settings.

Finally, note that in general if the function R is identical to 0 (tail independence), ${\widehat{\gamma }}_g$ has the same asymptotic behavior as $\widehat{\gamma }$. If $R(x,y)>0$, for some $x,y>0$ (tail dependence), then $R(1,1)>0$. In that case for $g=\gamma$ we obtain $R_g=R(1,1)$: the asymptotic variance of ${\widehat{\gamma }}_g$ is smaller than that of $\widehat{\gamma }$ and it decreases with R(1). In case of tail independence, simulations indicate that in finite-samples there can be a variance inflation up to about 10%; see Section D in the supplementary material. Therefore, it is recommended to apply the method in practice only when the data clearly exhibit tail dependence.

2.2 Estimation of an Extreme Quantile

Let p be a very small positive number. For the asymptotic theory this means that $p=p(n)$ with $\underset{n\to \infty }{\mathrm{\text{lim}}}p(n)=0$. Then the extreme quantile x_p is defined by $x_p=F_1^{-1}(1-p)=U_1(1/p).$ We first estimate γ and the scale $a(n/k)$ with the pseudo-MLEs $\widehat{\gamma }$ (as above) and $\widehat{\sigma }$, based on $X_{n-k:n},\dots ,X_{n:n}$, see, for example, sec. 3.4 in de Haan and Ferreira (2006). The natural estimator of the extreme quantile x_p is then given by (see Dekkers, Einmahl, and de Haan 1989; de Haan and Rootzén 1993) (9) $${\widehat{x}}_p=X_{n-k:n}+\widehat{\sigma } \frac{{\left(\frac{k}{\mathit{\text{np}}}\right)}^{\widehat{\gamma }}-1}{\widehat{\gamma }}.$$(9)

Using (2), we see that g and ${(\frac{n}{k})}^g$ are numbers that mimic an extreme value index and a scale. They are estimated with $\widehat{g}$ and ${\tilde{\sigma }}_g$ using the pseudo-MLEs based on ${\tilde{Y}}_{n-k:n},\dots ,{\tilde{Y}}_{n:n}$, for the same k as above. Similar to the derivation of ${\widehat{\gamma }}_g$, we find an adapted estimator of σ based on the joint asymptotic normality of $\widehat{\sigma }$ and ${\tilde{\sigma }}_g$. For this, write for $g\ne 0$ $$\begin{array}{c}{S}_g=(\gamma +2)(g+2)R(1,1)-\frac{g(\gamma +1)}{\gamma }{\int }_0^1\frac{R(s,1)}{s}\mathit{\text{ds}}\\ -\frac{\gamma (g+1)}{g}{\int }_0^1\frac{R(1,t)}{t}\mathit{\text{dt}}+\frac{{(\gamma +1)}^2(g+1)(2\gamma +1)}{\gamma g}\\ (\frac{1}{g+1}-\frac{1}{\gamma +1}+(g+1)\frac{g-\gamma }{\gamma +g+1}){\int }_0^1\frac{R(s,1)}{s^{1-\gamma }}\mathit{\text{ds}}\\ +\frac{(\gamma +1){(g+1)}^2(2g+1)}{\gamma g}(\frac{1}{\gamma +1}-\frac{1}{g+1}\\ +(\gamma +1)\frac{\gamma -g}{\gamma +g+1}){\int }_0^1\frac{R(1,t)}{t^{1-g}}\mathit{\text{dt}},\\ S_0=2(\gamma +2)R(1,1)+(3\gamma -1){\int }_0^1\frac{R(1,t)}{t}\mathit{\text{dt}}\\ +\gamma {\int }_0^1\frac{R(1,t)}{t}\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}t \mathit{\text{dt}}-{(2\gamma +1)}^2{\int }_0^1\frac{R(s,1)}{s^{1-\gamma }}\mathit{\text{ds}},\end{array}$$and estimate it with ${\widehat{S}}_g$ obtained by replacing γ with $\widehat{\gamma }$ and R with $\widehat{R}$. Now the SSE for the scale σ is defined by $${\widehat{\sigma }}_g=\widehat{\sigma }\left(1-\frac{{\widehat{S}}_g}{1+{(1+g)}^2}\left(\frac{{\tilde{\sigma }}_g}{{(\frac{n}{k})}^g}-1\right)\right).$$

The adapted version of the extreme quantile estimator in (9) is now obtained by plugging in the SSEs for the extreme value index and the scale instead of the pseudo-MLEs: (10) $${\widehat{x}}_{p_g}=X_{n-k:n}+{\widehat{\sigma }}_g\frac{{\left(\frac{k}{\mathit{\text{np}}}\right)}^{{\widehat{\gamma }}_g}-1}{{\widehat{\gamma }}_g}.$$(10)

In order to establish its asymptotic normality, we first present a theorem which states the joint asymptotic normality of the three random components involved in this definition. The proofs of the results in this section are deferred to the supplementary material, Section C.

Theorem 2.2.

Assume $\gamma >-\frac{1}{2}$ and choose $g>-\frac{1}{2}$. Assume that F₂ is continuous, (1, 3), and (4) hold, and $\sqrt{k}A\left(\frac{n}{k}\right)\to \lambda \in \mathbb{R},$ as $n\to \infty .$ Then as $n\to \infty ,$ (11) $$\begin{array}{c}\sqrt{k}\left({\widehat{\gamma }}_g-\gamma ,\frac{{\widehat{\sigma }}_g}{a\left(\frac{n}{k}\right)}-1,\frac{X_{n-k:n}-U_1\mathrm{}\left(\frac{n}{k}\right)}{a\left(\frac{n}{k}\right)}\right)\\ \hspace{1em}\stackrel{d}{\to }N\left(\left[\frac{\lambda (\gamma +1)}{(1-\rho )(1+\gamma -\rho )},\frac{-\rho \lambda }{(1-\rho )(1+\gamma -\rho )},0\right],K\right),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\end{array}$$(11) $$K=\left[\begin{array}{ccc}{(1+\gamma )}^2\left[1-(1-{\nu }^2)R_g^2\right]& -(1+\gamma )[1+(1-{\nu }^2)Q]& (1-{\nu }^2)M_1\\ -(1+\gamma )[1+(1-{\nu }^2)Q]& 1+{(1+\gamma )}^2-(1-{\nu }^2)\frac{S_g^2}{1+{(1+g)}^2}& \gamma +(1-{\nu }^2)M_2\\ (1-{\nu }^2)M_1& \gamma +(1-{\nu }^2)M_2& 1\end{array}\right],$$where $Q=\frac{R_g{S}_g}{1+{(1+g)}^2}+\frac{1+\gamma }{1+g}{R}_g{Q}_1+\frac{S_g}{1+{(1+g)}^2}{Q}_2,$ R_g is as in (5), $$\begin{array}{c}{Q}_1:={(1+g)}^2[\frac{-\gamma -2}{(g+1)(\gamma +1)}R(1,1)+\frac{\gamma }{\gamma +1}{\int }_0^1\frac{R(1,t)}{t}\frac{1-t^g}{g}\mathit{\text{dt}}\\ \hspace{1em}+(\frac{2(g-\gamma )((\gamma +1)(g+1)+g)}{(\gamma +1)(g+1)(\gamma +g+1)}+\frac{2g+1}{(g+1)(\gamma +g+1)})\\ \hspace{1em}\times {\int }_0^1\frac{R(1,t)}{t^{1-g}}\mathit{\text{dt}}+\frac{(2\gamma +1)(\gamma -g)}{(\gamma +g+1)(g+1)}{\int }_0^1\frac{R(s,1)}{s^{1-\gamma }}\mathit{\text{ds}}],\\ Q_2:=(1+\gamma )(1+g)[\frac{-g-2}{(g+1)(\gamma +1)}R(1,1)+\frac{g}{g+1}{\int }_0^1\frac{R(s,1)}{s}\frac{1-s^{\gamma }}{\gamma }\mathit{\text{ds}}\\ \hspace{1em}+(\frac{2(\gamma -g)((\gamma +1)(g+1)+\gamma )}{(\gamma +1)(g+1)(\gamma +g+1)}+\frac{2\gamma +1}{(\gamma +1)(\gamma +g+1)})\\ \hspace{1em}\times {\int }_0^1\frac{R(s,1)}{s^{1-\gamma }}\mathit{\text{ds}}+\frac{(2g+1)(g-\gamma )}{(\gamma +g+1)(\gamma +1)}{\int }_0^1\frac{R(1,t)}{t^{1-g}}\mathit{\text{dt}}],\\ M_1:=(1+\gamma )(1+g)R_g\\ \hspace{1em}\times [2{\int }_0^1\frac{R(1,s)}{s^{1-g}}\mathit{\text{ds}}-{\int }_0^1\frac{R(1,s)}{s}\frac{1-s^g}{g}\mathit{\text{ds}}-\frac{1}{g+1}R(1,1)],\\ M_2:=\frac{-(1+g)}{1+{(1+g)}^2}{S}_g\\ \hspace{1em}\times [(2g+3){\int }_0^1\frac{R(1,s)}{s^{1-g}}\mathit{\text{ds}}-{\int }_0^1\frac{R(1,s)}{s}\frac{1-s^g}{g}\mathit{\text{ds}}-\frac{g+2}{g+1}R(1,1)].\end{array}$$

In the expression of $Q_1,Q_2,M_1$, and M₂, when considering g = 0 or γ = 0, the terms such as $\frac{1-t^g}{g},\frac{1-s^{\gamma }}{\gamma }$ should be read as their corresponding limit.

From Theorem 2.2, we obtain the main result in this section, the asymptotic normality of the adapted extreme quantile estimator. Define $q_{\gamma }(t)={\int }_1^t{s}^{\gamma -1}\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s \mathit{\text{ds}}$, for t > 1.

Theorem 2.3.

Assume $\gamma >-\frac{1}{2}$ and choose $g>-\frac{1}{2}$. Assume that F₂ is continuous, (1, 3), and (4) hold, the second-order parameter ρ is negative, or zero with γ negative, and $\sqrt{k}A\left(\frac{n}{k}\right)\to \lambda \in \mathbb{R},$ as $n\to \infty$. Further assume $p=o(k/n)$ and $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(\mathit{\text{np}})=o\left(\sqrt{k}\right),$ as $n\to \infty .$ Then as $n\to \infty ,$ (12) $$\sqrt{k}\frac{{\widehat{x}}_{p_g}-x_p}{a\left(\frac{n}{k}\right)q_{\gamma }\mathrm{}\left(\frac{k}{\mathit{\text{np}}}\right)}\stackrel{d}{\to }N\left(\frac{\lambda (1+\gamma )}{(1-\rho )(\gamma -\rho +1)}b,v\right)$$(12) where b = 1 and $v={(1+\gamma )}^2\left[1-(1-{\nu }^2)R_g^2\right]$, for $\gamma \ge 0>\rho$, and $b=\rho (1+2\gamma )/(\gamma +\rho )$ and $v={(1+\gamma )}^2(1+2\gamma )-(1-{\nu }^2)[{(1+\gamma )}^2{R}_g^2+{\gamma }^2\frac{S_g^2}{1+{(1+g)}^2}-2\gamma (1+\gamma )Q-2{\gamma }^2{M}_1+2{\gamma }^3{M}_2]$, for $\gamma <0$.

Remark 2.2.

Note that the asymptotic bias is the same as that of the standard extreme quantile estimator ${\widehat{x}}_p$ in (9), whereas the asymptotic variance is different. In case $\gamma \ge 0$, the relative variance reduction when using the adapted estimator is $(1-{\nu }^2)R_g^2\ge 0$, which is the same as the relative variance reduction of the adapted extreme value index estimator. For $\gamma <0$ it is not clear if the variance reduction is nonnegative, but in the simulation study below we show substantial variance reductions, irrespective of the sign of γ. A joint variance reduction approach, instead of improving $\widehat{\gamma }$ and $\widehat{\sigma }$ separately, would lead to a nonnegative variance reduction for all $\gamma >-\frac{1}{2}$, but this is beyond the scope of this article.

3 Main Results: Multiple Covariates

In this section we consider the more general situation with d–1 covariates where d > 2. Consider a d-variate distribution F, with marginals $F_1,\dots ,F_d$. We assume again that (only) F₁ is in the max-domain of attraction of an extreme value distribution $G_{\gamma }$. Let $F_{-}$ be the distribution function of the last d–1 components of a random vector with distribution function F. Let $(X_1,Y_{1,2},\dots ,Y_{1,d}),\dots ,(X_n,Y_{n,2},\dots ,Y_{n,d})$ be a random sample of size n from F and let $(Y_{n+1,2},\dots ,Y_{n+1,d}),\dots ,(Y_{n+m,2},\dots ,Y_{n+m,d})$ be a random sample of size m from $F_{-}$, independent of the d-variate random sample of size n. This is the multivariate semi-supervised setting.

Then, for fixed $j=2,\dots ,d$, we use all data for the covariates ${\left\{Y_{i,j}\right\}}_{i=1}^{n+m}$ to obtain ${\left\{{\tilde{Y}}_{i,j}\right\}}_{i=1}^n$ as in (2), where we may choose a number $g>-\frac{1}{2},$ that mimics an extreme value index, as before. For $k\in \{1,\dots ,n-1\}$, let, similarly as in the previous section, $\widehat{\gamma }$ and ${\widehat{g}}_j,j=2,\dots ,d$, be the pseudo-MLEs of γ and (d–1 times) of g, respectively. Assume the existence of the tail copula R_ij of the ith and the jth component, (13) $$R_{\mathit{\text{ij}}}(x,y)=\underset{t\downarrow 0}{\mathrm{\text{lim}}}\frac{1}{t}\mathbb{P}\left(1-F_i(Y_{1,i})\le \mathit{\text{tx}},1-F_j(Y_{1,j})\le \mathit{\text{ty}}\right),$$(13)

where $(x,y)\in {[0,\infty ]}^2\setminus \{(\infty ,\infty )\},1\le i,j\le d.$ Here $Y_{1,1}$ is understood as X₁. Again, we first consider the joint asymptotic normality of $\widehat{\gamma },$ and ${\widehat{g}}_j,j=2,\dots ,d.$

Proposition 3.1.

Assume $\gamma >-\frac{1}{2}$ and choose $g>-\frac{1}{2}$. Assume that $F_j,j=2,\dots ,d$, is continuous, (3, 4), and (13) hold, and as $n\to \infty ,$ $\sqrt{k}A\left(\frac{n}{k}\right)\to \lambda \in \mathbb{R},$ then with probability tending to 1, there exist unique maximizers of the likelihood functions based on ${\left\{X_i\right\}}_{i=1}^n,{\left\{{\tilde{Y}}_{i,2}\right\}}_{i=1}^n,\dots ,{\left\{{\tilde{Y}}_{i,d}\right\}}_{i=1}^n$, denoted as $(\widehat{\gamma },{\widehat{g}}_2,\dots ,{\widehat{g}}_d)$, such that $$\begin{array}{c}\left(\sqrt{k}(\widehat{\gamma }-\gamma ),\sqrt{k}({\widehat{g}}_2-g),\dots ,\sqrt{k}({\widehat{g}}_d-g)\right)\\ \hspace{1em}\stackrel{d}{\to }N\left(\left[\frac{\lambda (\gamma +1)}{(1-\rho )(1+\gamma -\rho )},0,\dots ,0\right],{\Sigma }_d\right),\end{array}$$with ${\Sigma }_d=\Gamma {\Gamma }^T°H$ (“$°$” is the Hadamard or entrywise product), where $$\Gamma =\left[\begin{array}{c}1+\gamma \\ 1+g\\ .\\ .\\ 1+g\end{array}\right],\hspace{1em}H=\left[\begin{array}{ccccc}1& h_{12}& .& .& h_{1d}\\ h_{12}& 1-{\nu }^2& .& .& h_{2d}\\ .& .& .& & .\\ .& .& & .& .\\ h_{1d}& h_{2d}& .& .& 1-{\nu }^2\end{array}\right],$$

$h_{1i}=(1-{\nu }^2)[R_{1i}(1,1)+\frac{g-\gamma }{\gamma +g+1}[(2\gamma +1){\int }_0^1\frac{R_{1i}(s,1)}{s^{1-\gamma }}\mathit{\text{ds}}$ $-(2g+1){\int }_0^1\frac{R_{1i}(1,t)}{t^{1-g}}\mathit{\text{dt}}]],$ and $h_{\mathit{\text{ij}}}=(1-{\nu }^2)R_{\mathit{\text{ij}}}(1,1),\hspace{1em}i=2,\dots ,d,j=i+1,\dots ,d.$

Very similar to the bivariate case, let λ = 0 and derive the SSE of γ by using the approximate multivariate normal distribution of $(\widehat{\gamma },{\widehat{g}}_2-g,\dots ,{\widehat{g}}_d-g),$ with mean $[\gamma ,0,\dots ,0]$, and variance $\frac{1}{k}{\widehat{\Sigma }}_d=\frac{1}{k}\widehat{\Gamma }{\widehat{\Gamma }}^T°\widehat{H},$ where for the estimation of $R_{\mathit{\text{ij}}},$ ${\widehat{R}}_{\mathit{\text{ij}}}$ is defined like in (6). By maximizing the approximate likelihood function of $(\widehat{\gamma },({\widehat{g}}_2-g),\dots ,({\widehat{g}}_d-g))$ with respect to γ, we obtain the SSE in this multivariate setting: (14) $${\widehat{\gamma }}_g=\widehat{\gamma }+\frac{1+\widehat{\gamma }}{1+g}\sum _{j=2}^d\frac{{\widehat{H}}_{1j}^{-1}}{{\widehat{H}}_{11}^{-1}}({\widehat{g}}_j-g),$$(14) where ${\widehat{H}}_{\mathit{\text{ij}}}^{-1}$ is the entry in the ith row and jth column of the inverse of the matrix $\widehat{H}$. The following theorem shows the asymptotic behavior of the improved estimator ${\widehat{\gamma }}_g$.

Theorem 3.1.

Assume that H is invertible. Then under the conditions of Proposition 3.1, as $n\to \infty ,$ (15) $$\sqrt{k}({\widehat{\gamma }}_g-\gamma )\stackrel{d}{\to }N\left(\frac{\lambda (1+\gamma )}{(1-\rho )(1+\gamma -\rho )},{\sigma }^2\right),$$(15) where $$\begin{array}{c}{\sigma }^2={(1+\gamma )}^2(1+\frac{1}{{(H_{11}^{-1})}^2}[2\sum _{i=1}^d\sum _{j=i+1}^d{H}_{1i}^{-1}{H}_{1j}^{-1}{h}_{\mathit{\text{ij}}}\\ +(1-{\nu }^2)\sum _{j=2}^d{(H_{1j}^{-1})}^2]).\end{array}$$

Remark 3.1.

It can be shown that increasing the number of covariates leads to a lower (or the same) asymptotic variance. More specifically, the quality of the estimation improves if we replace the “one-covariate” estimator ${\widehat{\gamma }}_g$ in (7) by the present one. In particular, if we consider d = 3 and for simplicity $g=\gamma$ and $R_{12}(1,1)=R_{13}(1,1)$, the asymptotic relative variance reduction in this case is equal to $2(1-{\nu }^2)R_{12}^2(1,1)/(1+R_{23}(1,1))$, see the formula for the asymptotic relative variance reduction in Ahmed and Einmahl (2019), which applies in this case. Consequently, the maximal variance reduction is achieved when $R_{23}(1,1)$ takes its minimal attainable value. Decreasing the tail dependence of the covariates will make them, as a pair, more relevant and hence will make the variance reduction larger. See also the simulations in the trivariate case at the end of the next section.

4 Simulation Study

Here we perform mainly for the one-covariate setting a detailed simulation study, but at the end of this section we briefly consider two covariates.

First we investigate the finite sample performance of our novel SSE of γ and then we compare in detail the variances of the SSE with those of the pseudo-MLE based on $X_{n-k:n},\dots ,X_{n:n}$ only. In addition, we compare the variances of the adapted extreme quantile estimator with those of the standard estimator in (9).

We begin with simulating data from the bivariate Cauchy distribution restricted to the first quadrant. This Cauchy density is proportional to ${(1+x{S}^{-1}{x}^T)}^{-3/2}$ where S is a 2 × 2 scale matrix with 1 on the diagonal and s off-diagonal. For s we take alternatingly two values: 0 and 0.8. For the value 0.8 we have a rather strong tail dependence and for s = 0 the tail dependence is somewhat weaker. More precisely, the value 0 corresponds to an R(1)-value of 0.59, whereas for s = 0.8 this value is 0.76. We expect that stronger tail dependence leads to a larger variance reduction. These data are denoted by $(\check{X_i},Y_i)$. To obtain our data (X_i, Y _i), where the X_i have extreme value index γ, we transform the ${\check{X}}_i$, as follows (16) $$X_i=\{\begin{array}{cc}\frac{{(1-F_s({\check{X}}_i))}^{-\gamma }-1}{\gamma },& \mathit{\hspace{1em}\gamma }\ne 0,\\ -\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1-F_s({\check{X}}_i)),& \mathit{\hspace{1em}\gamma }=0,\end{array}$$(16) where F_s is the distribution function of $\check{X_i}$. Simulations are performed for values of γ that are negative, positive or 0.

First, we generate 500 samples of sizes $n=500,m=1000$, for s = 0.8, and estimate γ using the SSE and the pseudo-MLE for $n=500,m=1000$. We depict the root mean squared error (RMSE) based on these 500 samples as a function of k. We consider $\gamma =-0.25,0,$ and, $0.25,$ and take g = 0. The RMSE of the SSE (indicated by AMLE in Figure 1) is indeed substantially lower than that of the pseudo-MLE for the different values of $\gamma .$

Fig. 1 RMSE using the pseudo-MLE and the SSE-MLE. From left to right: $\gamma =-0.25,0,0.25$.

Next, we focus on the (relative) variance reduction of the SSE in comparison to the pseudo-MLE. We use the following values of n and m (and k):

• $n=1000$, m = 500 (less unlabeled than labeled data) and k = 250,

• $n=1000$, m = 1000 (equal number of unlabeled and labeled data) and k = 250,

• $n=500$, m = 1000 (more unlabeled than labeled data) and k = 125.

Table 1 shows the empirical percentages of variance reduction for different values of γ and g. The results are based on $10,000$ replications. We observe that the variance reduction ranges from 10% to more than 30%, hence, indeed the SSE has a substantially smaller variance than the pseudo-MLE. By comparing the three panels, we observe that the variance reduction increases substantially with the ratio of the number of unlabeled data m and the number of labeled data n, which is in line with the asymptotic theory. Observe that the actual choice of g does not have a large influence as long as it is somewhat close to γ, a choice that is in practice often feasible.

Table 1 Variance reduction for different extreme value indices.

(i) $(n,m)=(1000,500)$
g	γ				g	γ
	–0.3	–0.2	–0.1	0		0	0.1	0.2	0.3
–0.25	13.8%	14.0%	13.4%	12.0%	0	15.4%	15.6%	15.4%	14.8%
–0.125	12.6%	14.2%	14.8%	14.1%	0.125	15.3%	15.9%	16.1%	15.8%
0	10.4%	13.0%	14.6%	15.4%	0.25	14.7%	15.5%	16.1%	16.4%
(ii) $(n,m)=(1000,1000)$
	–0.3	–0.2	–0.1	0		0	0.1	0.2	0.3
–0.25	20.2%	21.1%	20.5%	18.8%	0	23.6%	23.5%	23.1%	22.2%
–0.125	18.9%	21.5%	22.4%	22.0%	0.125	23.4%	24.1%	24.3%	25.0%
0	16.3%	19.9%	22.3%	23.6%	0.25	22.4%	23.7%	24.6%	24.6%
(iii) $(n,m)=(500,1000)$
	–0.3	–0.2	–0.1	0		0	0.1	0.2	0.3
–0.25	22.4%	25.8%	25.7%	24.5%	0	30.8%	30.2%	29.9%	29.4%
–0.125	21.6%	26.2%	28.0%	28.0%	0.125	31.0%	30.4%	31.0%	30.8%
0	18.6%	24.6%	28.1%	30.8%	0.25	30.0%	30.1%	31.4%	31.5%

Next we investigate in more detail the sensitivity of the variance reduction to the choice of g using a wider range of values of g, including cases where γ and g have opposite sign. In these simulations, we take s = 0 for the bivariate Cauchy distribution. The results, based on $10,000$ replications, for the aforementioned values of n, m and k are presented in Figure 2. Generally, there is always variance reduction, but for $|g-\gamma |$ relatively large the reduction is lower than when g is closer to γ. Observe that for the present range of γ (–0.3 to 0.3) the choice g = 0 yields an almost maximal variance reduction.

Fig. 2 Variance reduction for various combinations of γ and g.

Finally we study in more detail the effect of the size of m, the number of unlabeled data, on the variance reduction; again we take s = 0. We consider the case where n = 500 and let m vary; we choose g = 0. The results are based on 500 replications. Table 2 shows that the variance reduction approximately doubles when m ranges from 500 to $10,000$.

Table 2 Variance reduction for different numbers of unlabeled data m.

m	γ
	–0.3	–0.2	–0.1	0	0.1	0.2	0.3
500	6.0%	8.3%	9.6%	11.4%	11.7%	11.9%	11.8%
1000	8.7%	12.3%	14.4%	15.6%	12.9%	13.0%	15.5%
5000	11.2%	16.4%	19.5%	22.5%	21.8%	22.0%	21.7%
10,000	11.9%	17.5%	20.8%	24.6%	23.0%	23.2%	22.9%

The last part of the one-covariate simulation study is devoted to extreme quantile estimation. We take again s = 0 for the Cauchy distribution and use $10,000$ replications. For $n=500,m=1000,$ and $k=125,$ we estimate the extreme quantile x_p for $p=1/n=\mathrm{0.002.}$ Table 3 shows the variance reductions when using the adapted extreme quantile estimator in (10) relative to the standard estimator in (9), for various values of γ and g. The variance reductions are substantial, also for $\gamma <0$, and range from 17% to $33\%.$

Table 3 Variance reduction for the extreme quantile estimator.

g	γ
	–0.3	–0.2	–0.1	0	0.1	0.2	0.3
−0.25	24.07%	25.12%	23.6%	22.52%	20.57%	19.10%	16.56%
−0.125	27.92%	30.07%	29.26%	28.03%	26.74%	25.54%	22.98%
0	27.57%	31.34%	31.89%	31.28%	30.78%	30.02%	27.81%
0.125	23.99%	29.85%	31.95%	32.44%	32.60%	32.42%	30.75%
0.25	20.57%	27.74%	30.95%	32.34%	33.09%	33.46%	32.42%

We also study the effect of using the SSEs for both the extreme value index and the scale leading to the extreme quantile estimator in (10), compared to using a hybrid approach that only uses the SSE ${\widehat{\gamma }}_g$ for the extreme value index but keeps the classical $\widehat{\sigma }$ when estimating the scale.

We take g = 0 and estimate the extreme quantile x_p for $p=1/n$. Figure 3 shows the remarkable effect on the variance reduction when using the adapted estimator in (10) instead of the hybrid approach. The variance reduction is doubled in most of the cases, meaning that replacing $\widehat{\sigma }$ by the SSE ${\widehat{\sigma }}_g$ leads to a considerably improved performance.

Fig. 3 Relative variance reduction for the extreme quantile estimators.

To conclude the bivariate setting we briefly compare the here obtained sizable variance reductions with the reductions obtained from the asymptotic theory. It turns out that for the present sample sizes, when estimating γ, the asymptotic theory yields even higher reductions than the ones here obtained, partly due to the variability in the necessary estimation of the tail copulas, which is not reflected in the asymptotic variance. Based on limited comparisons it seems that for positive γ, asymptotic theory and simulations match somewhat better than in case γ is negative. Theoretically, for $\gamma \ge 0$ the asymptotic variance reduction of the extreme quantile estimator ${\widehat{x}}_{p_g}$ is the same as that of ${\widehat{\gamma }}_g$. In contrast, Table 3, in conjunction with Table 2, indicates that for the practically relevant extreme quantile estimation, the variance reductions obtained in simulations can be even higher than those inferred from the asymptotic theory.

At the end of this simulation section we now briefly consider the two-covariates setting. We begin with simulating data from the trivariate Cauchy distribution restricted to the first octant. The density is proportional to ${(1+x{S}^{-1}{x}^T)}^{-2}$ where S is now a 3 × 3 scale matrix with 1 on the diagonal and s off-diagonal, but $S_{23}=S_{32}=r$. We take respectively the values $s=r=0,s=r=0.8$ and $s=0.8,r=0.3$. When $s=r=0.8$, we have rather strong tail dependence, with R(1)-values of 0.77, whereas for $s=r=0$ the R(1)-values are 0.59. The case $s=0.8,r=0.3$ yields R(1)-values of 0.81, except the one between the covariates, which is 0.63. Then similarly as in (16) in the bivariate case, we transform the first coordinate to obtain our final data $(X_i,Y_{i,2},Y_{i,3})$; we will use $\gamma =-0.2,0$ and 0.2. We take $n=500,m=1000,k=125$ and $10,000$ replications. Table 4 shows that the variance reduction can be as high as 50%. We observe that the case $s=r=0.8$ indeed yields larger improvements than the case $s=r=0$ and the one-covariate setting in Table 1. Interestingly the case $s=0.8,r=0.3$ yields even larger variance reductions than the case $s=r=0.8$. The smaller tail dependence between the covariates, makes them as a pair more informative and leads to the much improved performance of ${\widehat{\gamma }}_g$.

Table 4 Variance reduction with two covariates.

	$s=0,r=0$			$s=0.8,r=0.8$			$s=0.8,r=0.3$
g	γ			γ			γ
	–0.2	0	0.2	–0.2	0	0.2	–0.2	0	0.2
−0.25	19.8%	11.9%	0.9%	37.4%	28.4%	8.0%	48.9%	42.1%	21.0%
−0.125	20.6%	19.5%	16.2%	37.8%	37.2%	29.3%	48.4%	49.3%	41.9%
0	19.2%	21.8%	21.3%	35.1%	39.1%	36.3%	44.5%	49.8%	48.2%
0.125	16.6%	22.1%	23.0%	30.6%	38.7%	38.5%	38.6%	48.2%	49.8%
0.25	13.1%	21.5%	23.5%	24.8%	36.9%	38.9%	31.0%	45.1%	49.5%

5 Application

In this section, we demonstrate an application using the SSE for analyzing forecasted precipitation data.

The national French weather service, Météo France, produces daily forecasted precipitation (in mm) at very high resolution (${0.1}^{°}\times {0.1}^{°}$) covering the mainland of France, between 2012 and 2017. Although usually an ensemble is forecasted, we use one single forecast (an ensemble member, not a forecasted mean or median) at every grid point and every day. To improve the forecasting model, meteorologists want to check if the forecasted precipitation shares the same distribution as the observed precipitation at the same location, particularly in the right tail. Consequently, the goal of this study is to estimate quantities such as the extreme value index and extreme quantiles of the forecasted precipitation distribution. We focus on forecasting at grid points that are close to an actual weather station.

Besides the forecasted precipitation, Météo France records at 123 weather stations the actual daily precipitation, between 1980 and 2017. We pair each weather station with a forecasting grid point that is closest to the station, and regard the two as the same location. When focusing on the fall seasons (91 days per year), at the 123 locations, we have 38 years actual precipitation data (3458 observations) with the last 6 years paired with forecasting data (546 observations). For the last 6 years, the paired data are dependent since the forecasting data are made to forecast the precipitation at the same location on the same day. Part of this dataset has been employed in a study comparing the spatial dependence structure of extreme forecasted precipitation and extreme observed precipitation in southern France; see Oesting and Naveau (2020).¹

Since it is challenging to conduct extreme value analysis, such as extreme quantile estimation, for the forecasted precipitation based on only 546 observations, we make use of the available information in the actual precipitation to improve the estimation accuracy, exploiting the semi-supervised setting. Observe that the SSE uses the tail dependence between the two datasets, but not the (marginal) tail heaviness of the observed precipitation. Therefore, the SSE of the extreme value index of the forecasted precipitation does not incorporate properties of the distribution of the observed precipitation. This justifies comparing this SSE with the MLE of the extreme value index of the observed precipitation.

To validate that our proposed methodology can be applied to the dataset, we perform two pre-tests on the actual precipitation data (3458 observations) at each station. First, we test whether the actual precipitation at each station possesses the same distribution across time, using the test statistic T₂ proposed in Einmahl, de Haan, and Zhou (2016) (with k = 200). Second, we test whether the extreme precipitation at each station can be regarded as independent over time, based on testing whether the extremal index is significantly different from 1, using the sliding block estimator proposed in Berghaus and Bücher (2018) (with b = 80). To achieve the 5% significance level for the joint test, we exclude all stations for which any of the two tests rejects the null at the 2.5% significance level. Eventually, for 100 stations both null hypotheses are not rejected, and we apply our proposed method to these 100 stations.

We use the SSE ${\widehat{\gamma }}_g$ with g = 0 to estimate the extreme value index, and compare it with the pseudo-MLE $\widehat{\gamma }$. For both estimators, we take k = 136. In particular, we estimate both the tail dependence coefficient R(1, 1) and the variance reduction factor $(1-{\nu }^2)R_g^2$ to demonstrate the presence of tail dependence and to evaluate the improvement when using the SSE. In addition, we estimate the “once per 10 year” extreme rainfall, that is the quantile at the probability level $1-1/910$, by (9) and (10), to compare the impact of using the SSE on practically relevant quantities. Last but not least, we estimate the extreme value index and the same level high quantile based on the actual precipitation, using the pseudo-MLE. In the estimation, we use the top $k_1=600$ observations, which in line with (3) corresponds to a lower fraction $k_1/(n+m)$ of POT compared to k/n used for the forecasted precipitation.

Table 5 shows the results for three selected stations. We select the three stations from very distant areas: one from the south, one from the northwest, and one from the southwest. For the station from the south, Nîmes, the estimated extreme value index is positive indicating a heavy-tailed distribution. Taking ${\nu }^2=n/(n+m)$, the reduction in variance is estimated at 16%. The difference of the two estimates of the extreme value index leads to a substantial difference in the quantile estimates: the quantile estimated using the SSE exceeds the usual quantile estimate with roughly 50%. In contrast, for the station in the northwest, Dieppe, the estimated extreme value index is around zero. The difference between the two point estimates is small, with the SSE having 17.5% variance reduction. The two quantile estimates are about the same. Finally, for the station in the southwest, Ciboure, both estimators lead to negative estimates, although not significantly different from zero. The variance reduction is at a pronounced level: 28.4%. The quantile estimate using the SSE is somewhat lower than the usual one. Observe that the quantile estimates based on the forecasting data are overall lower than those based on the observed precipitation, also when using the SSE. This suggests that the forecasting model produces forecasts with a lighter tail than that of the actual observations. We added 90% confidence intervals for the extreme value indices and extreme quantiles. Here for the SSE approach, the “common” quantities in the asymptotic variance that have to be estimated consistently (in particular ${(1+\gamma )}^2$ and $a(n/k)q_{\gamma }(k/(\mathit{\text{np}}))$) are estimated with the minimum of the estimates obtained from the MLE and SSE procedures.

Table 5 Estimation results for three stations.

	Nîmes	Dieppe	Ciboure
(LAT, LON)	(43.86, 4.41)	(49.93, 1.09)	(43.39,–1.69)
$\widehat{\gamma }$ (MLE)	0.358	–0.017	–0.056
	(0.166,0.549)	(–0.156,0.121)	(–0.189,0.077)
${\widehat{\gamma }}_0$ (SSE)	0.517	0.005	–0.088
	(0.341,0.692)	(–0.121,0.130)	(–0.201,0.025)
$\widehat{R}(1,1)$	0.478	0.456	0.537
Reduction	16.0%	17.5%	28.4%
Quantile (MLE)	114.85	33.05	59.32
	(37.77,191.93)	(22.37,43.74)	(41.72,76.93)
Quantile (SSE)	167.82	32.48	56.36
	(110.37,225.26)	(23.35,41.61)	(43.64,69.07)
${\widehat{\gamma }}_{\mathit{\text{obs}}}$ (MLE)	0.421	0.052	0.003
	(0.326, 0.516)	(–0.019,0.122)	(–0.064,0.071)
Quantile (Obs)	175.29	42.13	81.08
	(119.39, 231.20)	(35.31,48.95)	(68.77, 93.38)

Finally we present a heatmap of the estimated variance reduction factors across all 100 locations in Figure 4. The variance reduction using the SSE compared to the pseudo-MLE ranges from 9% to 34.8%, and is on average 18.5%.

Fig. 4 Heatmap of the variance reduction factor across 100 stations.

6 Proofs

We first present proofs for the one-covariate (bivariate) case and then extend the proofs to the multivariate case. The asymptotic normality for $\sqrt{k}(\widehat{\gamma }-\gamma )$, the first component of the pair in Proposition 2.1, is established in Drees, Ferreira, and de Haan (2004), see also the supplementary material, Section A. However, we cannot directly use that proof, since we have to keep track of the joint behavior of $\widehat{\gamma }$ and $\widehat{g}$. Nevertheless, we mimic that proof for both $\widehat{\gamma }$ and $\widehat{g}$, with observing that $\widehat{g}$ is based on dependent observations. In this respect the proof has to be adapted substantially. We begin with various lemmas which are needed for the main proofs.

Let C be a copula corresponding to the distribution function of $(-X_1,-Y_1).$ Let $(V_{1,1},V_{1,2}),$ $(V_{2,1},V_{2,2})\dots ,(V_{n,1},V_{n,2})$ be a random sample of size n from C and $V_{n+1,2},\dots ,V_{n+m,2}$ be a random sample of size m from the uniform-(0, 1) distribution, independent of the random sample from C. Clearly all the $V_{i,j},i=1,\dots ,n,j=1,2$, have also a uniform-(0, 1) distribution. Write $X_i=F_1^{-1}(1-V_{i,1}),i=1,\dots ,n,$ and $Y_l=F_2^{-1}(1-V_{l,2}),l=1,\dots ,n+m$. Then (X_i, Y _i), $i=1,\dots ,n$, and $Y_{n+1},\dots ,Y_{n+m}$ have the distributions as specified in the beginning of Section 2.

Consider the following uniform empirical distribution functions: $$\begin{array}{c}{\Gamma }_{n,j}(s)=\frac{1}{n}\sum _{i=1}^n{1}_{[0,s]}(V_{i,j}),0\le s\le 1,j=1,2,\\ {\Gamma }_{n+m}(t)=\frac{1}{n+m}\sum _{l=1}^{n+m}{1}_{[0,t]}(V_{l,2}),0\le t\le 1.\end{array}$$

The corresponding uniform tail empirical processes are $$\begin{array}{c}{w}_{n,j}(s)=\sqrt{k}\left[\frac{n}{k}{\Gamma }_{n,j}\left(\frac{k}{n}s\right)-s\right],\hspace{1em}0\le s\le 1,j=1,2,\\ w_{n+m}(t)=\sqrt{\frac{(n+m)k}{n}}\left[\frac{n}{k}{\Gamma }_{n+m}\left(\frac{k}{n}t\right)-t\right],\hspace{1em}0\le t\le 1.\end{array}$$

Define the Gaussian vector of processes $(W_1,W_2,W_3)$, where $W_j,j=1,2,3,$ is a standard Wiener process on $[0,T],T>0$, with covariances: $$\begin{array}{c}\text{cov}(W_1(s),W_2(t))=R(s,t),\hspace{1em}\text{cov}(W_1(s),W_3(t))=\mathit{\nu R}(s,t),\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{cov}(W_2(s),W_3(t))=\nu (s\wedge t),\hspace{1em}0\le s,t\le T.\end{array}$$

Let I denote the identity function. Then we have on ${(D[0,T])}^3$, for all $0\le \delta <\frac{1}{2}$ (17) $$\left(\frac{w_{n,1}}{I^{\delta }},\frac{w_{n,2}}{I^{\delta }},\frac{w_{n+m}}{I^{\delta }}\right)\stackrel{d}{\to }\left(\frac{W_1}{I^{\delta }},\frac{W_2}{I^{\delta }},\frac{W_3}{I^{\delta }}\right),\text{ as n}\to \infty .$$(17)

The proof of (17) is given in Ahmed and Einmahl (2019); note that in there T = 1, but the proof for arbitrary T > 0 follows similarly. Now using a Skorohod construction we obtain from (17) that (18) $$\begin{array}{c}\underset{0<s\le T}{\mathrm{\text{sup}}}\frac{|w_{n,j}(s)-W_j(s)|}{s^{\delta }}\stackrel{a.s.}{\to }0,j=1,2,\hspace{1em}\text{and }\\ \underset{0<s\le T}{\mathrm{\text{sup}}}\frac{|w_{n+m}(s)-W_3(s)|}{s^{\delta }}\stackrel{a.s.}{\to }0.\end{array}$$(18)

The processes in (18) are different from those in (17) but we keep the same notation, since the new vector $(w_{n,1},w_{n,2},w_{n+m})$ has the same distribution as the old vector and also the new vector $(W_1,W_2,W_3)$ has the same distribution as the old vector. In the sequel the X_i and Y_i are transformations as above of the uniform-(0,1) random variables on which the $w_{n,j}$ are based. We continue with the processes satisfying (18).

For convenience we introduce the following notation. Let f_n, h_n be positive functions on $[l_n,u_n]$. Then we write, as $n\to \infty$, $$f_n\stackrel{\mathbb{P}}{\asymp }{h}_n{|}_{l_n}^{u_n},$$if both $f_n(s)=O_{\mathbb{P}}(h_n(s))$ and $h_n(s)=O_{\mathbb{P}}(f_n(s))$ hold uniformly for $s\in [l_n,u_n]$. This notation is useful for the following lemma, see Shorack and Wellner (2009, p. 419).

Lemma 6.1.

Let ${\Gamma }_{n,j}^{-1},j=1,2,$ be the empirical quantile functions corresponding to ${\Gamma }_{n,j},j=1,2$, respectively. Then, as $n\to \infty$, ${\Gamma }_{n,j}\stackrel{\mathbb{P}}{\asymp }I{|}_{{\Gamma }_{n,j}^{-1}(1/(2n))}^1\hspace{1em}\text{and}\hspace{1em}{\Gamma }_{n,j}^{-1}\stackrel{\mathbb{P}}{\asymp }I{|}_{1/(2n)}^1,\hspace{1em}j=1,2.$

The following lemma states the weighted convergence of the tail quantile processes corresponding to ${\Gamma }_{n,j}^{-1},j=1,2,$ to the processes $-W_1$ and $-W_2$ in (18).

Lemma 6.2.

Let ${\Gamma }_{n,j}^{-1}$, be the empirical quantile functions corresponding to $w_{n,j}$, j = 1, 2, in (18) and let W_j be as in (18), j = 1, 2. Then for any $\delta <\frac{1}{2},$ as $n\to \infty ,$ (19) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}\frac{|\sqrt{k}(\frac{n}{k}{\Gamma }_{n,j}^{-1}\left(\frac{k}{n}s\right)-s)+W_j(s)|}{s^{\delta }}\stackrel{\mathbb{P}}{\to }0.$$(19)

Proof.

Write $v_{n,j}(s)=\sqrt{k}(\frac{n}{k}{\Gamma }_{n,j}^{-1}(\frac{k}{n}s)-s),j=1,2.$ Theorem 2.3 in Einmahl (1992) yields, as $n\to \infty ,$ (20) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}\frac{|v_{n,j}(s)-W_{n,j}(s)|}{s^{\delta }}\stackrel{\mathbb{P}}{\to }0,$$(20) where $W_{n,j}$ is an appropriate sequence of standard Wiener processes. Let W be a standard Wiener process and let $\epsilon >0.$ It is well-known that there exist an $\eta >0$, such that (21) $$\mathbb{P}\left(\underset{0<s\le \eta }{\mathrm{\text{sup}}}\frac{\left|W(s)\right|}{s^{\delta }}\ge \frac{\epsilon }{2}\right)\le \frac{\epsilon }{2}.$$(21)

Combining (20) and (21) yields, for large n, (22) $$\mathbb{P}\left(\underset{\frac{1}{2k}\le s\le \eta }{\mathrm{\text{sup}}}\frac{|v_{n,j}(s)|}{s^{\delta }}\ge \epsilon \right)\le \epsilon .$$(22)

Combining (18, 21, 22), and Lemma 1 in Vervaat (1972), yields (19). □

The next lemma is very similar to Lemma 3.1 in Drees, Ferreira, and de Haan (2004), but the lemma therein cannot be used here because we need specifically the approximation with the present W₁ in order to obtain the joint behavior of $\widehat{\gamma }$ and $\widehat{g}$.

Lemma 6.3.

Let $\epsilon >0$. Assume that (3) and (4) hold and $\sqrt{k}A\left(\frac{n}{k}\right)=O(1)$, as $n\to \infty$. Then for suitably chosen functions a and A in (4), as $n\to \infty$, $$\begin{array}{c}\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{s}^{\gamma +1/2+\epsilon }|\sqrt{k}\left(\frac{X_{n-\left[\mathit{\text{ks}}\right]:n}-U_1\left(\frac{n}{k}\right)}{a\left(\frac{n}{k}\right)}-\frac{s^{-\gamma }-1}{\gamma }\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-s^{-\gamma -1}{W}_1(s)-\sqrt{k}A\left(\frac{n}{k}\right)\Psi (s^{-1})|\stackrel{\mathbb{P}}{\to }0.\end{array}$$

Proof.

From (4) we obtain inequality (2.3.17) in de Haan and Ferreira (2006): for any $\theta ,\delta >0$ to be specified later, there exists $t_0=t_0(\theta ,\delta )$ such that for all $t,\mathit{\text{tx}}\ge t_0$, $$\left|\frac{\frac{U_1(\mathit{\text{tx}})-U_1(t)}{a(t)}-\frac{x^{\gamma }-1}{\gamma }}{A(t)}-\Psi (x)\right|\le \theta x^{\gamma +\rho }\mathrm{\text{max}}(x^{\delta },x^{-\delta }).$$

We replace tx by $1/{\Gamma }_{n,1}^{-1}\left(\frac{k}{n}s\right)$ and t by $\frac{n}{k}$. Then we have, writing $\check{s}=\frac{n}{k}{\Gamma }_{n,1}^{-1}\left(\frac{k}{n}s\right)$, with probability tending to 1, as $n\to \infty$, (23) $$\begin{array}{c}\left|\frac{X_{n-\left[\mathit{\text{ks}}\right]:n}-U_1\left(\frac{n}{k}\right)}{a\left(\frac{n}{k}\right)}-\frac{{\check{s}}^{-\gamma }-1}{\gamma }-A\left(\frac{n}{k}\right)\Psi \left({\check{s}}^{-1}\right)\right|\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\le \left|A\left(\frac{n}{k}\right)\right|\theta {\check{s}}^{-\gamma -\rho }·\mathrm{\text{max}}({\check{s}}^{-\delta },{\check{s}}^{\delta }).\end{array}$$(23)

Define $f(s)=\frac{s^{-\gamma }-1}{\gamma }$. Then by a Taylor expansion for some ${\check{\Theta }}_n(s)$ between $\check{s}$ and s we have $f(\check{s})-f(s)=f\prime (s)(\check{s}-s)+\frac{1}{2}f″({\check{\Theta }}_n(s)){(\check{s}-s)}^2.$ Lemma 6.1 implies ${\check{\Theta }}_n\stackrel{\mathbb{P}}{\asymp }I{|}_{\frac{1}{2k}}^1$ and thus $f″({\check{\Theta }}_n)\stackrel{\mathbb{P}}{\asymp }{I}^{-\gamma -2}{|}_{\frac{1}{2k}}^1$. Next, by Lemma 6.2 and the fact that for all ${\delta }_1<\frac{1}{2},\underset{0\le s\le 1}{\mathrm{\text{sup}}}|W_1(s)|/s^{{\delta }_1}=O_{\mathbb{P}}(1),$ we have that $\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{(\check{s}-s)}^2/s^{2{\delta }_1}=O_{\mathbb{P}}(1/k),$ as $n\to \infty$. This and again Lemma 6.2 with $\delta ={\delta }_1$ yield, as $n\to \infty$, uniformly for all $\frac{1}{2k}\le s\le 1$, $$\begin{array}{c}f\left(\check{s}\right)-f(s)=-s^{-\gamma -1}\frac{1}{\sqrt{k}}\left(-W_1(s)+s^{{\delta }_1}{o}_{\mathbb{P}}(1)\right)\\ \hspace{1em}+s^{-\gamma -2+2{\delta }_1}{O}_{\mathbb{P}}\left(\frac{1}{k}\right).\end{array}$$

Choose δ₁ such that $\frac{1-\epsilon }{2}<{\delta }_1<\frac{1}{2}$. Then ${\delta }_1>\frac{1}{2}-\epsilon$ and $2{\delta }_1+\epsilon >1$. Hence, as $n\to \infty$, $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{s}^{-\frac{3}{2}+\epsilon +2{\delta }_1}\le \mathrm{\text{max}}\left(1,{(2k)}^{\frac{3}{2}-\epsilon -2{\delta }_1}\right)=o\left(\sqrt{k}\right).$$

Therefore, as $n\to \infty$, uniformly for all $\frac{1}{2k}\le s\le 1$, (24) $$\begin{array}{c}f(\check{s})=f(s)+\frac{1}{\sqrt{k}}{s}^{-\gamma -1}\\ \hspace{1em}\left(W_1(s)+s^{{\delta }_1}{o}_{\mathbb{P}}(1)+s^{-1+2{\delta }_1}{O}_{\mathbb{P}}\left(\frac{1}{\sqrt{k}}\right)\right)\\ =f(s)+\frac{1}{\sqrt{k}}{s}^{-\gamma -1}(W_1(s)+s^{1/2-\epsilon }(s^{{\delta }_1-1/2+\epsilon }{o}_{\mathbb{P}}(1)\\ \hspace{1em}+s^{-3/2+\epsilon +2{\delta }_1}{O}_{\mathbb{P}}\left(\frac{1}{\sqrt{k}}\right)))\\ =\frac{s^{-\gamma }-1}{\gamma }+\frac{1}{\sqrt{k}}{s}^{-\gamma -1}\left(W_1(s)+s^{1/2-\epsilon }{o}_{\mathbb{P}}(1)\right).\end{array}$$(24)

From the mean value theorem, for some ${\Theta }_n(s)$ between $\check{s}$ and s $$\Psi \left({\check{s}}^{-1}\right)=\Psi (s^{-1})-\Psi \prime (1/{\Theta }_n(s)){({\Theta }_n(s))}^{-2}\left(\check{s}-s\right).$$

As above, ${\Theta }_n\stackrel{\mathbb{P}}{\asymp }I{|}_{\frac{1}{2k}}^1$, which implies that as $n\to \infty$, uniformly for $\frac{1}{2k}\le s\le 1$, $$\left|\Psi \prime (1/{\Theta }_n(s)){({\Theta }_n(s))}^{-2}\right|=s^{-\gamma -\rho -1}(1+|\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s|)O_{\mathbb{P}}(1).$$

Hence, using Lemma 6.2 with $\delta ={\delta }_1$ (as above), we have uniformly for $\frac{1}{2k}\le s\le 1$, $$\begin{array}{c}A\left(\frac{n}{k}\right)\left(\Psi ({\check{s}}^{-1})-\Psi (s^{-1})\right)\\ \hspace{1em}=\frac{1}{\sqrt{k}}A\left(\frac{n}{k}\right)s^{-\gamma -\rho -1+{\delta }_1}(1+|\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s|)O_{\mathbb{P}}(1).\end{array}$$

With δ₁ chosen as above, we have that as $n\to \infty$, uniformly for $\frac{1}{2k}\le s\le 1$, (25) $$A\left(\frac{n}{k}\right)\Psi ({\check{s}}^{-1})=A\left(\frac{n}{k}\right)\Psi (s^{-1})+\frac{1}{\sqrt{k}}{s}^{-\gamma -\epsilon -\frac{1}{2}}{o}_{\mathbb{P}}(1).$$(25)

Next consider the right-hand side of (23), where we take $\delta <1/2$. Using Lemma 6.1, it can be bounded, uniformly for $\frac{1}{2k}\le s\le 1$, by (26) $$\begin{array}{c}\theta \left|A\left(\frac{n}{k}\right)\right|s^{-\gamma -\rho -\delta }{O}_{\mathbb{P}}(1)=\theta \sqrt{k}\left|A\left(\frac{n}{k}\right)\right|\frac{1}{\sqrt{k}}{s}^{-\gamma -\delta }{O}_{\mathbb{P}}(1)\\ \hspace{1em}=\theta \frac{1}{\sqrt{k}}{s}^{-\gamma -\epsilon -1/2}{s}^{\epsilon +1/2-\delta }{O}_{\mathbb{P}}(1)=\theta \frac{1}{\sqrt{k}}{s}^{-\gamma -\epsilon -1/2}{O}_{\mathbb{P}}(1).\end{array}$$(26)

Now, plugging (24)–(26) into inequality (23) and noting that $\theta >0$ can be chosen arbitrarily small, we obtain the statement in the lemma. □

Define $$Z_n(s)=\sqrt{k}\left(\frac{X_{n-\left[\mathit{\text{ks}}\right]:n}-X_{n-k:n}}{a\left(\frac{n}{k}\right)}-\frac{s^{-\gamma }-1}{\gamma }\right).$$

Then for functions a and A as in Lemma 6.3, for any $\epsilon >0$, uniformly for $\frac{1}{2k}\le s\le 1$, $$Z_n(s)=s^{-\gamma -1}{W}_1(s)-W_1(1)+\sqrt{k}A\left(\frac{n}{k}\right)\Psi (s^{-1})$$ $$\hspace{1em}+o_{\mathbb{P}}(1)s^{-\gamma -1/2-\epsilon }.$$

Hence, for $\gamma >-\frac{1}{2}$, (27) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{s}^{\gamma +1/2+\epsilon }|Z_n(s)|=O_{\mathbb{P}}(1).$$(27)

Proposition 6.1.

Under the conditions of Lemma 6.3, for $\gamma >-\frac{1}{2}$ and $\gamma \ne 0$, with probability tending to 1, there exists a unique maximizer of the likelihood function based on ${\left\{X_i\right\}}_{i=1}^n$ denoted as $\widehat{\gamma }$, such that as $n\to \infty$, $$\sqrt{k}(\widehat{\gamma }-\gamma )-\frac{{(\gamma +1)}^2}{\gamma }{\int }_0^1(s^{\gamma }-(2\gamma +1)s^{2\gamma })Z_n(s)\mathit{\text{ds}}=o_{\mathbb{P}}(1),$$and, for γ = 0, $\sqrt{k}\widehat{\gamma }+{\int }_0^1(2+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s)Z_n(s)\mathit{\text{ds}}=o_{\mathbb{P}}(1).$

Proof.

The existence of $\widehat{\gamma }$ follows from Theorem 4.1 in Zhou (2009). Then, using Lemma 6.3 and (27) above in conjunction with Lemma 3.2 in Drees, Ferreira, and de Haan (2004), the result is obtained following the same steps as in the proof of Proposition 3.1 in Drees, Ferreira, and de Haan (2004). A sketch of the main steps of the proof can be found in the supplementary material, Section A. □

To study the asymptotic behavior of $\widehat{g}$ we need the following result. Define $$\begin{array}{c}{\tilde{w}}_n(s)=\frac{n}{\sqrt{k}}\left({\Gamma }_{n+m}\left({\Gamma }_{n,2}^{-1}\left(\frac{k}{n}s\right)\right)-\frac{k}{n}s\right)\hspace{1em}\text{ and }\\ \tilde{W}(s)=\nu W_3(s)-W_2(s).\end{array}$$

Lemma 6.4.

Assume that F₂ is continuous and k satisfies (3), then for any $0\le \delta <\frac{1}{2}$, as $n\to \infty$, $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}\frac{|{\tilde{w}}_n(s)-\tilde{W}(s)|}{s^{\delta }}\stackrel{\mathbb{P}}{\to }0.$$

Proof.

We have $${\tilde{w}}_n(s)=\sqrt{\frac{n}{n+m}}{w}_{n+m}\left(\mathrm{}\frac{n}{k}{\Gamma }_{n,2}^{-1}\left(\mathrm{}\frac{k}{n}s\mathrm{}\right)\mathrm{}\mathrm{}\right)+\frac{n}{\sqrt{k}}\left(\mathrm{}{\Gamma }_{n,2}^{-1}\left(\mathrm{}\frac{k}{n}s\mathrm{}\right)-\frac{k}{n}s\mathrm{}\mathrm{}\right).$$

Define $\widehat{s}=\frac{n}{k}{\Gamma }_{n,2}^{-1}\left(\frac{k}{n}s\right)$. From Lemma 6.2 with j = 2, (3) and (21), we see that it suffices to show that, as $n\to \infty$, (28) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}\frac{|w_{n+m}(\widehat{s})-W_3(s)|}{s^{\delta }}\stackrel{\mathbb{P}}{\to }0.$$(28)

Let $s_0\in (0,1)$. We first handle the region $s\ge s_0$. Obviously we have $1/s^{\delta }\le 1/s_0^{\delta }.$ By Lemma 6.2, as $n\to \infty$, (29) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}\left|\widehat{s}-s\right|\stackrel{\mathbb{P}}{\to }0.$$(29)

Using this, (18), and the uniform continuity of W₃ we obtain, as $n\to \infty$, $$\underset{s_0\le s\le 1}{\mathrm{\text{sup}}}\frac{|w_{n+m}(\widehat{s})-W_3(s)|}{s^{\delta }}\stackrel{\mathbb{P}}{\to }0.$$

It remains to show that for $\epsilon >0$ there exists $s_0\in (0,1)$ such that for large n $$\mathbb{P}\left(\underset{\frac{1}{2k}\le s\le s_0}{\mathrm{\text{sup}}}\frac{|w_{n+m}(\widehat{s})-W_3(s)|}{s^{\delta }}\ge 3\epsilon \right)\le 3\epsilon .$$

Using again (21), for this it suffices to show that $\mathbb{P}(\underset{\frac{1}{2k}\le s\le s_0}{\mathrm{\text{sup}}}|w_{n+m}(\widehat{s})|/s^{\delta }\ge 2\epsilon )\le 2\epsilon .$ Using Lemma 6.1, the proof is complete if we show that for all $\epsilon >0,\kappa >0$ there exists $s_0\in (0,1)$ such that for large n $\mathbb{P}(\underset{\frac{1}{2k}\le s\le s_0}{\mathrm{\text{sup}}}|w_{n+m}(\widehat{s})|/{\widehat{s}}^{\delta }\ge 2\kappa )\le \epsilon .$ We have $$\begin{array}{c}\mathbb{P}\left(\underset{\frac{1}{2k}\le s\le s_0}{\mathrm{\text{sup}}}\frac{|w_{n+m}(\widehat{s})|}{{\widehat{s}}^{\delta }}\ge 2\kappa \right)\\ \hspace{1em}\le \mathbb{P}\left(\underset{0<t\le 2s_0}{\mathrm{\text{sup}}}\left|\frac{w_{n+m}(t)}{t^{\delta }}\right|\ge \kappa \right)+\mathbb{P}(\widehat{s}>2s_0).\end{array}$$

From (18) and (21), we have that for small enough $s_0\in (0,1)$ the first term on the right is bounded by $\epsilon /2$ for large n, and using (29) we obtain that the second term on the right also does not exceed $\epsilon /2$ for large n. □

In the following we prove a result for the tail quantile process based on ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$ instead of ${\left\{X_i\right\}}_{i=1}^n$. The proof of the next lemma uses Lemma 6.4 and is very similar to but easier than that of Lemma 6.3, and hence will be omitted.

Lemma 6.5.

Let $\epsilon >0$. Assume that F₂ is continuous and that (3) holds, then, as $n\to \infty$, $$\begin{array}{c}\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{s}^{g+1/2+\epsilon }|\sqrt{k}\left(\frac{\left({\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-\frac{{(\frac{n}{k})}^g-1}{g}\right)}{{\left(\frac{n}{k}\right)}^g}-\frac{s^{-g}-1}{g}\right)\\ +t{s}^{-g-1}\tilde{W}(s)|\stackrel{\mathbb{P}}{\to }0.\end{array}$$

Define $$H_n(s):=\sqrt{k}\left(\frac{{\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}}{{(\frac{n}{k})}^g}-\frac{s^{-g}-1}{g}\right).$$

Then for any $\epsilon >0$, uniformly for $s\in [\frac{1}{2k},1]$, $$H_n(s)=\tilde{W}(1)-s^{-g-1}\tilde{W}(s)+o_{\mathbb{P}}(1)s^{-g-1/2-\epsilon }.$$

Hence, for $g>-\frac{1}{2}$, (30) $$\underset{\frac{1}{2k}\le s\le 1}{\mathrm{\text{sup}}}{s}^{g+1/2+\epsilon }|H_n(s)|=O_{\mathbb{P}}(1).$$(30)

Next we show a version of Lemma 3.2 in Drees, Ferreira, and de Haan (2004) based on ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$.

Lemma 6.6.

Assume that F₂ is continuous and k satisfies (3). Let g_n be a sequence of random variables such that $g_n=g+O_{\mathbb{P}}(k^{-1/2}).$ Then, if $-1/2<g<0$ or g > 0, as $n\to \infty ,$ (31) $$\mathbb{P}\left(\mathrm{}\mathrm{}1+g_n\frac{{\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}}{{(\frac{n}{k})}^g}\ge C_n{s}^{-g},\text{ forall }\mathrm{s}\in \left[\mathrm{}\mathrm{}\frac{1}{\text{2k}},1\mathrm{}\right]\mathrm{}\right)\to 1,$$(31) for some random variables $C_n>0$ such that $1/C_n=O_{\mathbb{P}}(1).$

If g = 0, as $n\to \infty$, (32) $$\mathbb{P}\left(\mathrm{}1+g_n\left({\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}\mathrm{}\right)\ge \frac{1}{2},\text{ forall }\mathrm{s}\in \left[\mathrm{}\frac{1}{\text{2k}},1\mathrm{}\right]\mathrm{}\right)\to 1,$$(32) and (33) $$\underset{s\in [0,1]}{\mathrm{\text{sup}}}{\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}=O_{\mathbb{P}}(\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}k).$$(33)

Proof.

Consider first $-1/2<g<0$ or g > 0. Applying Lemma 6.1 to ${\Gamma }_{n+m}$ and ${\Gamma }_{n,2}^{-1}$ yields, as $n\to \infty$, $${\Gamma }_{n+m}\left({\Gamma }_{n,2}^{-1}\left(\frac{k}{n}I\right)\right)\stackrel{\mathbb{P}}{\asymp }\frac{k}{n}I{|}_{\frac{1}{2k}}^1.$$

Define $G_n(s)={\Gamma }_{n+m}({\Gamma }_{n,2}^{-1}(\frac{k}{n}s))+\frac{1}{2(n+m)},s\in (0,1]$. Hence, as $n\to \infty$, (34) $$G_n\stackrel{\mathbb{P}}{\asymp }\frac{k}{n}I{|}_{\frac{1}{2k}}^1.$$(34)

Observe that for $g\ne 0$ $$\begin{array}{c}{s}^g\left(1+g_n\frac{{\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}}{{(\frac{n}{k})}^g}\right)\\ \hspace{1em}=s^g\left(1+\frac{g_n}{g}\frac{[{(G_n(s))}^{-g}-{(G_n(1))}^{-g}]}{{(\frac{n}{k})}^g}\right)\\ \hspace{1em}=\frac{g_n}{g}{\left(\frac{G_n(s)}{(\mathit{\text{ks}})/n}\right)}^{-g}+s^g[\left(1-{\left(\frac{G_n(1)}{k/n}\right)}^{-g}\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-\left(\frac{g_n}{g}-1\right){\left(\frac{G_n(1)}{k/n}\right)}^{-g}]=:T_1(s)+s^g[T_2-T_3].\end{array}$$

From (34) and $g_n/g\stackrel{\mathbb{P}}{\to }1$, we have that $1/\underset{s\in [1/(2k),1]}{\mathrm{\text{inf}}}{T}_1(s)=O_{\mathbb{P}}(1)$, as $n\to \infty$. Lemma 6.4 for s = 1 yields that $T_2=O_{\mathbb{P}}(1/\sqrt{k})$ and hence, since $g>-1/2,\underset{s\in [1/(2k),1]}{\mathrm{\text{sup}}}{s}^g·T_2\stackrel{\mathbb{P}}{\to }0$. By the assumption on g_n and again (34) we obtain similarly $\underset{s\in [1/(2k),1]}{\mathrm{\text{sup}}}{s}^g·T_3\stackrel{\mathbb{P}}{\to }0$. This yields (31).

In case g = 0, for $1/(2k)\le s\le 1$, (35) $${\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}=-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{G}_n(s)+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{G}_n(1)\le 2\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{A}_n-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s,$$(35) with $$A_n=\mathrm{\text{max}}\left(\underset{s\in [\frac{1}{2k},1]}{\mathrm{\text{sup}}}\frac{G_n(s)}{\frac{k}{n}s},\underset{s\in [\frac{1}{2k},1]}{\mathrm{\text{sup}}}\frac{\frac{k}{n}s}{G_n(s)}\right).$$

If $g_n\ge 0$, then $1+g_n({\tilde{Y}}_{n-[\mathit{\text{ks}}]:n}-{\tilde{Y}}_{n-k:n})\ge 1.$ If $g_n<0$, then for $1/(2k)\le s\le 1$, $$1+g_n\left({\tilde{Y}}_{n-\left[\mathit{\text{ks}}\right]:n}-{\tilde{Y}}_{n-k:n}\right)\ge 1+g_n(2\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{A}_n+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}2+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}k).$$

Since, as $n\to \infty ,A_n=O_{\mathbb{P}}(1)$ and $g_n=O_{\mathbb{P}}(k^{-1/2})$, we obtain (32). Finally, the $\mathrm{\text{sup}}$ in (33) is attained at $s=1/(2k)$. Hence, (35) yields (33). □

Finally, the following proposition provides the asymptotic behavior of the pseudo-MLE based on ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$.

Proposition 6.2.

Assume that F₂ is continuous and k satisfies (3). For $g>-\frac{1}{2}$ and $g\ne 0$, with probability tending to 1, there exists a unique maximizer of the likelihood function based on ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$, denoted as $\widehat{g}$, such that, as $n\to \infty$, $$\sqrt{k}(\widehat{g}-g)-\frac{{(g+1)}^2}{g}{\int }_0^1(s^g-(2g+1)s^{2g})H_n(s)\mathit{\text{ds}}=o_{\mathbb{P}}(1),$$and, for g = 0, $\sqrt{k}\widehat{g}+{\int }_0^1(2+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}s)H_n(s)\mathit{\text{ds}}=o_{\mathbb{P}}(1).$

Proof.

The proof follows the same steps as that of Proposition 6.1; see supplementary material, Section A. The main difference is that ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$ are not iid observations. Nevertheless, Lemma 6.5 guarantees that statistics based on the tail quantile process of ${\left\{{\tilde{Y}}_i\right\}}_{i=1}^n$, for example the Hill estimator for g > 0, possess similar asymptotic behavior as in the iid case, with the only difference that the random limit is driven by a proper functional of $\tilde{W}$ instead. Such asymptotic expansions are sufficient to ensure that the proof can still be realized.

Then, following the steps explained in the proof of Proposition 6.1, by using Lemma 6.5, (30), and Lemma 6.6, we get the analogous result as in Proposition 6.1. Note that Lemma 6.6 plays an analogous role as Lemma 3.2 in Drees, Ferreira, and de Haan (2004). □

Proof of Proposition 2.1.

Combining (18), Propositions 6.1 and 6.2 we obtain, as $n\to \infty ,$ $$\left(\sqrt{k}(\widehat{\gamma }-\gamma ),\sqrt{k}(\widehat{g}-g)\right)\stackrel{d}{\to }\left(\Omega ,\tilde{\Omega }\right),$$where $$\begin{array}{c}\Omega =\frac{{(\gamma +1)}^2}{\gamma }{\int }_0^1(s^{\gamma }-(2\gamma +1)s^{2\gamma })(s^{-\gamma -1}{W}_1(s)\\ -W_1(1))\mathit{\text{ds}}+\frac{\lambda (\gamma +1)}{(1-\rho )(1+\gamma -\rho )}\end{array}$$

and $$\tilde{\Omega }=\frac{{(g+1)}^2}{g}{\int }_0^1(t^g-(2g+1)t^{2g})(\tilde{W}(1)-t^{-g-1}\tilde{W}(t))\mathit{\text{dt}}.$$

Since the Wiener processes involved have mean zero, we obtain immediately the mean of the limiting pair. Also the individual variances of the limiting pair follow readily, see Drees, Ferreira, and de Haan (2004). The proof is completed by calculating the covariance of the two terms, which is deferred to the supplementary material, Section B. □

Proof of Theorem 2.1.

From the uniform consistency of $\widehat{R}$ on ${[0,1]}^2$, it can be shown that ${\widehat{R}}_g\stackrel{\mathbb{P}}{\to }{R}_g.$ Using the latter convergence in combination with Proposition 2.1 we obtain that, as $n\to \infty ,$ (36) $$\sqrt{k}({\widehat{\gamma }}_g-\gamma )=\sqrt{k}(\widehat{\gamma }-\gamma )-\frac{1+\gamma }{1+g}{R}_g\sqrt{k}(\widehat{g}-g)+o_{\mathbb{P}}(1).$$(36)

Now Proposition 2.1 in conjunction with the continuous mapping theorem yields (8). □

The proof of Proposition 3.1 can be given along the same lines as that of Proposition 2.1 and will be omitted. Note that the lemmas and propositions needed for the proof of Proposition 2.1 are of univariate nature and that hence immediately very similar lemmas can be stated (and proved) in the more-covariates case. Once these results are given, only a straightforward covariance calculation remains; see Ahmed and Einmahl (2019) for the joint weak convergence of all the tail empirical processes involved.

Proof of Theorem 3.1.

From the uniform consistency of the tail copula estimators we obtain ${\widehat{H}}_{1j}^{-1}\stackrel{\mathbb{P}}{\to }{H}_{1j}^{-1},j=1,\dots ,d$, which in combination with Proposition 3.1 yields that, as $n\to \infty ,$ $$\sqrt{k}({\widehat{\gamma }}_g-\gamma )=\sqrt{k}(\widehat{\gamma }-\gamma )+\frac{1+\gamma }{1+g}\sum _{j=2}^d\frac{H_{1j}^{-1}}{H_{11}^{-1}}\sqrt{k}({\widehat{g}}_j-g)+o_{\mathbb{P}}(1).$$

Now Proposition 3.1 and the continuous mapping theorem yield (15). □

Supplementary Materials

The supplementary material provides additional proofs needed, additional simulation results, as well as all codes used for the simulations.

Acknowledgments

We thank the Editor, Associate Editor, and three Referees for a careful reading of the article and various helpful comments that greatly improved the article. John Einmahl holds the Arie Kapteyn Chair 2019–2022 and gratefully acknowledges the corresponding research support.

Disclosure Statement

The authors report there are no competing interests to declare.

Notes

1 We thank Meteo France for providing the dataset. The data are available upon request, with the consent of Meteo France.