Keywords

Statistical decision , Asymptotic approximation , Samples , Differentiable models , Logistic modes

1 . Introduction

As regards statistical decision for bivariate extremes, we will use an approach similar to the one used for univariate extremes: we will replace, for large samples, \(\mathrm{ F^{n}( x,y ) }\) by an asymptotic approximation

\(\mathrm{ L ( \frac{x- \lambda }{ \delta },\frac{x- \lambda ’}{ \delta ’}) }\),

where \(\mathrm{ L(x-y) }\) will have the expressions given before, thus reducing the difficulty of not knowing \(\mathrm{ F^{n}( x,y ) }\) to the simpler \(\mathrm{ ( ! ) }\) difficulties of not knowing the dependence structure ( \(\mathrm{ k ( . ) }\) or \(\mathrm{ A ( . ) }\) for instance) and four independent margin parameters. Recall that we have the results similar to the ones used to justify the substitution of \(\mathrm{ F^{n} ( x ) }\)  by \(\mathrm{ L ( \frac{x- \lambda }{ \delta } ) }\) in the univariate case:

The stability and the uniformity of convergence are valid as well as results similar to those of Watson and others, showing that even in the case of some dependence between the successive observations the limiting distributions are still the same, are also true.

In the asymptotic approximation of   \(\mathrm{ L ( \frac{x- \lambda }{ \delta },\frac{x- \lambda ’}{ \delta ’}) }\), to the actual and unknown \(\mathrm{ F^{n}( x,y ) }\), we will use, in principle, a dependence function, although Geffroy’s sufficient condition could suggest the independence model. The existence of a strong correlation between the maxima (or minima), a correlation not completely destroyed in some applications by the relatively large number of observations from which the maxima (or minima) are taken, suggests the use of an asymptotic bivariate extremes distribution function as a better approximation to the actual (unknown) distribution function of bivariate extremes. This stand is also justified by the fact that the discrepancy between the independence and the dependence cases may be very large in a probabilistic sense but not expressed by the usual indicators. As an example we can cite the fact that the uniform distance between the biextremal model with standard Gumbel margins, where the distribution function is \(\mathrm{ \Lambda ( x,y \vert \theta) =exp \{ -e^{-x}+e^{-y}-min ( e^{-x}, \theta e^{-y} ) \} }\)(note that \(\mathrm{ \Lambda ( x,y \vert 0 ) = \Lambda ( x ) \cdot \Lambda ( y ) }\) — independence — and that \(\mathrm{ \Lambda ( x,y \vert 1 ) = \Lambda ( min⁡( x,y ) ) }\), — the diagonal case—) , and the independence case is

\( \mathrm {\begin{array}{c}\\ \mathrm{sup} \\ \mathrm{x,y} \end{array} | \Lambda \left( x,y \vert \theta \right) - \Lambda \left( x,y \vert 0 \right) \vert = \theta \left( 1+ \theta \right) ^{-1-1/ \theta } }\);

consequently if we take the independence model as an approximation to a biextremal model with \(\mathrm{ \theta }\) close to 1, i.e., close to the diagonal case, we could obtain probability evaluations for quadrants with an error around \(\mathrm{ 1 /4 }\), the value of the distance for \(\mathrm{ \theta=1 }\), but that can differ by almost 1 if the set used is a strip containing the first diagonal!

Experience will finally, decide the fitness and accuracy of those approximations. But, for the moment, as the \(\mathrm{ (m,N) }\) natural model shows, we can have any finite number of parameters, forgetting the four of the margins. And it also can be shown that for no one-parameter models (with standard margins), do there exist sufficient statistics. An alternative is the intrinsic estimation of the dependence function \(\mathrm{ (k\,or\,A) }\), as will be seen.

For some models we do not have the best test, the best estimation procedure, etc. (we may have one, but not necessarily the best) although in other situations nothing at all is known. An example: the statistical choice of bivariate extreme models — so important for applications — is now beginning to be considered! In general, the few papers up to now have chosen one model from the beginning or compare two of them by the use of the Kolmogoroff-Smirnov statistic or another test; see, for instance, Gumbel and Golstein (1964).

As a general rule in statistical decision for bivariate maxima with Gumbel margins, the papers Tiago de Oliveira (1970), (1975a), (1975b), (1980), (1984) and (1989) contain results that are improved here.

As regards modelling, we have to deal with some models which have a planar density and with others, which appear naturally, that have a singular component.

Finally, let us recall that as the margins in both cases (Gumbel and exponential) have all moments, we also have crossed moments, etc., a fact that will be used.

Note that for all models to be considered, we only give the correlation coefficients and index of dependence when they have a compact form. The correlation coefficients can then easily be used for statistical decision. The non-parametric correlation coefficients and the index of dependence are the same for the corresponding \(\mathrm{ \Lambda \left( x,y \vert \theta \right) }\) and \(\mathrm{ S \left( x’,y’ \vert \theta \right) }\); the classical correlation coefficient is different.

It is important to note that separation between models is difficult for moderate samples. A graphical explanation is the essence of Exercise 3.23.

2 . The differentiable models: statistical decision

Until now only two models of bivariate extremes with a planar density (with a point exception in one case) have been considered and studied in some depth in the literature: the logistic model and the mixed model, both appearing in Gumbel (1962) and (1966). Both are symmetrical but nothing precludes the existence of asymmetric models, which is the case in the majority of non-differentiable models to be considered in the next section. Here, as well as for non-differentiable models, we will always study them under the form of bivariate maxima with standard Gumbel margins, i.e., with the distribution function \(\mathrm{ \Lambda ( x,y ) = ( \Lambda ( x ) \Lambda ( y ) ) ^{k ( y-x ) } }\) with the dependence function \(\mathrm{ k ( w ) }\); at the end of each case we will consider the corresponding forms for bivariate minima with standard exponential margins, i.e., with the survival function \(\mathrm{ S ( x’,y’ ) =exp⁡ \{ - \left( x’+y’ \right) A ( \frac{y’}{x’+y’}) \} ( x’+y’ \geq 0 ) }\) obtained by the transformation \(\mathrm{ X’=e^{-X},Y’=e^{-Y} }\) and with \(\mathrm{ A ( u ) =k ( log ( \frac{1-u}{u} ) ) }\), as known.

The logistic model has the distribution function \(\mathrm{ \Lambda ( x,y \vert \theta ) =exp⁡ \{ -e^{-x/ ( 1- \theta ) }+e^{-y/ ( 1- \theta ) } ) ^{1- \theta } \} ,0 \leq \theta <1 }\) (for \(\mathrm{ \theta =0 }\) we have independence and for \(\mathrm{ \theta =1^- }\) we get the diagonal case, which is the only situation without density). It is so called because

\(\mathrm{ Prob \{ Y-X \leq w \vert \theta \} =D ( w \vert \theta ) = ( 1+e^{-w/ ( 1- \theta ) } ) ^{-1} }\)

corresponding to the dependence function

\(\mathrm{ k ( w \vert \theta ) =\frac{ ( 1+e^{-w/( 1- \theta ) } ) ^{-1- \theta }}{1+e^{-w}} }\);

note that we must have \(\mathrm{ 0 \leq \theta \leq 1 }\), as comes from the inequality \(\mathrm{ \frac{max⁡ ( 1+e^{w} ) }{1+e^{w}} \leq k ( w ) \leq 1 }\).

The margins are exchangeable as shown by the expression of \(\mathrm{ \Lambda \left( x,y \vert \theta \right) }\) , by the fact that \(\mathrm{ k ( -w \vert \theta ) =k ( w \vert \theta ) }\) or even because \(\mathrm{ D ( w \vert \theta ) +D( -w \vert \theta ) =1 }\).

The equidistribution curves are very easy to draw. The median curve has the equation

\(\mathrm{ e^{-x/ \left( 1- \theta \right) }+e^{-y/ \left( 1- \theta \right) }= \left( log~2 \right) ^{1/ \left( 1- \theta \right) } }\)

which is similar, with similarity ratio \(\mathrm{ 1/ \left( 1- \theta \right) }\), to the equidistribution curve for the independence case.

The classical correlation coefficient is \(\mathrm{ \rho \left( \theta \right) = \theta \left( 2- \theta \right) }\), increasing from \(\mathrm{ \rho \left( 0 \right) =0~to~ \rho \left( 1 \right) =1 }\), the difference-sign correlation coefficient is \(\mathrm{ \tau \left( \theta \right) = \theta }\), increasing from \(\mathrm{ \tau \left( 0 \right) =0~to~ \tau \left( 1 \right) =1 }\), the medial correlation coefficient is \(\mathrm{ v \left( 0 \right) =4^{ \left( 1-2^{- \theta } \right) }-1 }\), increasing from \(\mathrm{ v \left( 0 \right) =0~to~v \left( 1 \right) =1 }\). Linear regression of \(\mathrm{ Y }\) in \(\mathrm{ X }\) is, evidently, \(\mathrm{ L_{y}( x \vert \theta ) = \gamma + \theta ( 2- \theta ) ( x- \gamma ) }\). It can be shown that the index of dependence

\( \mathrm {\begin{array}{c}\\ \mathrm{sup} \\ \mathrm{x,y} \end{array} |\Lambda \left( x,y \vert \theta \right) - \Lambda \left( x,y \vert 0 \right) \vert = \left( 1-2^{- \theta } \right) 2^{- \theta / ( 2^{ \theta }-1)} }\)

increases from \(\mathrm{ 0 ( at~ \theta =0 ) to~1/4=.25 ( at~ \theta =1 ) }\). It is, thus intuitive that the distance between the independence \(\mathrm{ (\theta =0 ) }\) and the assumed model for \(\mathrm{ \theta >0 }\) is small in general, and for small samples it will naturally be difficult to separate models.

Owing to the importance of independence, as a reference because in general we can expect association, we can formulate as confirmation of dependence an L.M.P. test of \(\mathrm{ \theta =0~vs~ \theta >0 }\). As for the sample \(\mathrm{\{ ( x_{1},y_{1} ) , \dots, ( x_{n},y_{n} ) \} }\) the likelihood is

\(\mathrm{ L \left( x_{i},y_{i} \vert \theta \right) = \begin{array}{c} \mathrm{n} \\ \mathrm{\pi}\\ \mathrm{1} \end{array}\frac{ \partial ^{2} \Lambda \left( x_{i}, y_{i} \vert \theta \right) }{ \partial \,x_{i}~ \partial\, y_{i}} }\),

the rejection region has the form

\(\mathrm{ \frac{ \partial }{ \partial \, \theta } ( \sum _{1}^{n}\frac{ \partial ^{2}\, \Lambda ( x_{i}, y_{i} \vert \theta ) }{ \partial \,x_{i} \partial \,y_{i}} ) \vert _{ \theta =0}>c_{n} }\)

or

\(\mathrm{ \sum _{1}^{n}v ( x_{i},y_{i} ) \geq c_{n} }\)

where

\(\mathrm{ v ( x,y ) =x+y+x~e^{-x}+y~e^{-y}+ ( e^{-x}+e^{-y}-2 ) ~log( e^{-x}+e^{-y} ) +\frac{1}{e^{-x}+e^{-y}}}\).

In the case of independence, \(\mathrm{ v ( X,Y ) }\)  has mean value zero but infinite variance; evidently \(\mathrm{ 1/n\cdot \sum _{1}^{n}v ( x_{i},y_{i} ) \stackrel{a.s.}\rightarrow0 }\), in the case of independence, by Khintchine’s theorem. Although then the usual Central Limit Theorem is not applicable, it was shown by Tawn (1988) that \(\mathrm{ \sum _{1}^{n}v ( x_{i},y_{i} ) /\sqrt[]{ ( n\,log\,n ) /2} }\) is asymptotically standard normal. Estimation by the maximum likelihood is possible and Cramér — regular for  \(\mathrm{ \theta \neq 0 }\), having for \(\mathrm{ \theta = 0 }\) the difficulties arising from the infinite variance of \(\mathrm{ v ( X,Y ) }\). Thus we have \(\mathrm{ c_{n} \sim \sqrt[]{ ( n\,log\,n ) /2}~N^{-1} ( 1- \alpha ) }\) where \(\mathrm{ N(.) }\) denotes the standard normal distribution function and \(\mathrm{ \alpha }\) the significance level of the one-sided tests.

As \(\mathrm{ \theta }\) is a dependence parameter, it is natural to use the correlation coefficients whose estimators are, as is well known, independent of the location and dispersion parameters of the margins, to test \(\mathrm{ \theta =0~vs.~ \theta >0 }\). It was shown, in Tiago de Oliveira (1964) and (1975), that the most efficient is the classical correlation coefficient. Recall that, in the case of independence, \(\mathrm{ \sqrt[]{n}~r }\) is asymptotically standard normal. As \(\mathrm{ \rho \left( \theta \right) = \theta \left( 2- \theta \right) }\), a point estimator of \(\mathrm{ \theta }\) is given by \(\mathrm{ r= \theta ^{*} ( 2- \theta ^{*} ) ~or~ \theta ^{*}=1-\sqrt[]{1-r} }\), which leads to no truncation of the estimator if \(\mathrm{ r \geq 0 }\) ; otherwise we will take \(\mathrm{ \theta ^{*}=0 }\) (independence).

A quick approach may be the quadrants method (see Annex 5); as the medians of the margins, in reduced form, are \(\mathrm{ \tilde{\mu} =-~log~log~2=0.3665129 }\) we have \(\mathrm{ p \left( \theta \right) =Prob \{ X \leq \tilde{\mu} ,Y \leq \tilde{\mu} \} =exp \left( -2^{1- \theta }~log~2 \right) }\), so the probability that both components are at the same side of the medians, i.e., \(\mathrm{ \{{ X \leq \tilde{\mu} ,Y \leq \tilde{\mu} }\}~ or~ \{ X> \tilde{\mu} ,Y>\tilde{\mu} \} }\) is \(\mathrm{ 2 p \left( \theta \right) }\). If \(\mathrm{ N }\) denotes the number of sample pairs with components both smaller or both larger than the margin medians \(\mathrm{ \tilde{\mu } }\), we know that \(\mathrm{ \sqrt[]{n}\frac{{N}/{n}-2~p \left( \theta \right) }{\sqrt[]{2~p \left( \theta \right) \left( 1-2~p \left( \theta \right) \right) }} }\) is asymptotically standard normal. The point estimator of \(\mathrm{ \theta }\) is then given by

\(\mathrm{ \theta ^{**}=1-\frac{~log~ \left( log \left( 2n/N \right) /log~2 \right) }{log~2} }\)

and the condition \(\mathrm{ 0 \leq \theta ^{*} \leq 1 }\) is equivalent to \(\mathrm{ 1/2 \leq N/n \leq 1 }\) ; as the RHS inequality is trivial, if the LHS one is violated (i.e., \(\mathrm{ n>2N }\) ) we get \(\mathrm{ \theta ^{**}<0 }\) and we must assume independence (if the model is valid).

Notice that \(\mathrm{ p ( \theta ) =exp \{ ( -2^{1- \theta }~log~2 ) \} }\)  increases from  \(\mathrm{ p ( 0 ) =1/4~to~p ( 1 ) =1/2 }\).

It is easy, by the use of the \(\mathrm{ \delta }\)-method, to show that \(\mathrm{ log~2\cdot\sqrt[]{\frac{n~N}{n-N}}~log\frac{2~n}{N} ( \theta ^{**}- \theta ) }\) is asymptotically standard normal.

The logistic model for bivariate minima with standard exponential margins is immediately obtained by the transformation \(\mathrm{ X’=e^{-X},Y’=e^{-Y} A ( u ) =k ( log ( \frac{1-u}{u} ) ) }\) has the following survival function:

\(\mathrm{ S( x’,y’ \vert \theta ) =P ( X’ \geq x’,Y’ \geq y’ ) =exp \{ - [ x’^{{1}/{ ( 1- \theta ) }}+y’^{{1}/{ ( 1- \theta ) }} ] ^{1- \theta } \} }\)

for \(\mathrm{ x’,y’ \geq 0 }\); also we get \(\mathrm{ A ( u \vert \theta ) = ( u^{{1}/{ ( 1- \theta ) }}+ ( 1-u ) ^{{1}/{ ( 1- \theta ) }} ) ^{1- \theta } }\) and the equisurvival median curve is \(\mathrm{ x’^{{1}/{ \left( 1- \theta \right) }}+y’^{{1}/{ \left( 1- \theta \right) }}= \left( log~2 \right) ^{{1}/{ \left( 1- \theta \right) }} }\) ; for \(\mathrm{ \theta =1 }\)  we obtain a straight lines. We also have \(\mathrm{ B( u \vert \theta ) =Prob \{ \frac{Y’}{X’+Y’} \leq u \} =1-D ( log\frac{1-u}{u} \vert \theta ) =\frac{u^{{1}/{ ( 1- \theta ) }}}{u^{{1}/{ ( 1- \theta ) }}+ ( 1-u ) ^{{1}/{ ( 1- \theta ) }}} }\). The function \(\mathrm{ v ( x’,y’ ) is~v ( x’,y’ ) =- ( log~x’+x’~log~x’+log~y’+y’~log~y’ ) + }\)

\(\mathrm{ ( x’+y’-2 ) log ( x’+y’ ) +1/ ( x’+y’ ) }\).

The classical correlation coefficient  \(\mathrm{ \rho = \int _{o}^{1}A^{-2} \left( u \vert \theta \right) ~d~u-1, }\) by the transformation \(\mathrm{ u= ( 1+ ( \frac{1-t}{t} ) ^{1- \theta } ) ^{-1} }\), takes the form \(\mathrm{ \rho \left( \theta \right) = \left( 1- \theta \right) \int _{0}^{1}t^{- \theta } \left( 1-t \right) ^{- \theta }d~t-1= \left( 1- \theta \right) B \left( 1- \theta ,1- \theta \right) -1=}\) \(\mathrm{ ( 1- \theta ) \frac{ \Gamma ^{2} ( 1- \theta ) }{ \Gamma ( 2 ( 1- \theta ) ) }-1 ( as~B ( .,. ) }\) is the beta function), which increases from \(\mathrm{ \rho \left( 0 \right) =0~to~ \rho \left( 1^{-} \right) =1~as~ \left( 1- \theta \right) \Gamma \left( 1- \theta \right) \rightarrow 1~for~ \theta \rightarrow 1^{-} }\). 

Let us now consider the mixed model. Its distribution function is \(\mathrm{ \Lambda ( x,y \vert \theta ) =exp( -e^{-x}-e^{y}+ \theta / ( e^{x}+e^{y} ) ) ,0 \leq \theta \leq 1 }\). When \(\mathrm{ \theta =0 }\) we have independence but for \(\mathrm{ \theta =1 }\)  it does not degenerate, in the diagonal case, as happens with the logistic model.

The model has the dependence function \(\mathrm{ k \left( w \vert \theta \right) =1- \theta \frac{w}{ \left( 1+e^{w} \right) ^{2}}=k \left( -w \vert \theta \right) }\), and so the model is exchangeable, as also follows from

\(\mathrm{ D \left( w \vert \theta \right) =\frac{e^{w}}{1+e^{w}}\cdot\frac{ \left( 1+e^{w} \right) ^{2}- \theta }{ \left( 1+e^{w} \right) ^{2}- \theta ~e^{w}}=1-D \left( -w \vert \theta \right) }\).

The dependence function \(\mathrm{ k \left( w \vert \theta \right) }\) is the \(\mathrm{( 1- \theta , \theta ) }\) mixture of the dependence functions \(\mathrm{ 1 }\) and \(\mathrm{ 1-\frac{e^{w}}{ ( 1+e^{w} ) ^{2}} }\); the inequality \(\mathrm{ \frac{max⁡ ( 1,e^{w} ) }{1+e^{w}} \leq k( w ) \leq 1 }\) shows that the domain \(\mathrm{ [ 0,1] }\), of cannot be enlarged.

The equidistribution curves are drawn as usual. The median curve takes the form \(\mathrm{e^{-x}+e^{-y}- \theta / ( e^{x}+e^{y})=log~2 }\).

The classical correlation coefficient has the expression \(\mathrm{ \rho ( \theta ) =\frac{6}{ \pi ^{2}} \cdot ( arc~cos ( 1-\frac{ \theta }{2} ) ) ^{2} }\) which increases from \(\mathrm{ \rho \left( 0 \right) =0~to~ \rho \left( 1 \right) =2/3 }\), the medial correlation coefficient is \(\mathrm{ v \left( \theta \right) =2^{ \theta /2}-1 }\) increasing from \(\mathrm{ v \left( 0 \right) =0~to~v \left( 1 \right) =1 }\); the model is, thus, not fit to express a strong dependence, as confirmed below by the index of dependence. Linear regression of \(\mathrm{ Y }\) is \(\mathrm{ X }\) is evidently \(\mathrm{ L_{y} ( x \vert \theta ) = \gamma +\frac{6}{ \pi ^{2}} ( arc~cos ( 1-{ \theta }/{2} ) ) ^{2} ( x- \gamma ) }\). We have also

\( \mathrm {\begin{array}{c}\\ \mathrm{sup} \\ \mathrm{x,y} \end{array} \vert \Lambda \left( x,y \vert \theta \right) - \Lambda \left( x \right) y \vert 0 ) \vert =\frac{ \theta }{4- \theta } \left( 1-{ \theta }/{4} \right) ^{4/ \theta }}\)

increasing from \(\mathrm{ 0 \left( at~ \theta =0 \right) to~3^{3}/4^{4}=.105~at=1 }\). The smaller variation of the correlation coefficient and of the index of dependence shows that the deviation from independence is smaller and more difficult to detect.

The LMP test of  \(\mathrm{ \theta =0~vs. \theta >0 }\) is given by the rejection region

\(\mathrm{ \sum _{1}^{n}v \left( x_{i},y_{i} \right) \geq c_{n} }\)

with

\(\mathrm{ v \left( x,y \right) =2\frac{e^{2x} \cdot e^{2y}}{ \left( e^{x}+e^{y} \right) ^{3}}-\frac{e^{2x}+e^{2y}}{ \left( e^{x}+e^{y} \right) ^{2}}+\frac{1}{e^{x}+e^{y}} }\).

As before, the situation is marred by the fact that for independence, although \(\mathrm{ v ( X,Y ) }\) has mean value zero, the variance is infinite. The usual Central Limit Theorem cannot be applied but, as before, \(\mathrm{ 1/n\cdot \sum _{1}^{n}v ( x_{i},y_{i} ) \stackrel{a.s.}\rightarrow0 }\). But as shown also by Tawn (1988), \(\mathrm{ \sum _{1}^{n}v \left( x_{i},y_{i} \right) /\sqrt[]{ \left( n\,log\,n \right) /15} }\) is also asymptotically standard normal; then \(\mathrm{ c_{n} \sim \sqrt[]{ \left( n\,log\,n \right) /15} )~ N^{-1} \left( 1- \alpha \right) }\). Estimation by the maximum likelihood is also possible and Cramér-regular for \(\mathrm{ \theta \neq 0 }\) with analogous difficulties for \(\mathrm{ \theta =0 }\), as before.

Correlation coefficients, independent of the margin parameters, are then to be used to test \(\mathrm{ \theta =0~vs. \theta >0 }\). Once more the usual correlation coefficient is the best. The point estimator of \(\mathrm{ \theta }\), if \(\mathrm{ 0 \leq r \leq 2/3 }\), is given by \(\mathrm{ \theta ^{*}=2 ( 1-cos\sqrt[]{\frac{ \pi ^{2}}{6}}r ) }\) with the natural truncation otherwise.

A quick approach may also use the quadrants method. As the medians of the margins are \(\mathrm{ \tilde{\mu }=-~log~log~2 }\) we have \(\mathrm{ p \left( \theta \right) =Prob \{ X \leq \tilde{\mu} ,Y \leq \tilde{\mu} \} =2^{ \theta /2}/4 }\) which varies from \(\mathrm{ p \left( \theta \right) =1/4~to~p \left( 1 \right) =\sqrt[]{2}/4 }\). If \(\mathrm{ N/n }\) denotes the fraction of sample pairs both smaller or both larger than the margin medians, we know that \(\mathrm{ \sqrt[]{n}~\frac{{N}/{n}-2~p \left( \theta \right) }{\sqrt[]{2~p \left( \theta \right) \left( 1-2~p \left( \theta \right) \right) }} }\) is asymptotically standard normal. If \(\mathrm{ 1/2 \leq N/n \leq \sqrt[]{2}/4 }\) the estimator of \(\mathrm{ \theta }\) is then given by \(\mathrm{ \theta ^{**}=2 ( 1+\frac{log\,N/n}{log\,2}) }\) being truncated as before when the double inequality does not apply and, once more by the \(\mathrm{ \delta }\)-method, we see that \(\mathrm{ \frac{log~2}{2}\sqrt[]{\frac{n~N}{n-N}} ( \theta ^{**}- \theta ) }\) is asymptotically standard normal.

The mixed model for bivariate minima with standard exponential margins immediately has the survival function

\(\mathrm{ S ( x’,y’ \vert \theta ) =exp \{ - ( x’+y’ ) +\frac{ \theta x’+y’}{x’+y’} \} ,x’+y’ \geq 0 }\);

we have \(\mathrm{ A \left( u \vert \theta \right) =1- \theta ~u \left( 1-u \right) }\)and the equisurvival median curve is \(\mathrm{ x’+y’- \theta \frac{x’y’}{x’+y’}=log~2 }\). We have \(\mathrm{ B \left( u \vert \theta \right) =1-D ( log~\frac{1-u}{u} \vert \theta ) =\frac{u- \theta ~u \left( 1-u \right) ^{2}}{1- \theta ~u \left( 1-u \right) } }\). The function \(\mathrm{ v( x’,y’ ) }\) is \(\mathrm{ v \left( x’,y’ \right) =2\frac{x’ \cdot y’}{ \left( x’+y’ \right) ^{3}}-\frac{x’^{2}+y’^{2}}{ \left( x’+y’ \right) ^{2}}+\frac{x’ \cdot y’}{x’+y’} }\). The classical correlation coefficient is simply

\(\mathrm{ \rho \left( \theta \right) =\frac{8}{ \left( 4- \theta \right) ~\sqrt[]{4~ \theta - \theta ^{2}}}tg^{-1}\sqrt[]{\frac{ \theta }{4- \theta }}-\frac{2- \theta }{4- \theta } }\)

increasing from \(\mathrm{ \rho ( 0 ) =0~to~ \rho ( 1 ) =\frac{1}{3} ( \frac{4~ \pi }{3\sqrt[]{3}}-1 ) =.4727997 }\).

Finally, we discuss the statistical choice between the two non-separated models (logistic and mixed), non-separated because they have in common the independence case ( \(\mathrm{ \theta =0 }\)  in both cases). As developed in Tiago de Oliveira (1984) and generalized in Tiago de Oliveira (1985), the decision rule (with \(\mathrm{ v_L }\) and \(\mathrm{ v_M }\) denoting the \(\mathrm{ v }\) functions for the logistic and mixed models, given before) is:

decide for independence if            \(\mathrm{ \sum _{1}^{n}v_{L} ( x_{i},y_{i} ) \leq a_{n}, \sum _{1}^{n}v_{M} ( x_{i},y_{i} ) \leq b_{n} }\);  

decide for the logistic model if       \(\mathrm{ \sum _{1}^{n}v_{L} ( x_{i},y_{i}) >a_{n} }\)

and                                                 \(\mathrm{ \sum _{1}^{n}\frac{v_{L} ( x_{i},y_{i} ) }{a_{n}}> \sum _{1}^{n}\frac{v_{M} ( x_{i},y_{i} ) }{b_{n}} }\);

decide for the mixed model if         \(\mathrm{ \sum _{1}^{n}v_{M} ( x_{i},y_{i} ) >b_{n} }\)

and                                                 \(\mathrm{ \sum _{1}^{n}\frac{v_{M} ( x_{i},y_{i} ) }{b_{n}}>\sum _{1}^{n}\frac{v_{L} ( x_{i},y_{i} ) }{a_{n}} }\).

It was shown, in Tawn (1988), that the asymptotic distribution of the pair \(\mathrm{ ( \sum _{1}^{n}v_{L} ( x_{i},y_{i} ) /\sqrt[]{ ( n~log~n ) /2}, \sum _{1}^{n}v_{M} ( x_{i},y_{i} ) /\sqrt[]{ ( n~log~n ) /15 ) } }\) is asymptotically binormal independent with standard margins, and so we should take, with the usual meanings, \(\mathrm{ a_{n} \sim \sqrt[]{ ( n~log~n ) /2}~N^{-1} ( 1- \alpha /2 ) ~and~b_{n} \sim \sqrt[]{ ( n~log~n ) /15}~N^{-1} ( 1- \alpha /2 ) }\).

Let us now make some comments.

It must be strongly emphasized that the tests and estimators based on the correlation coefficients are independent of the values of the margin parameters. However, in the other cases we have to “estimate” reduced values by using \(\mathrm{ \hat{x}_{i}= ( x_{i}- \hat{\lambda} _{x} ) / \hat{\delta }_{x}~and~ \hat{y}_{i}= ( y_{i}- \hat{\lambda} _{y} ) / \hat{ \delta} _{y} }\) in the Gumbel margins case, the effect of which is not known.

It is natural to consider, as a first step in the statistical analysis of random pairs of maxima or minima, a test for independence, although in many cases we may expect its rejection, maxima or minima being almost always dependent, as observed above.

We can now summarize the results for a test of independence given in Tiago de Oliveira (1965).

The comparison of the usual estimators of the classical, difference-sign, grade, and medial correlation coefficients for one-sided tests, owing to the non-negativity of the values of those coefficients for the bivariate extremal distributions, shows, if we are in the model with bivariate margins, that:

a) the difference-sign and grade correlation are asymptotically equivalent as tests for dependence with respect to any (reasonable) system of alternatives considered, and as the Spearman correlation is always asymptotically equivalent to the grade correlation, the three are equivalent in this case; note that as the grade correlation can be optimal for testing for independence — see Tiago de Oliveira (1985) — with \(\mathrm{ F ( x,y \vert \theta ) }\) such that \(\mathrm{ F \left( x,+ \infty \vert \theta \right) =F_{1} \left( x \right) ,F \left( + \infty,y \vert \theta \right) =F_{2} \left( y \right) ,F \left( x,y \vert 0 \right) =F_{1} \left( x \right) F_{2} \left( x \right) }\) only if \(\mathrm{ \frac{ \partial ~F \left( x,y \vert \theta \right) }{ \partial ~ \theta } \vert _{ \theta =0}=F_{1} \left( x \right) \left( 1-F_{1} \left( x \right) \right) F_{2} \left( y \right) \left( 1-F_{2} \left( y \right) \right) }\) which is not the case; no one of the correlation coefficients tests above can be optimal;

b) for the two differentiable models it was shown that the classical correlation coefficient

\(\mathrm{ r= \sum _{}^{} ( x_{i}-\bar{x} ) ( y_{i}-\bar{y} ) / \{ \sum _{}^{} ( x_{i}-\bar{x} ) ^{2} ( y_{i}-\bar{y} ) ^{2} \} ^{1/2} }\)

gives the best test of the ones above based on correlation coefficients; as \(\mathrm{ \sqrt[]{n}~r }\) is asymptotically standard normal we will accept the hypothesis of independence if \(\mathrm{ \sqrt[]{n}~r \leq \lambda _{ \alpha } }\) where \(\mathrm{ N \left( \lambda _{ \alpha } \right) =1- \alpha }\);

c) the Gumbel and Mustafi (1967 test is equivalent to use of medial correlation;

d) the strip method devised by Posner, Rodemich, Ashlock and Lurie (1969) also has an efficiency smaller than that of the use of the classical product-moment method.

Evidently, some of those results — the non-parametric ones — are immediately extended.

There exist asymmetric extensions of the logistic model; see Tawn (1988). Also there exists another asymmetric model under study, presented in Tiago de Oliveira (1989c), with high correlations \(\mathrm{ \left( \rho >.77 \right) }\); it has a non-differentiable extension such that \(\mathrm{ 0 \leq \rho \leq 1 }\). It gave a good fit to Fox River.

3 . The non-differentiable models: statistical decision

We know essentially three forms of the non-differentiable models — i.e., with a singular component, thus not having a planar density: the Gumbel model, the biextremal model, and the natural model, the last one being recently generalized by the author (see the next chapter for some details); but no statistical results have appeared so far except when it degenerates into the natural model.

The Gumbel model has the distribution function is \(\mathrm{ \Lambda ( x,y \vert \theta ) =exp\,\{-(e^{-x}-e^{-y} -\theta \,\,min(e^{-x},e^{-y}) \},0\leq\theta\leq1 }\) ; its dependence function is \(\mathrm{ k ( w \vert \theta ) =1- \theta \frac{min⁡ ( 1,e^{-w} ) }{1+e^{-w}} }\) which is the \(\mathrm{ \left( 1- \theta , \theta \right) }\) mixture of \(\mathrm{ k_{1}( w ) =1 }\) (independence) and \(\mathrm{ k_{D} ( w ) =\frac{max⁡ ( 1,e^{-w}) }{1+e^{-w}} }\) (diagonal case). The model is exchangeable as \(\mathrm{ k \left( w \vert \theta \right) =k \left( -w \vert \theta \right) }\) and we have \(\mathrm{ D \left( w \vert \theta \right) =\frac{1- \theta }{1- \theta +e^{-w}} }\) if \(\mathrm{ w<0 }\), \(\mathrm{ D \left( w \vert \theta \right) =\frac{e^{w}}{1- \theta +e^{w}} }\)  if \(\mathrm{ w>0 }\),  with a jump of \(\mathrm{ \frac{ \theta }{2- \theta } }\), at \(\mathrm{ w=0 }\), which is probability \(\mathrm{ P ( Y=X ) }\); if \(\mathrm{ w \neq 0 }\) we have \(\mathrm{ D \left( w \vert \theta \right) +D \left( -w \vert \theta \right) =1 }\). This model was introduced in Tiago de Oliveira (1971) and is a conversion, to Gumbel margins, of a bivariate model with exponential margins, due to Matshall and Olkin (1967).

The median equidistribution curve is, evidently, \(\mathrm{ e^{-x}+e^{-y}- \theta\,min \left( e^{-x},e^{-y} \right) =log\,2 }\). The classical correlation coefficient can take the form \(\mathrm{ \rho \left( \theta \right) =\frac{12}{ \pi ^{2}} \int _{0}^{ \theta }\frac{log \left( 2-t \right) }{1-t}d~t }\), the difference-sign correlation coefficient is \(\mathrm{ \tau \left( \theta \right) = \theta / \left( 2- \theta \right) }\), the grade correlation coefficient is \(\mathrm{ \chi \left( \theta \right) =\frac{3~ \theta }{4- \theta } }\), and the medial correlation coefficient is \(\mathrm{ v \left( \theta \right) =2^{ \theta }-1 }\), all the coefficients increasing from \(\mathrm{ 0 }\) to \(\mathrm{ 1 }\) as \(\mathrm{ \theta }\) varies from \(\mathrm{ 0 }\) to \(\mathrm{ 1 }\). The index of dependence is

\( \mathrm {\begin{array}{c} \\ \mathrm{sup} \\ \mathrm{x,y} \end{array} \vert \Lambda \left( x,y \vert \theta \right) - \Lambda \left( x \right) y \vert 0 ) =\frac{ \theta }{2- \theta } \left( 1- \theta /2 \right) ^{2/ \theta } }\)

increasing also from \(\mathrm{ 0 }\) to \(\mathrm{ 1/4 }\) as \(\mathrm{ \theta }\)  increases from \(\mathrm{ 0 }\) to \(\mathrm{ 1 }\).

The use of formulae given for curvilinear regression leads to

\(\mathrm{ \bar{y} \left( x \vert \theta \right) = \gamma +log \left( 1- \theta \right) + \int _{ \left( 1- \theta \right) ^{e^{-x}}}^{+ \infty}{e^{-t}}/{t} \cdot d~t- \left( 1- \theta \right) exp \left( \theta ~e^{-x} \right) ~ \int _{e^{-x}}^{+ \infty}{e^{-t}}/{t} \cdot d~t }\)

\(\mathrm{ =x+ \int _{0}^{ \left( 1- \theta \right) ^{e^{-x}}}\frac{1-e^{-t}}{t}d~t- \left( 1- \theta \right) exp \left( \theta ~e^{-x} \right) \int _{e^{-x}}^{+ \infty}\frac{e^{-t}}{t}d~t }\).

The form of \(\mathrm{ \bar{y} \left( x \vert \theta \right) }\) is very different from the linear regression \(\mathrm{ \gamma + \rho \left( \theta \right) \left( x- \gamma \right) }\), but the non-linearity index NL oscillates between \(\mathrm{ 0 }\) (from \(\mathrm{ \theta =0 }\) to \(\mathrm{ \theta =1 }\) ) and .006; the improvement is very small in the passage from the (simple) linear regression to the (more complex) curvilinear regression: in fact, a large difference between the two regressions only happens for very small or very large values of \(\mathrm{ X }\), whose probability is very small: see Tiago de Oliveira (1974).

Let us now consider statistical decision for this model.

The estimation of the margin parameters is avoided using the correlation coefficients. The test for independence can be done, as before, based on the fact that if \(\mathrm{ \theta =0 }\) (independence) \(\mathrm{ \sqrt[]{n}~r }\) is asymptotically standard normal.

The method of quadrants can be used, once more with \(\mathrm{ p \left( \theta \right) =Prob \{{ X \leq \tilde{\mu} ,Y \leq \tilde{\mu} }\} ={2^{ \theta }}/{4}=2^{ \theta -2} }\) , increasing from \(\mathrm{ p \left( 0 \right) =1/4~to~p \left( 1 \right) =1/2 }\) . Then \(\mathrm{ \sqrt[]{n}~\frac{{N}/{n}-2~p \left( \theta \right) }{\sqrt[]{2~p \left( \theta \right) \left( 1-2~p \left( \theta \right) \right) }} }\) is asymptotically standard normal and \(\mathrm{ \theta ^{**}=\frac{log~N/n}{log\,2}+1 }\) (supposing that \(\mathrm{ N/n \geq 1/2 }\) ) is an estimator of \(\mathrm{ \theta }\) , with truncation if necessary. Also it is immediate to show that \(\mathrm{ \sqrt[]{n}~\frac{log2}{\sqrt[]{2^{1- \theta }-1}} \left( \theta ^{**}- \theta \right) }\) and \(\mathrm{ \left( log~2 \right) \sqrt[]{\frac{n~N}{n-N}} \left( \theta ^{**}- \theta \right) }\) are asymptotically standard normal.

The Gumbel model for bivariate minima with standard exponential margins was presented, in a different approach, as said before, by Marshall and Olkin (1967). We have \(\mathrm{ S ( x’,y’ \vert \theta ) =e^{-x’-y’+ \theta\, min⁡ ( x’,y’ ) },0 \leq \theta \leq 1 }\); we have \(\mathrm{ A \left( u \vert \theta \right) =1- \theta ~min \left( u,1-u \right) a~mixture \left( \theta ,1- \theta \right) }\)of independence \(\mathrm{ \left( A \left( u \right) =1 \right) }\) and diagonal case \(\mathrm{ \left( A \left( u \right) =max \left( u,1-u \right) \right) }\) as before, naturally; and the equisurvival median curve is \(\mathrm{ x’+y’- \theta ~min \left( x’,y’ \right) =log~2 }\).

We have \(\mathrm{B \left( u \vert \theta \right) =1-D \left( log\frac{1-u}{u} \vert \theta \right) =\frac{ \left( 1- \theta \right) u}{1- \theta ~u}if~u<\frac{1}{2},B \left( u \vert \theta \right) =\frac{u}{1- \theta ~\left( 1-u \right) }if~u \geq \frac{1}{2} }\)with a jump of \(\mathrm{ \theta / \left( 2- \theta \right) }\) at \(\mathrm{ u=1/2 }\), as could be expected from the result for \(\mathrm{ D \left( w \vert \theta \right) }\).

The classical correlation coefficient is \(\mathrm{ \rho \left( \theta \right) = \theta / \left( 2- \theta \right) }\) increasing from \(\mathrm{ \rho \left( 0 \right) =0~to~ \rho \left( 1 \right) =1 }\).

Let us now consider the biextremal model. It appears naturally in the study of extremal processes; see Tiago de Oliveira (1968).

It can defined directly as follows: consider \(\mathrm{ ( X,Z ) }\) independent reduced Gumbel random variables and form the new pair \(\mathrm{( X,Y ) }\) with \(\mathrm{ Y=max ( X-a,Z-b ) }\), \(\mathrm{ a }\) and \(\mathrm{ b }\) such that \(\mathrm{ Prob \{{ Y \leq y }\} =Prob \{ max ( X-a,Z-b ) \leq y \} = \Lambda ( y ) }\), which implies \(\mathrm{ e^{-a}+e^{-b}=1 }\). Putting \(\mathrm{ e^{-a}= \theta ,e^{-b}=1- \theta \left( 0 \leq \theta \leq 1 \right) }\), we have \(\mathrm{ Y=max \left( X+log \,\theta ,Z+log \left( 1- \theta \right) \right) }\). This was the initial use of the max-technique to generate new models and is a particular case of a natural model, to be defined below, with \(\mathrm{ a=0,b=+ \infty,c=-log \,\theta \,\,and~d=-log ( 1- \theta ) }\). The distribution function is \(\mathrm{ \Lambda \left( x,y \vert \theta \right) =exp \{{ -max \left( e^{-x}+ \left( 1- \theta \right) e^{-y},e^{-y} \right) }\} , \left( 0 \leq \theta \leq 1 \right) }\)and the dependence function is immediately \(\mathrm{ k \left( w \vert \theta \right) =1-{min \left( \theta ,e^{w} \right) }/{ \left( 1+e^{w} \right) }= \left( 1- \theta +max \left( \theta ,e^{w} \right) \right) / \left( 1+e^{w} \right) }\)and we have \(\mathrm{ D \left( w \vert \theta \right) =0~if~w<log~ \theta ~and~ \left( 1+ \left( 1- \theta \right) e^{-w} \right) ^{-1} }\) if \(\mathrm{ w \geq log~ \theta }\) with a jump of \(\mathrm{ \theta }\) at \(\mathrm{ w = log~ \theta }\) , which is equal to the \(\mathrm{ Prob \{ Y=X+log~ \theta \} }\), and so the singular part is concentrated at the line \(\mathrm{ y=x+log~ \theta }\), with probability \(\mathrm{ \theta }\). Evidently \(\mathrm{ Y \geq X+log~ \theta ~ }\)with probability one. For \(\mathrm{ ~ \theta =0~and~ \theta =1 }\) we obtain independence and the diagonal case. Notice that  \(\mathrm{ k \left( w \vert \theta \right) \neq k \left( -w \vert \theta \right) }\) and so the model is not exchangeable.

The classical correlation coefficient of the biextremal model is

\(\mathrm{ \rho \left( \theta \right) =\frac{6}{ \pi ^{2}} \int _{0}^{ \theta }\frac{log~t}{1-t}d~t=1+\frac{6}{ \pi ^{2}}L \left( \theta \right) }\),

\(\mathrm{ L \left( \theta \right) = \int _{1}^{ \theta }\frac{log~w}{w-1}d~w }\) being the Spence integral; \(\mathrm{ p \left( \theta \right) }\) increases from \(\mathrm{ p \left( 0 \right) =0~to~p \left( 1 \right) =1 }\). Linear regression is thus \(\mathrm{ L_{y} \left( x \right) = \gamma + \rho \left( \theta \right) \left( x- \gamma \right) }\).

The grade correlation, difference-sign correlation, and medial correlation coefficients are easily seen to be \(\mathrm{ \chi \left( \theta \right) =3 \theta / \left( 2+ \theta \right) , \tau \left( \theta \right) = \theta ~and~v \left( \theta \right) =2^{ \theta }-1 }\) , all increasing from \(\mathrm{ 1 }\) to \(\mathrm{ \theta }\)  as \(\mathrm{ 0 }\) varies from \(\mathrm{ 1 }\).

The index of dependence is

\( \mathrm {\begin{array}{c}\\ \mathrm{sup} \\ \mathrm{x,y} \end{array} \vert \Lambda \left( x,y \vert \theta \right) - \Lambda \left( x,y \vert 0 \right) \vert = \theta / \left( 1+ \theta \right) ^{{ \left( 1+ \theta \right) }/{ \theta }} }\),

increasing from \(\mathrm{ 0~to~1/4~as~ \theta }\)  increases from \(\mathrm{ 0~to~1 }\).

The curvilinear regressions — here we do not have exchangeability — are

\(\mathrm{ \bar{y} \left( x \vert \theta \right) =x+log \, \theta + \int _{0}^{{ \left( 1- \theta \right) }/{ \theta \cdot e^{-x}}}\frac{1-e^{-t}}{t}d~t= }\)

\(\mathrm{ =x+log \left( 1- \theta \right) + \int _{\frac{1- \theta }{ \theta }e^{-x}}^{+ \infty}\frac{e^{-t}}{t}d~t }\)

and

\(\mathrm{ \bar{x} \left( y \vert \theta \right) =y-log\, \theta - \left( 1- \theta \right) exp \left( \theta\, e^{- \bar{x}} \right) \int _{ \theta \,e^{-x}}^{+ \infty}\frac{e^{-t}}{t}d~t }\). 

The comparison with the linear regressions \(\mathrm{ L_{y} \left( x \vert \theta \right) = \gamma + \rho \left( \theta \right) \left( x-y \right) }\) and \(\mathrm{ L_{x} \left( y \vert \theta \right) = \gamma + \rho \left( \theta \right) \left( x- \gamma \right) }\)shows that the non-linearity index varies from \(\mathrm{ 0~to~0.068~ }\)in the first curvilinear regression and from \(\mathrm{ 0~to~0.007 }\) in the second one; see Tiago de Oliveira (1974). We do not, therefore, have a substantial improvement in using curvilinear regression, although the comparison between \(\mathrm{ \bar{y} \left( x \vert \theta \right) }\) and \(\mathrm{ L_{y} \left( x \vert \theta \right) ~ }\)and between \(\mathrm{ \bar{x} \left( y \vert \theta \right) }\) and \(\mathrm{ L_{x} \left( y \vert \theta \right) }\) shows a great difference. The explanation is as before: the small probability mass attached to the regions where the regressions (linear and curvilinear) differ strongly.

The test of independence \(\mathrm{ \left( \theta =0~vs.~ \theta >0 \right) }\) can be made by using the classical correlation coefficient as usual.

The quadrants method can be applied as \(\mathrm{ p \left( \theta \right) =Prob \{ X \leq \tilde{\mu} ,Y \leq \tilde{\mu } \} =2^{ \theta }/4 }\), equal to the expression for the Gumbel model; we can use the techniques described before.

As regards the estimation of \(\mathrm{ \theta ,as~y_{i}-x_{i} \geq log \,\theta }\) and \(\mathrm{ 0 \leq \theta \leq 1 }\), the best estimator of the left-end point \(\mathrm{ log~ \theta }\) is \(\mathrm{ log \,\tilde{\theta }=min \left( y_{i}-x_{i},0 \right) }\) and so \(\mathrm{ \tilde{\theta} =e^{min \left( y_{i}-x_{i},0 \right) }=min \left( e^{y_{i}-x_{i}},1 \right) }\). We have, for \(\mathrm{ z<1,Prob \{{ \tilde{ \theta} \leq z \vert \theta }\} =Prob \{ \begin{array}{c} \mathrm{n} \\ \mathrm{min} \\ \mathrm{1} \end{array} \left( e^{y_{i}-x_{i}} \right) \leq z \vert \theta \} }\) \(\mathrm{ =1-Prob \{ \begin{array}{c} \mathrm{n} \\ \mathrm{min} \\ \mathrm{1} \end{array} \left( y_{i}-x_{i} \right) >log\,z \vert \theta \} =1- \left( 1-D \left( log\,z \vert \theta \right) \right) ^{n} }\) and so have \(\mathrm{ Prob \{{ \tilde{\theta} \leq z \vert \theta }\} =0~if~z< \theta ,Prob \{{ \tilde{\theta} \leq z \vert \theta }\} =1- ( \frac{z}{1- \theta +z} ) ^{n}~if~ \theta \leq z<1 }\) and \(\mathrm{ Prob \{{ \tilde{\theta} \leq z \vert \theta }\} =1~if~z \geq 1 }\); the jumps have the values \(\mathrm{ 1- \left( 1- \theta \right) ^{n}~at~z=0~and~ ( \frac{1- \theta }{2- \theta } ) ^{n}~at~z=1;if~ \theta \neq 1 }\)we see that \(\mathrm{ \tilde{\theta} {\stackrel {P} {\rightarrow \theta }}}\) ; see Tiago de Oliveira (1970) and (1975b) for other details.

The biextremal model for minima has the survival function \(\mathrm{ S \left( x’,y’ \vert \theta \right) =e^{-max \left( x’+ \left( 1- \theta \right) y’,y’ \right) },0 \leq \theta \leq 1,x’,y’ \geq 0 }\); we obtain \(\mathrm{ A \left( u \vert \theta \right) =max \left( u,1- \theta ~u \right) ~and~B \left( u \vert \theta \right) =1-D \left( log\frac{1-u}{u} \vert \theta \right) =\frac{ \left( 1- \theta \right) u}{1- \theta ~u}if~u<\frac{1}{1+ \theta } }\), \(\mathrm{ B \left( u \vert \theta \right) =1~if~u \geq \frac{1}{1+ \theta } }\) with a jump of \(\mathrm{ \theta ~at~u= {1}/{ \left( 1+ \theta \right) } }\). The equisurvival median curve is \(\mathrm{ max \left( x’+ \left( 1- \theta \right) y’,y’ \right) =y’+max \left( x’- \theta y’,0 \right) =log~2 }\). The classical correlation coefficient takes the value \(\mathrm{ \rho \left( \theta \right) = \theta }\).

Let us now consider the natural model. Using the max-technique, from a pair \(\mathrm{ \left( \bar{X},\bar{Y} \right) }\) of independent reduced Gumbel random variables we get a new (dependent) random pair \(\mathrm{ X=max \left( \bar{X}-a,\bar{Y}-b \right) ,Y=max \left( \bar{X}-c,\bar{Y}-d \right) }\) with Gumbel margins if and only if \(\mathrm{ e^{-a}+e^{-b}=1~and~e^{-c}+e^{-d}=1 }\). The dependence function is \(\mathrm{ k \left( w \vert a,b,c,d \right) = ( e^{-min \left( a,c+w \right) }+e^{-min \left( b,d+w \right) } ) \left( 1+e^{w} \right) }\), whose tails coincide with those of the diagonal case \(\mathrm{ k_{D} \left( w \right) =\frac{max \left( 1,e^{w} \right) }{1+e^{w}} }\). It is easy to see that \(\mathrm{ W=Y-X=min ( a+\bar{W},b ) -min( c+\bar{W},d ) }\) is such that \(\mathrm{ a-c \leq W \leq b-d~if~a+d<b+c~and~b-d \leq W \leq a-c~if~b+c<a+d }\). The case \(\mathrm{ a+b=b+c }\) is irrelevant, because \(\mathrm{ W }\) is then an almost sure random variable with \(\mathrm{ W=a-c=b-d }\), and as \(\mathrm{ W }\) must be zero (its mean value is zero) we have \(\mathrm{ a=c,b=d,X=Y }\).

From now on we will suppose \(\mathrm{ a+d<b+c }\), the case \(\mathrm{ b+c<a+d }\) being dealt with by exchange of \(\mathrm{ X }\) and \(\mathrm{ Y }\). Let us put \(\mathrm{ \alpha =a-c,b-d= \beta }\); by the fact that \(\mathrm{ \bar{W } }\) has a mean value of zero we have \(\mathrm{ \alpha <0< \beta }\). Using the relations  \(\mathrm{ e^{-a}+e^{-b}=1~and~e^{-c}+e^{-d}=1 }\), with \(\mathrm{ a-c= \alpha ~and~b-d= \beta }\) we get the final expression

\(\mathrm{ k \left( w \vert \alpha , \beta \right) =\frac{ \left( e^{ \beta }-1 \right) ~max \left( 1,e^{ \alpha -w} \right) + \left( 1-e^{ \alpha } \right) ~max \left( 1,e^{ \beta -w} \right) }{ \left( e^{ \beta }-e^{ \alpha } \right) ~ \left( 1+e^{-w} \right) } }\)

or, in more detail,

\(\mathrm{ k \left( w \vert \alpha , \beta \right) =\frac{1}{1+e^{w}}~if~w< \alpha }\),

\(\mathrm{ =\frac{e^{ \beta }-1+ \left( 1-e^{ \alpha } \right) ~e^{ \beta -w}}{ \left( e^{ \beta }-e^{ \alpha } \right) ~ \left( 1+e^{-w} \right) } }\)      if   \(\mathrm{ \alpha \leq w \leq \beta , }\)

  \(\mathrm{ =\frac{1}{1+e^{w}} }\)                   if   \(\mathrm{ \beta <w }\);

as said before for \(\mathrm{ ​​​w< \alpha }\) and for \(\mathrm{ w> \beta }\) we have \(\mathrm{ ​​​​k \left( w \vert \alpha , \beta \right) =k_{D} \left( w \right) }\).

Note that \(\mathrm{k \left( w \vert - \beta ,- \alpha \right) =k \left( -w \vert \alpha , \beta \right) }\); we get an exchangeable model if \(\mathrm{\alpha =- \beta }\) . As said before note also that for \(\mathrm{\alpha =log~ \theta , \beta =+ \infty }\), we get the biextremal model and \(\mathrm{ \alpha =- \infty, \beta =-log \,\theta }\) its dual.

Using now the parameters \(\mathrm{\left( \alpha , \beta \right)}\) ( \(\mathrm{ \alpha <0< \beta , \alpha \leq W \leq \beta }\) with probability one), we see that \(\mathrm{\Lambda \left( x,y \vert \alpha , \beta \right) =exp⁡ \{ -[ \left( e^{ \beta }-1 \right) max \left( e^{-x},e^{ \alpha -y} \right) + \left( 1-e^{ \alpha } \right) max~ \left( e^{-x},e^{ \beta -y} \right) ] /}\)  \(\mathrm{\left( e^{ \beta }-e^{ \alpha } \right) \} , - \infty< \alpha \leq 0 \leq \beta <+ \infty }\). \(\mathrm{D \left( w \vert \alpha , \beta \right) }\) is given by:

\(\mathrm{D \left( w \vert \alpha , \beta \right) =0 }\)               if   \(\mathrm{ w< \alpha }\)

\(\mathrm{ =\frac{e^{ \beta }-1}{e^{ \beta }-1+ \left( 1-e^{ \alpha } \right) ~e^{ \beta -w}} }\)           if   \(\mathrm{ \alpha \leq w \leq \beta }\), 

\(\mathrm{ =1 }\)                                if   \(\mathrm{ \beta \leq w }\);

with jumps of \(\mathrm{ \frac{ ( e^{ \beta }-1 ) e^{ \alpha }}{e^{ \beta }-e^{ \alpha }} }\) at \(\mathrm{ w= \alpha }\) and of \(\mathrm{ \frac{ \left( 1-e^{ \alpha } \right) }{e^{ \beta }-e^{ \alpha }} }\) at \(\mathrm{ w= \beta ; \left[ \alpha , \beta \right] }\) is, thus, the support of \(\mathrm{ W=Y-X,R= \beta - \alpha }\) the range and so when we have \(\mathrm{ - \infty< \alpha \leq 0 \leq \beta <+ \infty }\) there exists, as \(\mathrm{ \alpha \leq W=Y-X \leq \beta }\), a strong stochastic connection between \(\mathrm{X }\) and \(\mathrm{ Y }\); this model thus has flexibility to express various situations.

The classical correlation coefficient can be written is

\(\mathrm{ \rho \left( \alpha , \beta \right) =1-\frac{6}{ \pi ^{2}} \int _{- \infty}^{+ \infty}log\frac{k \left( w \vert \alpha , \beta \right) }{k_{D} \left( w \right) }d~w= }\)

\(\mathrm{ =1-\frac{6}{ \pi ^{2}} \int _{ \alpha }^{ \beta }log\frac{e^{ \beta }-1+ \left( 1-e^{ \alpha } \right)\, e^{ \beta -w}}{ \left( e^{ \beta }-e^{ \alpha } \right)\, max \,\left( 1+e^{-w} \right) }d~w }\);

it is immediate that \(\mathrm{\rho \left( \alpha , \beta \right) \rightarrow 1 }\) as  \(\mathrm{ \alpha , \beta \rightarrow 0 }\) and that \(\mathrm{ \rho \left( \alpha , \beta \right) \rightarrow 0~as- \alpha , \beta \rightarrow + \infty }\). 

The difference-sign correlation coefficient is

\(\mathrm{ \tau \left( \alpha , \beta \right) =1- \int _{ \alpha }^{ \beta }D \left( w \vert \alpha , \beta \right) ( 1-D \left( w \vert \alpha , \beta \right)) d~w= }\)

\(\mathrm{ =1-\frac{ \left( e^{ \beta }-1 \right) \left( 1-e^{ \alpha } \right) }{ \left( e^{ \beta }-e^{ \alpha } \right) }=\frac{e^{ \beta + \alpha }-2~e^{ \alpha }+1}{e^{ \beta }-e^{ \alpha }} }\);

the medial correlation coefficient takes the form

\(\mathrm{ v \left( \alpha , \beta \right) =2^{ \left( e^{ \beta + \alpha }-2~e^{ \alpha }+2 \right) / \left( e^{ \beta }-e^{ \alpha } \right) }-1 }\).

The equidistribution median curve is

\(\mathrm{ ( e^{ \beta }-1) max ( e^{-x},e^{ \alpha -y}) +( 1-e^{ \alpha } ) max ( e^{-x},e^{ \beta -y} ) =log~2 \cdot ( e^{ \beta }-e^{ \alpha } ) }\).

The index of dependence is

\(\mathrm{ \left( 1-k_{0} \right) exp \{ \frac{k_{0}}{1-k_{0}}log~k_{0} \}~ where~k_{0}=min ( \frac{1}{e^{ \alpha }+1},\frac{e^{ \beta }}{e^{ \beta }+1} ) }\).

The linear regressions with reduced margins are given by the straight lines \(\mathrm{ L_{y} \left( x \right) = \gamma + \rho \left( \alpha , \beta \right) \left( x- \gamma \right) ~and~L_{x} \left( y \right) = \gamma + \rho \left( \alpha , \beta \right) \left( y- \gamma \right) }\). For the curvilinear regression, as the exchange of \(\mathrm{X }\) and \(\mathrm{ Y }\) is equivalent to the exchange of \(\mathrm{ -\alpha }\) and \(\mathrm{\beta }\), it is sufficient to compute the curvilinear regression

\(\mathrm{ ​​​ \bar{y} \left( x \vert \alpha , \beta \right) =x+ \beta -\frac{ \left( e^{ \beta }-1 \right) }{e^{ \beta }-e^{ \alpha }}exp ( \frac{1-e^{ \alpha }}{e^{ \beta }-e^{ \alpha }}e^{-x} ) \cdot \int _{\frac{1-e^{ \alpha }}{e^{ \beta }-e^{ \alpha }}e^{-x}}^{\frac{e^{ \beta }}{e^{ \alpha }}\frac{1-e^{ \alpha }}{e^{ \beta }-e^{ \alpha }}e^{-x}}\frac{e^{-t}}{t}d~t }\).

Analogously to the situation for the biextremal model, we may expect that regression curves are a good approximation, also in the same sense, to the curvilinear regression.

As \(\mathrm{\alpha \leq w_{i} \leq \beta }\), with \(\mathrm{\alpha <0< \beta }\), we can take as estimators of \(\mathrm{ \alpha }\)  and \(\mathrm{\beta }\) the statistics \(\mathrm{\alpha ^{*}=min \left( 0,y_{i}-x_{i} \right) and~ \beta ^{*}=min \left( 0,y_{i}-x_{i} \right) }\). We have exchangeability \(\mathrm{\left( k \left( w \vert \alpha , \beta \right) =k \left( -w \vert \alpha , \beta \right) \right) }\) iff

\(\mathrm{ k \left( w \vert \alpha , \beta \right) =\frac{1}{1+e^{w}} }\)        if       \(\mathrm{ w<- \beta }\)

\(\mathrm{ =\frac{e^{ \beta }}{1+e^{ \beta }} }\)                         if      \(\mathrm{ - \beta \leq w \leq \beta }\)  

\(\mathrm{ =\frac{1}{1+e^{-w}} }\)                        if      \(\mathrm{ \beta \leq w }\)

and, consequently

  \(\mathrm{ D \left( w \vert \beta \right) ~=0 }\)    if     \(\mathrm{ w< \beta }\)

\(\mathrm{ =\frac{1}{1+e^{-w}} }\)              if     \(\mathrm{ - \beta \leq w \leq \beta }\)  

\(\mathrm{ =1 }\)                     if     \(\mathrm{ \beta \leq\alpha }\),

with jumps of \(\mathrm{ \left( 1+e^{ \beta } \right) ^{-1} }\) at \(\mathrm{ - \beta }\) and \(\mathrm{ \beta }\).

Also we have

\(\mathrm{ \rho \left( \beta \right) =1-\frac{12}{ \pi ^{2}} \int _{0}^{ \beta }\frac{w}{1+e^{-w}}d~w }\).

\(\mathrm{ \tau \left( \beta \right) =\frac{2}{1+e^{ \beta }} }\).

\(\mathrm{ v \left( \beta \right) =2^{2/ \left( 1+e^{ \beta } \right) }-1 }\),

\(\mathrm{{ \bar{y} {\left( x \vert \beta \right)}} =x+ \beta -\frac{e^{ \beta }}{1+e^{ \beta }}exp\frac{e^{-x}}{1+e^{ \beta }} \int _{\frac{e^{-x}}{1+e^{ \beta }}}^{\frac{e^{ \beta -x}}{1+e^{ \beta }}}\frac{e^{-t}}{t}~d~t }\).

and \(\mathrm{ \beta ^{*}=max \vert x_{i} \vert }\) as estimator of \(\mathrm{ \beta }\).

Passing to standard exponential margins and introducing the new parameters \(\mathrm{ \alpha ’=e^{- \alpha } \geq 1 \geq \beta ’=e^{- \beta } }\), we get

\(\mathrm{ S \left( x’,y’ \vert \alpha ’, \beta ’ \right) =exp⁡ \{ -\frac{ \left( 1- \beta ’ \right) max \left( \alpha ’x’,y’ \right) + \left( \alpha ’-1 \right) max \left( \beta ’x’,y’ \right) }{ \alpha ’- \beta ’} \} }\);

the dependence function

\(\mathrm{ A \left( u \vert \alpha ’, \beta ’ \right) =\frac{ \left( 1- \beta ’ \right) max \left( \alpha ’ \left( 1-u \right) ,u \right) + \left( \alpha ’-1 \right) max \left( \beta ’ \left( 1-u \right) ,u \right) }{ \alpha ’- \beta ’} }\), i.e.,

\(\mathrm{ A \left( u \vert \alpha ’, \beta ’ \right) =1-u~if~0 \leq u \leq \frac{ \beta ’}{ \beta ’+1} }\)  

\(\mathrm{ =\frac{ \alpha ’ \left( 1- \beta ’ \right) + \left( \alpha ’ \beta ’-1 \right) u}{ \alpha ’- \beta ’} }\)   if  \(\mathrm{\frac{ \beta ’}{ \beta ’+1} \leq u \leq \frac{ \alpha ’}{ \alpha ’+1} }\)

\(\mathrm{ =u~for\frac{ \alpha ’}{ \alpha ’+1} \leq u \leq 1 }\).

\(\mathrm{ B \left( u \vert \alpha ’, \beta ’ \right) }\) is equal to \(\mathrm{ 0~if~u< \beta ’/ \left( 1+ \beta ’ \right) }\), to \(\mathrm{ \frac{ \left( \alpha ’-1 \right) u}{ \left( 1- \beta ’ \right) \alpha ’- \left( 1- \alpha ’ \beta ’ \right) u} }\) if  \(\mathrm{ \beta ’/ \left( 1+ \beta ’ \right) \leq u< \alpha ’/ \left( 1+ \alpha ’ \right) }\)and to \(\mathrm{ 1 }\) when \(\mathrm{ \alpha ’/ \left( 1+ \alpha ’ \right) \leq u \leq 1 }\); the jumps at \(\mathrm{ u= \beta ’/ \left( 1+ \beta ’ \right) ~and~u= \alpha ’/ \left( 1+ \alpha ’ \right) }\)are respectively \(\mathrm{ \left( \alpha ’-1 \right) \beta ’/ \left( \alpha ’- \beta ’ \right) ~and~ \left( 1- \beta ’ \right) / \left( \alpha ’- \beta ’ \right) }\). The support for \(\mathrm{ \frac{Y’~}{X’+Y’} }\) is thus \(\mathrm{[ \frac{ \beta ’}{ \beta ’+1},\frac{ \alpha ’}{ \alpha ’+1} ] }\) or equivalently the angle \(\mathrm{ \beta ’X’ \leq Y’ \leq \alpha ’X’ }\) and not all the first quadrant as in the other examples.

Let us make a quick reference to a generalization of the natural model, the \(\left( 2,N \right) \) natural model, a particular case of the \(\left( m,N \right) \) natural model (for \(\mathrm{ m=2 }\) ) described briefly in the next Chapter, still under study.

Let \(\mathrm{ Z_{1},\dots,Z_{n} }\) be \(\mathrm{ N }\) independent reduced Gumbel random variables and define \(\mathrm{ X= \begin{array}{c} \mathrm{n} \\ \mathrm{max} \\ \mathrm{1} \end{array} \left( Z_{j}-p_{j} \right) ,Y= \begin{array}{c} \mathrm{n} \\ \mathrm{max} \\ \mathrm{1} \end{array} \left( Z_{j}-p_{j} \right) }\).  

The random pair \(\mathrm{ \left( X,Y \right) }\) has reduced Gumbel margins iff \(\mathrm{ \sum _{1}^{N}e^{-p_{j}}=1 }\) and \(\mathrm{ \sum _{1}^{N}e^{-q_{j}}=1 }\) then we have the distribution function

\(\mathrm{ \Lambda \left( x,y \vert p,q \right) =Prob~\{ X \leq x,Y \leq y \} =exp⁡ \{ \sum _{j=1}^{N}max⁡ \left( e^{-x-p_{j}},e^{-y-q_{j}} \right) \} }\),

with the dependence function

\(\mathrm{ k \left( w \vert p,q \right) =\frac{ \sum _{1}^{N}max \left( e^{-p_{j}},e^{-q_{j^{-w}}} \right) }{1+e^{-w}} }\).

The model thus has \(\mathrm{ 2(N-1) }\) parameters. The corresponding model for bivariate minima with standard exponential margins has the survival function

\(\mathrm{ S ( x’,y’ \vert p’,q’ ) =exp⁡ \{ - \sum _{1}^{N}max( x’p_{j}^{’},y’q_{j}^{’} ) \} }\), where \(\mathrm{p_{j}^{’}=e^{-p_{j}},q_{j}^{’}=e^{-q_{j}}, \sum _{1}^{N}p_{j}^{’}=1, \sum _{1}^{N}q_{j}^{’}=1 }\), with the dependence function

\(\mathrm{ A \left( u \vert p,q \right) = \sum _{1}^{N}max ( p_{j}^{’} \left( 1-u \right) ,q_{j}^{’}u ) }\).

4 . Intrinsic estimation of the dependence functions

We have given the necessary and sufficient conditions on \(\mathrm{ k \left( w \right) ~and~A \left( u \right) for~ \Lambda \left( x,y \right) =exp \{ - \left( e^{-x}+e^{-y} \right) k \left( y-x \right) \} }\)to be a distribution function of bivariate maxima with reduced Gumbel margins and for \(\mathrm{ S \left( x’,y’ \right) =exp \{ - \left( x’+y’ \right) A ( \frac{y’}{x’+y’} ) \} ,x’,y’ \geq 0 }\) to be the survival function of bivariate minima with standard exponential margins. Let us recall that, for \(\mathrm{ A \left( u \right) ,u\in \left[ 0,1 \right] }\), we have \(\mathrm{ A \left( 0 \right) =A \left( 1 \right) =1,0 \leq -A’ \left( 0 \right) ,A’ \left( 1 \right) \leq 1~and~A \left( u \right) }\) convex.

We will call an intrinsic estimator, sometimes also called non-parametric, of \(\mathrm{ A \left( u \right) \left( k \left( w \right) \right) }\) an estimator (i.e., a statistical function) converging in probability to the function \(\mathrm{ A \left( u \right) \left( k \left( w \right) \right) }\), for every \(\mathrm{ u \in \left[ 0,1 \right] \left( - \infty \leq w \leq + \infty \right) }\) that satisfies the necessary and sufficient conditions to be an \(\mathrm{ A \left( u \right) \left( k \left( w \right) \right) }\). In Tiago de Oliveira (1989b) and in Deheuvels and Tiago de Oliveira (1989) a first, and heavy, approach was made. Later, in Tiago de Oliveira (1991), a new and simpler approach was obtained. The procedure presented in Tiago de Oliveira (1989b), whose study was continued in Deheuvels and Tiago de Oliveira (1989), is more complicated than this one; also the estimator given there does not have all moments for any size sample, as opposed to what happens with the new one.

Let us consider, for simplicity, the intrinsic estimation of \(\mathrm{ A \left( u \right)}\) from an i.i.d. sample of \(\mathrm{ n }\) pairs \(\mathrm{ \{ \left( x_{i},y_{i} \right) \} }\) with survival function \(\mathrm{ S \left( x,y \right) }\), satisfying all necessary and sufficient conditions given above; the conversion to the estimation of \(\mathrm{ k \left( w \right) }\) is immediate.

For convenience of the exposé we will use \(\mathrm{ \Omega \left( u \right) =1-A \left( u \right) }\); we have the necessary and sufficient conditions: \(\mathrm{ \Omega \left( 0 \right) = \Omega \left( 1 \right) =0,0 \leq \Omega ^{’} \left( 0 \right) ,- \Omega ^{’} \left( 1 \right) \leq 1\,and~ \Omega \left( u \right) }\)concave and also from  \(\mathrm{ max \left( u,1-u \right) \leq A \left( u \right) \leq 1 }\) we obtain \(\mathrm{ 0 \leq \Omega \left( u \right) \leq min \left( u,1-u \right) =\mathit{l} \left( u \right) }\).

Let us consider to estimate \(\mathrm{ \Omega \left( u \right) }\), the statistical function

\(\mathrm{ \bar{\Omega} _{n} \left( u \right) = \bar{\Omega} _{n} ( u \vert x_{i}^{’},y_{i}^{’}; \alpha ) =\frac{c \left( \alpha \right) }{n} \sum _{1}^{n}min( \frac{1-u}{ \alpha +x_{i}^{’}},\frac{u}{ \alpha +y_{i}^{’}} ) }\).

where we will take \(\mathrm{ \alpha \left( = \alpha _{n} \right) \rightarrow 0 }\).

As \(\mathrm{ min ( \frac{1-u}{ \alpha +x_{i}^{’}},\frac{u}{ \alpha +y_{i}^{’}} ) \leq \frac{\mathit{l }\,( u ) }{ \alpha } }\), we see that each summand has all the moments. We will seek function \(\mathrm{ \left( s \right) c=c \left( \alpha \right) }\) and \(\mathrm{ \alpha = \alpha _{n} }\) sufficient to guarantee that \(\mathrm{ \bar{\Omega} _{n} \left( u \right) {\stackrel{m.s.}\rightarrow }\Omega \left( u \right) }\). From  the properties of each summand we can prove that \(\mathrm{ \bar{\Omega} _{n} \left( 0 \right) = \bar{\Omega} _{n} \left( 1 \right) =0 }\) and \(\mathrm{ \bar{\Omega} _{n} \left( u \right) }\) is concave, but the two other conditions on \(\mathrm{ {\Omega'} \left( u \right) }\) are not necessarily verified (the conditions on the derivatives are valid for each summand only if \(\mathrm{ 0<c \left( \alpha \right) \leq \alpha }\), which will not be the case). As a consequence we will use the truncation \(\mathrm{ \Omega _{n}^{*} \left( u \right) =min \left( \bar{ \Omega} _{n} \left( u \right) ,\mathit{l} \left( u \right) \right) }\), which is a simple intrinsic estimator of \(\mathrm{ {\Omega} \left( u \right) }\).

Let us begin by considering \(\mathrm{ \bar{\Omega} _{n} \left( u \right) }\) in the complete dependence case:

In the diagonal case, where \(\mathrm{ X’=Y’ }\) with probability one, \(\mathrm{ X’ \left( =Y’ \right) }\) has the standard exponential distribution.

Then we have

\(\mathrm{ M ( \bar{ \Omega} _{n} ( u ) ) =M ( c ( \alpha ) min⁡ ( \frac{1-u}{ \alpha +x’},\frac{u}{ \alpha +y’}) ) = ( c ( \alpha ) \int _{0}^{+ \infty}\frac{e^{-x’}}{ \alpha +x’}d~x ) \mathit{l} ( u ) = }\)

\(\mathrm{ = ( c \left( \alpha \right) e^{ \alpha } \int _{ \alpha }^{+ \infty}\frac{e^{-t}}{t}d~t ) \cdot \mathit{l} \left( u \right) =c \left( \alpha \right) J_{1} \left( \alpha \right) \mathit{l} \left( u \right) = }\)

and

\(\mathrm{ V( \bar{\Omega} _{n} ( u ) ) =\frac{1}{n}\,V ( c ( \alpha ) min⁡ ( \frac{1-u}{ \alpha +x’},\frac{u}{ \alpha +y’} ) ) =\frac{c^{2} \left( \alpha \right) }{n} [ J_{2} \left( \alpha \right) -J_{1}^{2} \left( \alpha \right) ] \mathit{l}^{2}⁡ \left( u \right) }\)

where \(\mathrm{ J_{1} \left( \alpha \right) =e^{ \alpha } \int _{ \alpha }^{+ \infty}\frac{e^{-t}}{t}d~t~and~J_{2} \left( \alpha \right) =e^{ \alpha } \int _{ \alpha }^{+ \infty}\frac{e^{-t}}{t^{2}}d~t }\); note that \(\mathrm{ J_{1}^{’} \left( \alpha \right) =J_{1} \left( \alpha \right) -1/ \alpha }\) and \(\mathrm{ J_{2}^{’} \left( \alpha \right) =J_{2} \left( \alpha \right) -1/ \alpha ^{2} }\). It is immediate that  \(\mathrm{ J_{1} \left( \alpha \right) \sim -log \,\alpha }\) and \(\mathrm{ J_{2} \left( \alpha \right) \sim -1/ \alpha }\) as \(\mathrm{\alpha \rightarrow 0 }\).

Clearly the unbiased estimator of \(\mathrm{ \mathit{l} \left( u \right) }\) leads to the choice \(\mathrm{ c \left( \alpha \right) ~J_{1} \left( \alpha \right) =1 }\) and the condition that \(\mathrm{ V \left( \bar{\Omega} _{n} \left( u \right) \right) \rightarrow 0 }\) takes the form \(\mathrm{ \frac{1}{n}( \frac{J_{2} \left( \alpha \right) }{J_{1}^{2} \left( \alpha \right) }-1 ) \rightarrow 0~or,as~\frac{J_{1}^{2} \left( \alpha \right) }{J_{2} \left( \alpha \right) } \rightarrow 0~if~ \alpha \rightarrow 0,\frac{1}{n}\frac{J_{2} \left( \alpha \right) }{J_{1}^{2} \left( \alpha \right) } \sim \frac{1}{n~ \alpha \left( log \,\alpha \right) ^{2}} \rightarrow 0~as~ \alpha = \alpha _{n} \rightarrow 0}\).

Let us analyse other choices of \(\mathrm{ c \left( \alpha \right) }\) which may be simpler and more convenient. Evidently we must have \(\mathrm{ c \left( \alpha \right) ~J_{1} \left( \alpha \right) \rightarrow 1 }\)  and \(\mathrm{ \frac{c^{2} \left( \alpha \right) }{n} ( J_{2} \left( \alpha \right) -J_{1}^{2} \left( \alpha \right) ) \rightarrow 0~or~c \left( \alpha \right) J_{1} \left( \alpha \right) \rightarrow 1 }\)and, also, \(\mathrm{ \frac{1}{n~ \alpha \left( log\, \alpha \right) ^{2}} \rightarrow 0 }\). The mean-square error of the estimator is thus

\(\mathrm{ MSE= ( c \left( \alpha \right) J_{1} \left( \alpha \right) -1 ) ) ^{2}\mathit{l}^{2}⁡ \left( u \right) +\frac{c^{2} \left( \alpha \right) }{n} ( J_{2} \left( \alpha \right) -J_{1}^{2} \left( \alpha \right) ) \mathit{l}^{2}⁡ \left( u \right) }\),

which converges to zero under the conditions stated above.

We could seek the function \(\mathrm{ c \left( \alpha \right) }\) minimizing the MSE but this has no practical interest.

As \(\mathrm{ \int _{0}^{ \alpha }\frac{1-e^{-t}}{t}d~t- \int _{ \alpha }^{+ \infty}\frac{e^{-t}}{t}d~t= \gamma +log~ \alpha }\) we see immediately that \(\mathrm{ J_{1} \left( \alpha \right) +log~ \alpha \rightarrow - \gamma ~as~ \alpha \rightarrow 0 }\). But the choice \(\mathrm{ c \left( \alpha \right) =1/ \left( - \gamma -log~ \alpha \right) c \left( \alpha \right) J_{1} \left( \alpha \right) \rightarrow 1 )~ has~c \left( \alpha \right) \leq 0 }\) for values of \(\mathrm{ \alpha \approx 1 }\). Let us take then \(\mathrm{ c \left( \alpha \right)} =\frac{1}{a-{\mathrm{log~ \alpha}}}\). The bias is then \(\mathrm{ \left( c \left( \alpha \right) J_{1} \left( \alpha \right) -1 \right) ~\mathit{l} \left( u \right)} = ( \frac{\mathrm{J_{1}} \left( \alpha \right) }{a-{\mathrm{log~ \alpha }}}-1 ) ~\mathit{l} \left({\mathrm{ u}} \right) \sim \frac{ \gamma +a}{\mathrm{log~ \alpha }}~\mathit{l} \left({\mathrm{ u}} \right) \); a must be chosen such that \({ a-\mathrm{log~ \alpha >0 }}\) and as we are going to deal with small values of \(\mathrm{ \alpha }\), we can suppose \(\mathrm{ \alpha \leq 1 }\) and then take \({ a>0 }\); note that this is equivalent to using \(\mathrm{ \alpha ’=e^{-a}~ \alpha }\) instead of \(\mathrm{ \alpha }\).

We can then compute asymptotically the mean-square error of the estimator. We have

\(\mathrm{ MSE}= \{ [ ( \frac{a+ \gamma }{a-{\mathrm{log~ \alpha }}} ) ^{2}+0 (\frac{1}{ \left( a-{\mathrm{log~ \alpha }}\right) ^{2}} ) ] + [ \frac{1}{\mathrm{n}~ \alpha \left( a-{\mathrm{log⁡ \,\alpha}} \right) ^{2}}+0 ( \frac{1}{\mathrm{n}~ \alpha \left( a-{\mathrm{log}}⁡ \,\alpha \right) ^{2}} ) ] \} \mathit{l}^{2}⁡ \left( {\mathrm{u}} \right) \)

where the first term is the square of the bias and the second corresponds to the asymptotic evaluation of the variance  \({ \frac{\mathrm{J}_{2} \left( \alpha \right) -{\mathrm{J}}_{1}^{2} \left( \alpha \right) }{\mathrm{n}~ \left( a-{\mathrm{log}}\,⁡ \alpha \right) ^{2}} }\)  because \(\mathrm{ J_{2} \left( \alpha \right) \sim 1/ \alpha }\) and \(\mathrm{ J_{1} \left( \alpha \right) \sim -log~ \alpha }\) , as stated previously.

As \(\mathrm{ \alpha \rightarrow 0 }\), and so \(\mathrm{ -log~ \alpha + \infty }\), the condition on \(\mathrm{ \alpha = \alpha _{n} }\)  is to be such that \({ \mathrm{n}~ \alpha \left( a-{\mathrm{log}} \,\alpha \right) ^{2} }\)  or \({ \mathrm{n}~ \alpha \left( a-{\mathrm{log}} \,\alpha \right) ^{2}\rightarrow + \infty }\).    

The discussion of the dominant terms in the MSE shows immediately that the squared bias is dominant if \(\mathrm{ n~ \alpha _{n} \rightarrow \infty }\), is of the same order the variance if \(\mathrm{ n~ \alpha _{n} \rightarrow H \left( 0<H<+ \infty \right) }\), and the variance term is dominant if  \(\mathrm{ n~ \alpha _{n} \rightarrow 0 }\).

The simplest form for \(\mathrm{ \alpha = \alpha _{n} }\) is to take \(\mathrm{ \alpha _{n}=Rn^{-p} \left( p>0 \right) }\). The condition \(\mathrm{ n~ \alpha \left( log \,\alpha \right) ^{2} \rightarrow + \infty }\) immediately imposes \(\mathrm{ 0<p \leq 1 }\) and the square of the bias and variance are asymptotic to, respectively, \({ \frac{ \left( a+ \gamma \right) ^{2}}{\mathrm{p}^{2} \left( {\mathrm{log~n}} \right) ^{2}}} \) and \(\mathrm{ \frac{1}{R~p^{2}~n^{1-p} \left( log~n \right) ^{2}} }\). Their ratio is asymptotic to \({ \left( a+ \gamma \right) ^{2}~{\mathrm{R~n^{1-p} }}}\) having the squared bias (asymptotic) dominance if \(\mathrm{ p<1 }\) and (asymptotic) equivalence if \(\mathrm{ p=1 }\) . In both cases the MSE is of the order of \(\mathrm{ \left( log~n \right) ^{-2} }\) , and it is simpler to take \(\mathrm{ p=1 }\) and so  \(\mathrm{ \alpha _{n}=R/n }\). Then \(\mathrm{ c \left( \alpha _{n} \right)} =\frac{1}{a-{\mathrm{log\,R-log\,n}}} \) and the condition \(\mathrm{ c \left( \alpha _{n} \right) >0 }\) implies \({ a-{\mathrm{log~R>0 }}}\). We can take \(\mathrm{ R=1 \left( \alpha _{1}=1 \right) }\) and dispose of \({ a}\). Taking for simplicity \(a=1\)  we get the simple expression \(\mathrm{ \alpha _{n}=1/n }\) and \(\mathrm{ c \left( \alpha _{n} \right) =\frac{1}{1-log \,\alpha _{n}}=\frac{1}{1+log\,n} }\), and then the squared bias is \(\mathrm{ \frac{ \left( 1+ \gamma \right) ^{2}}{ \left( 1+log~n \right) ^{2}}+0 ( \frac{1}{ \left( 1+log⁡\,n \right) ^{2}} ) }\), the variance is \(\mathrm{ \frac{1}{ \left( 1+log⁡\,n \right) ^{2}}+0 ( \frac{1}{ \left( 1+log⁡\,n \right) ^{2}} })\) and the MSE is \(\mathrm{ \frac{1+ \left( 1+ \gamma \right) ^{2}}{ \left( 1+log~n \right) ^{2}}+0 ( \frac{1}{ \left( 1+log\,⁡n \right) ^{2}} ) }\). The convergence, being of the order of \(\mathrm{ 1/ \left( log⁡\,n \right) ^{2} }\), is thus very slow. Notice that the condition \(\mathrm{ c \left( \alpha \right) < \alpha }\)  which was needed for the properties of \(\mathrm{ \Omega ’ \left( 0 \right) }\) and \(\mathrm{ \Omega ’ \left( 1 \right) }\) is not verified.

Clearly as \(\mathrm{ \bar{\Omega} _{n} \left( u \right)\stackrel {m.s.}\rightarrow \Omega \left( u \right) =\mathit{l} \left( u \right) ~then~ \bar{\Omega} _{n} \left( u \right) \stackrel{P}\rightarrow\,\mathit{l} \left( u \right) }\) and so  \(\mathrm{ \Omega _{n}^{*} \left( u \right) =min⁡ \left( \bar{\Omega} _{n} \left( u \right) ,\mathit{l} \left( u \right) \right) \stackrel{P}\rightarrow\mathit{l} \left( u \right) }\).

As a further note, consider two choices: \(\mathrm{ \alpha _{n}=R~n^{-p}~and~ \alpha _{n}^{’}=R’n^{-p’} ( 0<p,p’ \leq 1 ) }\). The ratio of the variances corresponding to \(\mathrm{ \alpha _{n} }\) and \(\mathrm{ \alpha _{n}^{'} }\) is

\(\mathrm{\sim \frac{{1}/{R~p^{2}n^{1-p} \left( log\,n \right) ^{2}}}{{1}/{R’}p’^{2}\,n^{1-p^{’}} \left( log\,n \right) ^{2}}=\frac{R^{’}p’^{2}}{R~p^{2}}n^{p-p^{’}} }\),

which converges to zero (variance of smaller asymptotic order) if  \(\mathrm{ n^{p-p^{’}} \rightarrow 0~or~p<p^{’} }\), showing that although the variance would be of smaller order with smaller \(\mathrm{ p }\) this effect is outweighed by the bias effect and so the choice of \(\mathrm{ \alpha _{n}=\frac{1}{1+log\,n} }\) is reasonable.

Consider now the general case

\(\mathrm{ \bar{\Omega} _{n} \left( u \right) =\frac{c \left( \alpha \right) }{n} \sum _{1}^{n}min ( \frac{1-u}{ \alpha +x_{i}^{’}},\frac{u}{ \alpha +y_{i}^{’}} ) }\)

where, as seen before, we will take \(\mathrm{ \alpha _{n}=1/n }\) and \(\mathrm{ c \left( \alpha _{n} \right) =\frac{1}{1-log~ \alpha _{n}}=\frac{1}{1+log~n} }\).

We have

\(\mathrm{ Q_{ \alpha ^{ \left( z \right) }}=Prob \{ c \left( \alpha \right) ~min ( \frac{1-u}{ \alpha +X^{’}},\frac{u}{ \alpha +Y^{’}} ) \geq z \} = }\)

\(\mathrm{ =Prob \{ X^{’} \leq c \left( \alpha \right) \frac{ \left( 1-u \right) }{z}- \alpha ,Y’ \leq c \left( \alpha \right) \frac{u}{z}- \alpha \} }\)

\(\mathrm{ =1-S ( \frac{c \left( 1-u \right) }{z}- \alpha ,0 ) -S ( 0,\frac{c,u}{z}- \alpha ) +S ( \frac{c \left( 1-u \right) }{z}- \alpha ,\frac{c \cdot u}{z}- \alpha) }\).

Then the mean value of   \(\mathrm{ c \left( \alpha \right) min ( \frac{1-u}{ \alpha +X^{’}},\frac{u}{ \alpha +Y^{’}} ) }\) is given by \(\mathrm{ \int _{0}^{+ \infty}z~d \left( 1-Q_{ \alpha } \left( z \right) \right) }\)

\(\mathrm{ = \int _{0}^{+ \infty}z~d~S ( \frac{c \left( 1-u \right) }{z}- \alpha ,0)+ \int _{0}^{+ \infty}z~d~S ( \frac{c~u}{z}- \alpha ,0 ) - \int _{0}^{+ \infty}z~d~S ( \frac{c \left( 1-u \right) }{z}- \alpha ,\frac{c~u}{z}- \alpha ) }\)

\(\mathrm{ =I_{1} \left( \alpha \right) +I_{2} \left( \alpha \right) -I_{3} \left( \alpha \right) }\).

Let us then compute the summands. We have

\(\mathrm{ I_{1} \left( \alpha \right) = \int _{0}^{+ \infty}z~d~S( \frac{c \left( 1-u \right) }{z}- \alpha ,0 ) = \int _{0}^{+ \infty}~z~d ( e^{- \left( \left( {c \left( 1-u \right) }/{z} \right) - \alpha \right) _{+}} ) = }\)

\(\mathrm{ \int _{ \left( {c \left( 1-u \right) }/{z} \right) - \alpha >0}^{+ \infty}z~d \left( e^{- \left( \left( {c \left( 1-u \right) }/{z} \right) - \alpha \right) } \right) =e^{ \alpha }c \left( 1-u \right) \int _{ \alpha }^{+ \infty}\frac{e^{- \xi }}{ \xi }d~ \xi }\)

and

\(\mathrm{ I_{2} \left( \alpha \right) = \int _{0}^{+ \infty}z~d~S ( 0,\frac{c~u}{z}- \alpha ) = \int _{0}^{+ \infty}z~d ( e^{- ( ( ( {c~u ) }/{z} ) - \alpha ) _{+}} ) = }\)

\(\mathrm{ \int _{ \left( ( {c~u ) }/{z} \right) - \alpha >0}^{+ \infty}~z~d ( e^{- ( ( ( {c ( 1-u ) }/{z} ) - \alpha ) } ) =e^{ \alpha }c~u \int _{ \alpha }^{+ \infty}\frac{e^{- \xi }}{ \xi }e~ \xi }\)

and so \(\mathrm{ I_{1} \left( \alpha \right) +I_{2} \left( \alpha \right) =e^{ \alpha }~c \left( \alpha \right) J_{1} \left( \alpha \right) \rightarrow 1 }\).

Now compute \(\mathrm{ I_{3} \left( \alpha \right) }\). We will show that \(\mathrm{ I_{3} \left( \alpha \right) \rightarrow A \left( u \right) }\) and so \(\mathrm{ M \left( \bar{\Omega} _{n} \left( u \right) \right) =I_{1} \left( \alpha \right) +I_{2} \left( \alpha \right) -I_{3} \left( \alpha \right) \rightarrow 1-A \left( u \right) = \Omega \left( u \right) }\).

We have

\(\mathrm{ I_{3} \left( \alpha \right) = \int _{0}^{+ \infty}z~d [ exp⁡ \{ - [ ( \frac{c \left( 1-u \right) }{z}- \alpha ) _{+}+ ( \frac{c~u}{z}- \alpha ) _{+} ] A ( \frac{ ( \frac{c~u}{z}- \alpha ) _{+}}{ ( \frac{c ( 1-u ) }{z}- \alpha ) _{+}+ \left( \frac{c~u}{z}- \alpha \right) _{+}} ) \}] }\).

Let us suppose, for convenience, that \(\mathrm{ 0 \leq u \leq 1-u \leq 1 }\). Then \(\mathrm{ ( \frac{c \left( 1-u \right) }{z}- \alpha ) _{+} }\) is different from zero when \(\mathrm{ c \left( 1-u \right) > \alpha ~z }\)  and \(\mathrm{ \left( \frac{c~u}{z}- \alpha \right) _{+} }\) is different from zero when \(\mathrm{ c~u> \alpha ~z }\). Then  \(\mathrm{ I_{3} \left( \alpha \right) }\) can be written as

\(\mathrm{ I_{3} \left( \alpha \right) = \int _{0}^{cu/ \alpha }z~d [ e^{- \left( {c}/{z}-2~ \alpha \right) A \left( \left( {cu}/{z} \right) - \alpha \right) / \left( \left( {c}/{z} \right) -2~ \alpha \right) ) } ] + \int _{cu/ \alpha }^{c \left( 1-u \right) / \alpha }z~d ( e^{- \left( \left( {c \left( 1-u \right) }/{z} \right) - \alpha \right) } ) }\) .

The second integral is

\(\mathrm{ c(\alpha)e^\alpha (1-u) \int _{ \alpha }^{ ( ( {1-u ) }/{u} ) \alpha }\frac{e^{- \xi }}{ \xi }d~ \xi =\frac{e^{ \alpha }}{1-log~ \alpha } ( 1-u ) \int _{ \alpha }^{ ( ( {1-u ) }/{u} ) \alpha }\frac{e^{- \xi }}{ \xi }d~ \xi \rightarrow 0~as~ \alpha \rightarrow 0 } \).

Let us consider the first integral. It can be written as

\(\mathrm{ - \int _{0}^{c~u/ \alpha }z~d( 1-e^{- \left( {c}/{z}-2~ \alpha \right) A \left( \left( c{u}/{z}- \alpha \right) / \left( {c}/{z}-2~ \alpha \right) \right) } ) }\),

which, by integration by parts, takes the form

\(\mathrm{ - [ z ( 1-e^{- \left( {c}/{z}-2~ \alpha \right)\, A\, \left( \left( c\,{u}/{z}- \alpha \right) / \left( {c}/{z}-2~ \alpha \right) \right) } ) \vert \begin{array}{c} \mathrm{c\,u/\alpha} \\ \mathrm{0} \end{array} ] + }\)

\(\mathrm{ + \int _{0}^{c\,u/ \alpha } ( 1-e^{- ( {c}/{z}-2\, \alpha ) \,A\, (( c\,{u}/{z}- \alpha ) / ( {c}/{z}-2\, \alpha )) } ) d\,z }\)

\(\mathrm{ =-\frac{c ( \alpha ) u}{ \alpha } ( 1-e^{- ( {1}/{u}-2) \alpha } ) + \int _{0}^{c \,u/ \alpha } ( 1-e^{- ( {c}/{z}-2\, \alpha ) A( ( c\, {u}/{z}- \alpha) / ( {c}/{z}-2\, \alpha)) } ) d\,z }\).

The first summand, as \(\mathrm{ c \left( \alpha \right) =\frac{1}{1-log~ \alpha } }\), converges to zero. Let us, then, study the integral. Now put \(\mathrm{ \xi =c/z }\) .The integral takes the form

\(\mathrm{ = \int _{ \alpha /u}^{+ \infty}c \left( \alpha \right) ( 1-e^{- \left( \xi -2\, \alpha \right) \,A\, (( { \xi u- \alpha) }/{ (\xi }-2\, \alpha ) ) }\,\frac{d\, \xi }{ \xi ^{2}} }\).

But as \(\mathrm{ u \leq 1-u \left( u \leq 1/2 \right) }\) we have \(\mathrm{ \frac{ \xi \, u- \alpha }{ \xi -2 \,\alpha }<u }\) and as \(\mathrm{ A \left( u \right) /u }\) is non-increasing we obtain

\(\mathrm{ \frac{A ( \frac{ \xi\, u- \alpha }{ \xi -2\, \alpha } ) }{\frac{ \xi\, u- \alpha }{ \xi -2\, \alpha }} \geq \frac{A \left( u \right) }{u}~or~ ( \xi - { \alpha }/{u} ) \,A\, \left( u \right) \leq \left( \xi -2\, \alpha \right) \,A\, ( \frac{ \xi\, u- \alpha }{ \xi -2\, \alpha } ) }\);

and from \(\mathrm{ \frac{A \left( u \right) }{1-u} }\)non-decreasing we get \(\mathrm{ \left( \xi -2~ \alpha \right) ~A ( \frac{ \xi u- \alpha }{ \xi -2~ \alpha } ) \leq ( \xi -\frac{ \alpha }{1-u} ) A \left( u \right) }\)and, thus, the integral is bounded above and below by

\(\mathrm{ \int _{ \alpha /u}^{+ \infty}c \left( \alpha \right) ( 1-e^{- \left( \xi -\alpha/ \left( { 1-u} \right) \right) ~A \left( u \right) } ) \frac{d~ \xi }{ \xi ^{2}}~and \int _{ \alpha /u}^{+ \infty}c \left( \alpha \right) ( 1-e^{- \left( \xi - \left( { \alpha }/{u} \right) \right) ~A \left( u \right) } ) \frac{d~ \xi }{ \xi ^{2}} }\).

Putting \(\mathrm{ \xi ~A \left( u \right) =t }\) we get, respectively, the following bounds above and below:

\(\mathrm{ A \left( u \right) c \left( \alpha \right) \int _{ \alpha\, A/u}^{+ \infty} \left( 1-e^{ \alpha \,A/ \left( 1-u \right) -t} \right) \frac{d~t}{t^{2}}~and~A \left( u \right) c \left( \alpha \right) \int _{ \alpha \,A/u}^{+ \infty} \left( 1-e^{ \left( \alpha \,A/u \right) -t} \right) \frac{d~t}{t^{2}} }\).

But \(\mathrm{ \int _{p \alpha }^{+ \infty} ( 1-e^{q \alpha -t} ) \frac{d~t}{t^{2}} }\) , with \(\mathrm{ p,q>0 }\), by integration by parts, has the expression \(\mathrm{ \frac{1-e^{ \left( q-p \right) \alpha }}{p^{ \alpha }}+e^{q \alpha } \int _{p \alpha }^{+ \infty}\frac{e^{-t}}{t}d~t }\) and with \(\mathrm{ c \left( \alpha \right) =\frac{1}{1-log~ \alpha } }\) we see immediately that \(\mathrm{ c \left( \alpha \right) \int _{p \alpha }^{+ \infty} \left( 1-e^{q \alpha -t} \right) \frac{d~t}{t^{2}} \rightarrow 1 }\), and so both bounds converge to \(\mathrm{ A \left( u \right) }\). Thus, as said initially, \(\mathrm{ I_{3} \left( \alpha \right) \rightarrow A \left( u \right) }\), and so the mean value \(\mathrm{ I_{1} \left( \alpha \right) +I_{2} \left( \alpha \right) -I_{3} \left( \alpha \right) \rightarrow 1-A \left( u \right) = \Omega \left( u \right) }\).  The proof when \(\mathrm{ 1-u \leq u }\) leads to the same result.

By the same technique we can show that

\(\mathrm{ V( c ( \alpha) min ( \frac{1-u}{ \alpha +X^{’}},\frac{u}{ \alpha +Y^{’}} ) ) \sim \frac{c^{2} ( \alpha ) }{ \alpha }C ( u ) }\),

when \(\mathrm{ C \left( u \right) }\) is a bounded function of \(\mathrm{ u }\) .

The condition of convergence is, once more, \(\mathrm{ V ( \bar{\Omega} _{n} \left( u \right) ) \sim \frac{C \left( u \right) }{n~ \alpha \left( log~ \alpha \right) ^{2}} \rightarrow 0~or~n~ \alpha \left( log~ \alpha \right) ^{2} \rightarrow 0 }\); as before for the diagonal case the previous analysis applies.

The final estimator is, as before, \(\mathrm{ \Omega _{n}^{*} \left( u \right) =min \left( \bar{\Omega} _{n} \left( u \right) ,\mathit{l} \left( u \right) \right) \stackrel{P}\rightarrow \Omega \left( u \right) =1-A \left( u \right) }\).

We have, then, for bivariate minima with standard exponential margins, as \(\mathrm{ \alpha _{n}=1/n }\) ,

\(\mathrm{ \bar{A}_{n} \left( u \right) =1-\frac{1}{1+log~n} \sum _{1}^{n}min⁡ ( \frac{1-u}{1+n~x_{i}^{’}},\frac{u}{1+n~y_{i}^{’}} ) }\),

and, more generally,

\(\mathrm{ A_{n}^{*} \left( u \right) =1-min⁡\, ( \mathit {l} ( u ) ,\frac{1}{1+log~n} \sum _{1}^{n}\,min\,⁡ ( \frac{1-u}{1+n~x_{i}^{’}},\frac{u}{1+n~y_{i}^{’}} )) }\),

Recall that \(\mathrm{ \bar{A}_{n} \left( u \right) }\) can be used if \(\mathrm{ 0 \leq -\bar{A}_{n}^{’} \left( 0 \right) ,\bar{A}_{n}^{’} \left( 1 \right) \leq 1 }\) that is iff

\(\mathrm{ 0 \leq \sum _{1}^{n}\frac{1}{1+n~x_{i}^{’}}, \sum _{1}^{n}\frac{u}{1+n~y_{i}^{’}} \leq 1+log~n }\),

and, in that case, we get

\(\mathrm{ \bar{S}_{n}( x,y ) =e^{- ( x+y ) \bar{A}_{n} ( y/ ( x+y ) ) }=e^{- ( x+y ) +1( 1+log~n) } \cdot \sum _{1}^{n}min⁡ ( x/ ( 1+n~x_{i}^{’} ) ,y/( 1+n~y_{i}^{’} ) ) }\).

As \(\mathrm{ \bar{A}_{n} \left( u \right)\stackrel {P}\rightarrow A \left( u \right) }\) and \(\mathrm{ A_{n}^{*} \left( u \right) =max \left( \bar{A}_{n} \left( u \right) ,1-\mathit{l} \left( u \right) \right) \stackrel {P}\rightarrow A \left( u \right) }\), we see that \(\mathrm{ \bar{S}_{n}( x’,y’ ) =exp⁡ \{ - ( x’+y’) \bar{A}_{n} ( \frac{y’}{x’+y’} ) \}\stackrel{P}\rightarrow S ( x’,y’ ) }\) and also  \(\mathrm{ S_{n}^{*}( x’,y’ ) =exp⁡ \{ - ( x’+y’ )~A_{n}^{*} ( \frac{y’}{x’+y’} ) \}\stackrel {P} \rightarrow S ( x’,y’ ) }\).

For bivariate maxima with reduced Gumbel margins we have \(\mathrm{ \bar{k}_{n} \left( w \right) =\bar{A}_{n} ( \frac{1}{1+e^{w}}) }\)  and \(\mathrm{ k_{n}^{*} \left( w \right) =A_{n}^{*} ( \frac{1}{1+e^{w}}) }\). Denoting by \(\mathrm{ \{ ( x_{i},y_{i} ) \} }\) the sample we have (with the transformation \(\mathrm{ \left( x,y \right) ~to~ \left( x’,y’ \right) = \left( e^{-x},e^{-y} \right) }\) as stated,

\(\mathrm{ \bar{k}_{n} \left( w \right) =1-\frac{1}{ \left( 1+log~n \right) \left( 1+e^{w} \right) } \sum _{1}^{n}min⁡ ( \frac{e^{w}}{1+e^{-n~x_{i}}},\frac{1}{1+e^{-n~y_{i}}} ) }\) and

\(\mathrm{ \bar{ A}_{n} \left( x,y \right) =exp⁡ \{ - \left( e^{-x}+e^{-y} \right) +\frac{1}{1+log~n} \sum _{1}^{n}min⁡ ( \frac{e^{-x}}{1+n~e^{-x_{i}}},\frac{1}{1+n~e^{-y_{i}}} ) \} }\).

Pickands (1981) gave a much more complex estimator, converging with probability one but with a very tedious calculation.

5 . Miscellaneous statistical results:

A ready-to-wear model:

In Ramberg, Dudewicz, Tadikamalle and Mytyka (1979) a family of quantile functions \(\mathrm{ Q \left( p \right) =F^{-1} \left( p \right) }\) is proposed that the authors say fit the data well; this family, of the form (where \(\mathrm{ 0^{0}=1 }\))

\(\mathrm{ Q \left( p \right) = \lambda + \delta ( p^{ \alpha }- \left( 1-p \right) ^{ \beta } ) }\),

extends Tukey’s \(\mathrm{ \lambda }\)-functions and so contains the logistic distribution. As the latter plays a central role for \(\mathrm{ D \left( w \right) }\) functions — corresponding to independence — it seems natural to study the quantile functions \(\mathrm{ Q \left( p \right) }\) that verify the conditions to be a \(\mathrm{ D \left( w \right) }\) function. First, as \(\mathrm{ D \left( w \right) }\) has all moments \(\mathrm{ \mu _{k}^{’}= \int _{0}^{1}Q^{k} \left( p \right) d~p,k \geq 0 }\), we see immediately that we must have \(\mathrm{ \alpha \geq 0 }\) and \(\mathrm{ \beta \geq 0 }\); in the contrary case for \(\mathrm{ k }\) integer, such that \(\mathrm{ k~ \alpha +1<0 \leq \left( k-1 \right) \alpha +1~or~k~ \beta +1<0 \leq \left( k-1 \right) \beta +1 }\), the integrals \(\mathrm{ \int _{0}^{1}p^{k \alpha }d~p~or \int _{0}^{1} \left( 1-p \right) ^{k \beta }d~p }\) would diverge. Then for \(\mathrm{ Q \left( p \right) }\) to be non-decreasing we must have \(\mathrm{ \delta >0 }\) (the case \(\mathrm{ \delta =0 }\) leads to the almost sure random variable at \(\mathrm{ \lambda }\), a case that can also be obtained with \(\mathrm{ \alpha=0 }\) and \(\mathrm{ \beta=0 }\) ). It is immediate that the conditions \(\mathrm{ \int _{- \infty}^{+ \infty}w~d~D \left( w \right) =0~and~D’ \left( w \right) \geq D \left( w \right) ( 1-D \left( w \right) ) }\) are equivalent to \(\mathrm{ \int _{0}^{1}Q \left( p \right) d~p=0 }\) and \(\mathrm{ \left( 0 \leq \right) p \left( 1-p \right) Q’ \left( p \right) \leq 1 }\). The condition \(\mathrm{ \int _{0}^{1}Q \left( p \right) d~p=0 }\) gives

\(\mathrm{ \lambda + \delta ( \frac{1}{ \alpha +1}-\frac{1}{ \beta +1}) =0 }\),

and so

\(\mathrm{ Q ( p ) = \delta \{ ( p^{ \alpha }-\frac{1}{ \alpha +1} ) + (( 1-p ) ^{ \beta }-\frac{1}{ \beta +1} ) \} }\),

\(\mathrm{ Q ( 0) = \delta \{ ( 0^{ \alpha }-\frac{1}{ \beta +1} ) - ( 1-\frac{1}{ \beta +1} ) \} }\),

\(\mathrm{ Q( 1) = \delta \{( 1-\frac{1}{ \alpha +1}) -( 0^{ \beta }-\frac{1}{ \beta +1} ) \} }\).

The support is \(\mathrm{ \left[ Q \left( 0 \right) ,Q \left( 1 \right) \right] }\) and the range \(\mathrm{ R= \delta \left( 2-0^{ \alpha }-0^{ \beta } \right) ,i.e.,R=2~ \delta ~if~ \alpha >0~and~ \beta >0 }\), \(\mathrm{ R= \delta ~if~ \alpha >0~and~ \beta =0~or~ \alpha =0~and~ \beta =1~and~R=0 }\) in the degenerated case where \(\mathrm{ W \left( =Y-X \right) =0 }\) with probability \(\mathrm{ 1 }\).

Let us now consider the condition \(\mathrm{ p \left( 1-p \right) Q’ \left( p \right) \leq 1 }\) for each \(\mathrm{ p \in \left[ 0,1 \right] }\) which is equivalent to \(\mathrm{ \delta \begin{array}{c}\\ \mathrm{sup} \\ \mathrm{0\leq\,p\leq1} \end{array}[ \alpha\, ( p^{ \alpha } ( 1-p ) + \beta\, ( 1-p ) ^{ \beta } ] \leq 1 }\) relating the three parameters \(\mathrm{ \alpha \geq 0, \beta \geq 0 }\) and \(\mathrm{ \delta > 0 }\).

We will consider only a limiting case when \(\mathrm{ \delta \rightarrow + \infty }\) that already includes the logistic distribution; note that, as we have shown with the natural model, we can have cases where the range of \(\mathrm{ W=Y-X }\) is a finite one. As \(\mathrm{ \frac{p^{ \alpha }-1/ \left( \alpha +1 \right) }{ \alpha } }\) has a limit \(\mathrm{\left ( 1+log~p \right) if~ \alpha \rightarrow 0 }\), then \(\mathrm{ \delta ( p^{ \alpha }-\frac{1}{ \alpha +1}) }\) has a limit \(\mathrm{ A \left( 1+log~p \right) ~iff~ \delta\, \alpha \rightarrow A ( 0 \leq A<+ \infty) }\); analogously we shall have \(\mathrm{ \delta \beta \rightarrow B ( 0 \leq A<+ \infty ) }\) and the limiting quantile function is \(\mathrm{ A \left( 1+log~p \right) -B \left( 1+log \left( 1-p \right) \right) ,A \geq 0,B \geq 0 }\),

Note that \(\mathrm{ A= 0 }\) gives the exponential \(\mathrm{ p=1-e^{- \left( 1+{w}/{B} \right) }~if~w \geq -B~and~p=0~for~w \leq -B }\);\(\mathrm{ for~B=0~we~obtain~p=e^{ \left( {w}/{A} \right) -1}~if~w \leq A~and~p=1~if~w \geq A;for~A=B }\) we get the logistic distribution \(\mathrm{ p= ( 1+e^{-{w}/{A}} ) ^{-1} }\).

The condition \(\mathrm{ 0 \leq p \left( 1-p \right) Q’ \left( p \right) \leq 1 }\) takes the form \(\mathrm{ 0 \leq A \left( 1-p \right) +B\,p \leq 1 }\) and so we get the conditions on \(\mathrm{ A~and~B,0 \leq A,B \leq 1 }\).

The estimation of \(\mathrm{ A }\) and \(\mathrm{ B }\) can easily be made by least squares, i.e., associating the ordered statistics, \(\mathrm{ w’_{i,n} }\) to the plotting positions \(\mathrm{ p_{i,n} }\), and thus minimizing \(\mathrm{ \sum _{1}^{n} [ w’_{i,n}- A \left( 1+log~p_{i,n} \right) -B \left( 1+log \left( 1-p_{i,n} \right) \right) ] ^{2} }\)and obtaining systematic (i.e., linear on the ordered statistics) estimators. There is no special reason to use this method, which also would introduce the question of plotting positions.

Let us compute some correlation coefficients. We have

\(\mathrm{ \rho =1-\frac{3}{ \pi ^{2}} \int _{- \infty}^{+ \infty}w^{2}d~D \left( w \right) =1-\frac{3}{ \pi ^{2}} \int _{0}^{1}Q^{2} \left( p \right) d~p }\)

\(\mathrm{=1-\frac{3}{ \pi ^{2}} \{ A^{2}+B^{2}-2 ( 1-\frac{ \pi ^{2}}{6} ) A~B \} }\)

as

\(\mathrm{\int _{0}^{1}log~p~d~p= \int _{0}^{1}log \left( 1-p \right) d~p=-1, \int _{0}^{1} \left( 1+log~p \right) ^{2}d~p\mathrm{ = \int _{0}^{1} \left( 1+log⁡ \left( 1-p \right) \right) ^{2}d~p=1} }\)

and

\(\mathrm{ \int _{0}^{1}log~p~log~ \left( 1-p \right) d~p=\frac{ \partial ^{2}~B \left( a,b \right) }{ \partial ~a~ \partial ~b~} \vert _{a=b=1}=2- \pi ^{2}/6 }\)

We have then \(\mathrm{\rho ( A,B) =1-AB-\frac{3}{ \pi ^{2}} ( A-B ) ^{2}=1+ ( \frac{12}{ \pi ^{2}}-1 ) AB-\frac{3}{ \pi ^{2}} ( A+B) ^{2} }\). The logistic model \(\mathrm{ \left( A=B \right) }\) has then \(\mathrm { \rho =1-A^{2}}\) ( \(\mathrm{ \theta =1-A }\) was the parameter previously used) oscillating between \(\mathrm{ 0 \left( A=1 \right) }\) and \(\mathrm{ 1 \left( A=0 \right) }\).

Also \(\mathrm{\tau=1- \int _{-\infty}^{+ \infty}w~d~D^{2} \left( w \right) =1- \int _{0}^{1}Q \left( p \right) d~p^{2}=1-\frac{A+B}{2} }\). Then taking \(\mathrm{ \tilde{\rho} ^{*}=max \left( 0, \rho ^{*} \right) }\) and \(\mathrm{ \tilde{ \tau}^{*}=max \left( 0, \tau^{*} \right) }\) where \(\mathrm{ \rho^* }\) and  \(\mathrm{ \tau ^* }\) are the usual estimators of \(\mathrm{ \rho }\) and \(\mathrm{ \tau }\), and the truncation by zero results from the fact that  \(\mathrm{ \rho \geq 0 }\) and  \(\mathrm{ \tau \geq 0 }\), the equations \(\mathrm{ \tilde{\rho} ^{*}= \rho \left( A^{*},B^{*} \right) }\) and \(\mathrm{ \tilde{\tau}^{*}= \tau \left( A^{*},B^{*} \right) }\) do not define uniquely \(\mathrm{ A^{*} }\) and \(\mathrm{ B^{*} }\) because of the symmetry of \(\mathrm{ \rho \left( A,B \right) }\) and \(\mathrm{ \tau \left( A,B \right) }\) in \(\mathrm{ A }\) and \(\mathrm{ B }\), but give two alternatives for \(\mathrm{ \left( A^{*},B^{*} \right) }\) corresponding to the choice for \(\mathrm{ A^{*} }\) of the smallest or the largest of a second degree equation and dually for \(\mathrm{ B^{*} }\).

A test of the logistic model as \(\mathrm{ \rho =1-A^{2} }\) and  \(\mathrm{ \tau=1-A }\) is a test of  \(\mathrm{ \rho = \tau \left( 2- \tau \right) }\).

The use of quantiles can estimate \(\mathrm{ A }\) and \(\mathrm{ B }\) in an unique way.

As \(\mathrm{ A }\) and \(\mathrm{ B }\) play symmetrical roles it is natural to use quantiles of probabilities \(\mathrm{ p }\) and \(\mathrm{ q=1-p \left( >p \Leftrightarrow p \leq 1/2 \right) }\) and solve the equations

\(\mathrm{ Q^{*} \left( p \right) =A^{*} \left( 1+log~p \right) -B^{*} \left( 1+log~q \right) }\)

\(\mathrm{ Q^{*} \left( q \right) =A^{*} \left( 1+log~q \right) -B^{*} \left( 1+log~p \right) }\).

We get

\(\mathrm{ A^{*}=\frac{ \left( 1+log~q \right) ~Q^{*} \left( q \right) - \left( 1+log~p \right) Q^{*} \left( p \right) }{ \left( 1+log~q \right) ^{2}- \left( 1+log~p \right) ^{2}} }\)

\(\mathrm{ B^{*}=\frac{ \left( 1+log~q \right) ~Q^{*} \left( p \right) - \left( 1+log~p \right) Q^{*} \left( q \right) }{ \left( 1+log~q \right) ^{2}- \left( 1+log~p \right) ^{2}} }\)

\(\mathrm{ \left( A^{*},B^{*} \right) }\) as linear combinations of the quantiles \(\mathrm{ Q^{*} \left( p \right) }\) and \(\mathrm{ Q^{*} \left( q \right) }\) have an asymptotically binomial distribution which can be obtained by the \(\mathrm{ \delta }\)-method. It does not seem worthwhile to advance this line of thought until it is felt that this “ready-to-wear” approach is a good one.

6 . The estimation of the association between the margins

The fact that we may use an intrinsic estimator of \(\mathrm{ A \left( u \right) }\) or \(\mathrm{ k \left( w \right) }\), even supposing that we are using \(\mathrm{ \bar{A}_{n} \left( w \right) }\) or \(\mathrm{ \bar{k}_{n} \left( w \right) }\), and not \(\mathrm{ A_{n}^{*} \left( u \right) }\) or \(\mathrm{ k_{n}^{*} \left( w \right) }\), does not facilitate the evaluation of the association. It is sufficient to recall that  \(\mathrm{ \rho \left( k \right) =-\frac{6}{ \pi ^{2}} \int _{- \infty}^{+ \infty}log~k \left( w \right) dw~or~ \rho \left( A \right) = \int _{0}^{1}A^{-2} \left( w \right) du-1 }\)to see that they are awkward to obtain except for numerical computation.

As \(\mathrm{ max \left( u,1-u \right) \leq A \left( u \right) \leq 1 }\) and  \(\mathrm{ \frac{max⁡ \left( 1,e^{w} \right) }{1+e^{w}} \leq k \left( w \right) \leq 1 }\), we can use as association indicators the ratios of areas between \(\mathrm{ A \left( u \right) }\) or \(\mathrm{ k(u) }\) and independence suitably normalized:

\(\mathrm{ \xi \left( A \right) =\frac{ \int _{0}^{1} \left( 1-A \left( u \right) \right) d~u}{ \int _{0}^{1} \left( 1-max \left( u,1-u \right) \right) d~u}=4 \int _{0}^{1} \left( 1-A \left( u \right) \right) d~u }\)

and

\(\mathrm{\xi \left( k \right) =\frac{ \int _{- \infty}^{+ \infty} \left( 1-k \left( w \right) \right) d~w}{ \int _{- \infty}^{+ \infty} ( 1-\frac{max \left( 1,e^{w} \right) }{1+e^{w}} ) d~w}=\frac{1}{2~log~2} \int _{- \infty}^{+ \infty} \left( 1-k \left( w \right) \right) d~w }\);

note that the last coefficient could be suggested by the use of the approximation \(\mathrm{ 1-k \left( w \right) to-log~k \left( w \right) }\)which, as  \(\mathrm{ 1/2 \leq k \left( w \right) \leq 1 }\), has the following bounds

\(\mathrm{ 0 \leq -log~k \left( w \right) - \left( 1-k \left( w \right) \right) \leq -log~ \left( 1/2 \right) - \left( 1/2 \right) =.1931471 }\),

and so the renormalization substitutes the coefficient \(\mathrm{ \frac{6}{ \pi ^{2}}=.6079271 }\) by the coefficient \(\mathrm{ \frac{1}{2~log~2}=.7213475 }\); it is immediate, by integrating the relation above, that \(\mathrm{ 0 \leq 2~log~2 \cdot \xi \left( k \right) \leq \frac{6}{ \pi ^{2}} \rho \left( k \right) }\).

We get the estimators

\(\mathrm{ \bar{\xi} _{n} \left( A \right) =4 \int _{- \infty}^{+ \infty} ( 1-\bar{A}_{n} \left( w \right) \left( u \right) ) d~u }\)

\(\mathrm{ =\frac{4}{1+log~n} \sum _{1}^{n} \int _{0}^{1}min⁡ ( \frac{1-u}{1+n~x_{i}^{’}},\frac{u}{1+n~y_{i}^{’}} ) d~u= }\)

\(\mathrm{ =\frac{2}{1+log~n} \sum _{1}^{n}\frac{1}{2+n ( x_{i}^{’}+y_{i}^{’} ) } }\)  

and

\(\mathrm{ \bar{ \xi} _{n} \left( k \right) =\frac{1}{2~log~2} \int _{- \infty}^{+ \infty} ( 1-\bar{k}_{n} \left( w \right) ) d~w }\)

\(\mathrm{ =\frac{1}{2~log~2 \cdot \left( 1+log~n \right) } \sum _{1}^{n} \int _{- \infty}^{+ \infty}\frac{1}{1+e^{w}}min⁡ ( \frac{e^{w}}{1+e^{-x_{i}}},\frac{1}{1+e^{-y_{i}}} ) d~w }\)

\(\mathrm{ =\frac{1}{2~log~2 \cdot \left( 1+log~n \right) } \sum _{1}^{n} [ \frac{1}{1+n~e^{-x_{i}}}log ( 1+\frac{1+n~e^{-x_{i}}}{1+n~e^{-y_{i}}} ) + }\)

\(\mathrm{ \frac{1}{1+n~e^{-y_{i}}}~log~ ( 1+\frac{1+n~e^{-y_{i}}}{1+n~e^{-x_{i}}} ) ] }\)

7 . A worked example:

Consider the data for the flood discharges of the Ocmulgee River in Georgia, USA at Hawkinsville (upstream) and Macon (dowstream), for the 40 years 1910/1949, contained in Table 1 taken from Gumbel and Goldstein (1964).

The margins data can be accepted to have a Gumbel distribution and the ML estimators of the parameters are \(\mathrm{ \hat{\lambda} _{H}=23.709, \hat{\delta} _{H}=15.057, \hat{\lambda} _{M}=26.378~and~ \hat{\delta} _{M}=17.042 }\). The empirical Kolmogoroff-Smirnov statistics for the margins are \(\mathrm{ KS_{H}=.09442~and~KS_{M}=~.06228~with\sqrt[]{n}~KS_{H}=.597~and~\sqrt[]{n}~KS_{M}=.394 }\) accepting thus the Gumbel model for the margins. The better analysis is to use the \(\mathrm{ \hat{V}_{n} }\) statistics for each margin. The medians are \(\mathrm{ \tilde{\mu} _{H}=30.15~and~ \tilde{\mu} _{M}=31.80; }\) note that the years ranked \(\mathrm{ 20 }\) and \(\mathrm{ 21 }\) in both stations are not the same. The number of concordant pairs (for the quadrants method) is \(\mathrm{ N=34 }\); the value given in Gumbel and Goldstein (1964) is \(\mathrm{ N=33 }\) because they estimate the medians through an inefficient estimation of the margin parameters. The correlation coefficient is \(\mathrm{ \rho ^{*}=.942 }\).

The estimation equation \(\mathrm{ 2p^{**}=\frac{N}{n}=\frac{N_{1}+N_{2}}{n} }\) imposes, as \(\mathrm{ \frac{1}{4} \leq p^{**} \leq \frac{1}{2},\frac{1}{2} \leq \frac{N}{n} \leq 1 }\) which is acceptable as \(\mathrm{ {N}/{n}={34}/{40}=.85~and~p^{**}=N/2\,n=.425 }\). Clearly this value of \(\mathrm{ p^{*} }\) rejects independence as can be seen.

As for the mixed model  \(\mathrm{ p_{M} \left( \theta \right) =2^{-2+ \theta /2},max\,p_{M} \left( \theta \right) =p_{M} \left( 1 \right) =2^{-3/2}=.354<.425 }\), the mixed model should be rejected or, accepting it, take \(\mathrm{ \theta =1 }\); the distance between \(\mathrm{ .354 }\) and \(\mathrm{ .425 }\) suggests avoiding the truncation. Although Posner, Rodemich, Ashlock and Lurie (1969) did find an application of the mixed model — see the data and the strip method in their paper — we may expect the model to be useful only for weak correlations.

Consider now the logistic model. As \(\mathrm{ p_{L} \left( \theta \right) =exp \{ -2^{1- \theta }~log~2 \} }\) has the maximum \(\mathrm{ p_{L} \left( 1 \right) =1/2 }\), the model can be accepted. In that case we have \(\mathrm{ \theta ^{**}=.696 }\).  As \(\mathrm{ p_{L} \left( \theta \right) = \theta \left( 2- \theta \right) }\) we get \(\mathrm{ \theta ^{*}=1-\sqrt[]{l- \rho ^{*}}=.759 }\)  close to  \(\mathrm{ \theta ^{**} }\) as \(\mathrm{ \theta ^{*}- \theta ^{**}=.063 }\). But better estimators of \(\mathrm{ \theta }\)  must be found.

From the practical point of view it seems that, for the differentiable models, the logistic one seems to give a good fit in general. But it must be recalled that Tiago de Oliveira (1987) has shown that the logistic and the natural models can be very close and difficult to separate for moderate samples. The latter is easily generalizable for asymmetric situations and even for more dimensions, which is not the case with the logistic model.

Tables

Figures

References