.

Modeling and Parametric Estimation of Bivariate Extremes

Enrique Villa-Diharce

Department of Probability and Statistics Centro de Investigación en Matemáticas A. C. Apdo. Postal 402, Guanajuato, Gto. 36000, México e-mail: villadi@cimat.mx

Abstract

Bivariate extreme value distributions arise as the limiting distributions of renormalizad component-wise maxima. There are considerable restrictions on the dependence structure between the marginals X and Y, which is determined by a dependence function defined as function of the quotient W = Y/(X + Y). This paper considers a characterization of the dependence function in terms of the distribution function of the quotient W. This allows us to study some statistical issues for the bivariate problem using a simpler univariate setting. The maximum likelihood, moments as well as the Bayesian estimators for the parameter of the proposed model are studied, and ilustrated with two examples.

1 Introduction

In the analysis of environmental extreme value data there is often need to study inter-site extreme behavior. Examples are joint flooding at various sea ports and extreme joint concentration of a pollutant at various monitoring sites.

In modelling multivariate extreme data, several asymptotic multivariate distributions with univariate margins in the domains of attraction have been suggested, these include Frechet, Weibull and Gumbel.

The possible forms of the bivariate extreme value distribution are classified according to the chosen margins (see for example, Galambos (1978), Haan and Resnick (1977), Joe (1993), Tawn (1988), Deheuvels (1983), and Resnick (1987)), the parametric family in all cases being infinite dimensional. In this paper we deal with bivariate extreme value distributions having standard exponential margins.

Let (X, Y) be a random vector with exponential standard margins and joint survival function S(x, y) = Pr(X > x, Y > y). It is well-known (see Pickands, 1981) that (X, Y) follows an extreme value bivariate distribution if and only if S(x, y) is of the form

S(x, y) = exp \{−(x + y)α (\frac{y}{x + y} )\} , x > 0, y > 0,  (1)

where

α(w) = ∫_{0}^{1} max{(1 − w)q, w(1 − q)}dH(q),  (2)

with H being a positive finite measure on [0, 1] such that

∫ _{0}^{1} qdH(q) = ∫ _{0}^{1}(1 − q)dH(q) = 1.  (3)

The dependence function α is convex within the region max(w, 1−w) ≤ α(w) ≤ 1, 0 ≤ w ≤ 1.

The dependence structure between X and Y is determined by α(w), where limit cases of independence and total dependence are respectively, α(w) = 1, 0 ≤ w ≤ 1 and α(w) = max(1 − w, w), 0 ≤ w ≤ 1.

Since extremes X and Y are invloved in the dependence function through the ratio W = Y/(X + Y), it is interesting to find out relationship between structure dependence of (X, Y) and the behaviour of W. In the limit cases of independence and total dependence between X and Y, W is Uniform(0, 1), and P (W = 1/2) = 1 respectively.

There is an interesting relationship between the dependence function α(w) and the distribution function F_{W} (w) given by

F_{W} (w) = w + w(1 − w)α′(w)/α(w).  (4)

The relation (4) provides an alternative representation for the dependence function α(w) in terms of F_{W} (w), namely,

α(w) = exp \{∫_{0}^{w} \frac{F (s) − s}{s(1 − s)}ds\} .  (5)

It is possible the use of F_{W} for the modeling and statistical analysis of bivariate extreme value distributions.

The distribution function of the quotient W corresponding to some parametric models of the symmetric dependence function can be obtained from relationship (4). For example we have

Mixed Model (Gumbel, 1960)

α(w) = θw^{2} − θw + 1, 0 ≤ θ ≤ 1,  (6)

F (w) = \frac{−θw^{3} + 2θw^{2} + (1 − θ)w}{θw^{2} − θw + 1}

Logistic Model (Gumbel, 1960)

α(w) = {(1 − w)^{1/θ} + w^{1/θ}}^{θ} , 0 ≤ θ ≤ 1,  (7)

independence and total dependence correspond to θ = 1 and θ → 0, respectively. In this case

F (w) = \frac{w^{1/θ}}{(1 − w)^{1/θ} + w^{1/θ}} = \frac{1}{1 + (\frac{1−w}{w} )^{1/θ}} .

Independence corresponds to

α(w) = 1, 0 ≤ w ≤ 1  (8)

and

F (w) = w, 0 ≤ w ≤ 1.

In the total dependence case

α(w) = max(w, 1 − w), 0 ≤ w ≤ 1  (9)

and

F (w) = \{ \begin{matrix}0 & 0 ≤ w < 0.5 \\ 1 & 0.5 ≤ w ≤ 1 \end{matrix}

Figures 1(a, b) show dependence α(w) and distribution F (w) functions corresponding to independence, total dependence as well as logistic model with mild dependence.

A measure of asymmetry of the dependence between X and Y is given by

E(W) = ∫_{0}^{1} [1 − F (s)] ds = 0.5 − ∫_{0}^{1} s(1 − s)α′(s)/α(s)ds.

Fig. 1.(a) Dependence functions for independence, total dependence and logistic dependence.

Fig. 1.(b) Distribution functions for independence, total dependence and logistic dependence

2 Proposed Parametric model

Consider the setting of defining structure dependence between X and Y through the distribution of the quotient W. In this section we present a new parametric model with one parameter to model dependence.

For the symmetric dependence case the distribution function proposed is

F (w) = \{ \begin{matrix}(1/2)(2w)^{a}, 0 ≤ w ≤ 1/2 \\ 1 − (1/2)(2(1 − w))^{a}, 1/2 < w ≤ 1, 0 < a < ∞ \end{matrix}  (10)

with corresponding density funcion

f(w) = \{ \begin{matrix}a (2w)^{a−1} , 0 ≤ w ≤ 1/2 \\ a (2(1 − w))^{a−1} , 1/2 < w ≤ 1, 0 < a < ∞ \end{matrix}  (11)

This model considers negative dependence between X and Y when 0 < a < 1, as well as positive dependence when 1 < a < ∞. The limit cases a → 0 and a → ∞, correspond to total negative and positive dependence. The case of independence corresponds to a = 1.

If we are sure that dependence between X and Y is positive, then a is in the range [1,∞).

2.1 Measure of Asociation

An usual measure of asociation between random variables X and Y is Kendal’s τ. Association between extreme values (X, Y), measured by τ, is obtained through the dependence function A(w), by Ghoudi, Khoudraji and Rivest (1998),

τ = ∫_{0}^{1} \frac{t(1 − t)}{A(t)}dA′(t).  (12)

It is possible to define τ in terms of distribution function F (w) W, for the quotient W, by Villa (1999)

τ = ∫_{0}^{1} [\frac{F (t)}{t} − 1] [\frac{F (t)}{1 − t} − 1] dt.  (13)

For the proposed model, defined by quotient distribution this association coeficient is given in terms of hipergeometric function F (α, β, γ, w) by (Villa, 1999),

τ = (1 − \frac{1}{a}) + \frac{1}{4a}F (1, 2a, 2a + 1, 1/2) − \frac{1}{a + 1}F (1, a + 1, a + 2, 1/2),

or the series representation of hipergeometric function,

τ = (1 − \frac{1}{a}) + \frac{1}{2} ∑_{i=0}^{∞} \frac{(\frac{1}{2})i}{2a + i} − ∑_{i=0}^{∞} \frac{(\frac{1}{2})i}{a + 1 + i}.  (14)

Fig. 2.1. Relationship between Kendall’s τ and parameter a of proposed

model.

Table 2.1. Relationship between Kendall’s τ and parameter a of proposed model.

\begin{matrix}a & 1.0 & 1.1 & 1.5 & 2 & 3 & 4 & 5 & 10 & 20 & 50 & 100 \\ τ & .0 & .078 & .291 & .445 & .611 & .699 & .754 & .871 & .933 & .973 & .986 \end{matrix}

In Figure (2.1) we show relationship between Kendall’s τ and parameter a of proposed model. We note small increments in a near of zero, correspond to large changes in τ. This relationship is in Table (2.1) for several values of a and τ.

2.2 Dependence function

Using distribution function proposed here, in representation (5) we obtain dependence function,

α(w) = \{ \begin{matrix}(1 − w) exp (∑_{n=0}^{m} \frac{w^{(n+a)}}{(n+a)} + O(|a + 1|^{−m−1})) & 0 ≤ w ≤ 1/2 \\ w exp (∑_{n=0}^{m} \frac{(1−w)^{(n+a)}}{(n+a)} + O(|a + 1|^{−m−1})) & 1/2 ≤ w ≤ 1 \end{matrix} (15)

where the number of terms in the series can be chosen according to the magnitude order O(|a + 1|^{−m−1}), for example, especifying accuracy 10^{−10} , the number of terms depends on the value of parameter a :

\begin{matrix}a & 1.001 & 1.01 & 1.1 & 1.5 & 2 & 3 & 4 & 5 & 6 \\ m & 32 & 31 & 30 & 25 & 20 & 17 & 14 & 12 & 11 \end{matrix}

3 Inference

At consider the dependence structure characterized by distribution function F (w), sample information is on quotients W_{1}, W_{2}, ..., W_{n}. Investigation of dependence of X and Y having a parametric model, is possible using likelihood, moments or Bayesian methods. These procedures are easy because of the simplicity of probability density function of proposed model,

f(w) = \{ \begin{matrix}a2^{a−1}w^{a−1}, 0 ≤ w ≤ 1/2 \\ a2^{a−1}(1 − w)^{a−1}, 1/2 < w ≤ 1. \end{matrix}  (16)
= a [2 min (w, 1 − w)]^{a−1} , 0 ≤ w ≤ 1,

where parameter a is in interval [1,∞).

Our distribution family belongs to exponential family because probability density function (16) can be expressed in exponential form,

f(w; a) = g(a)b(w) exp {ϕ(a)h(w)} ,  (17)

where, g(a) = a, b(w) = 1, ϕ(a) = −(a − 1), h(w) = − log [2 min (w, 1 − w)] .

It is nice that this distribution be in exponential family because share known properties for this extended family.

Therefore exists a complete sufficient statistic for parameter a. If information about a is in random sample {W_{1}, W_{2}, ..., W_{n}},

T (W_{1}, ..., W_{n}) = − ∑_{i=1}^{n} log [2 min (W_{i}, 1 − W_{i})] .  (18)

This statistic T has Gamma distribuion with shape and scale parameters n and a respectively. An inmediate consequence is that expectance and variance for T are E(T) = n/a and V (T) = n/a^{2}, moreover transformation S = (n−1)/T is an unbiased estimator for a. According to Lhemann-Scheffé theorem, the estimator

S(W_{1}, W_{2}, ..., W_{n}) = (n − 1)/T (W_{1}, W_{2}, ..., W_{n})
= −\frac{n − 1}{∑_{i=1}^{n} log [2 min (W_{i}, 1 − W_{i})]},

is an UMBU for parameter a.

Another result of the distribution of T, is that starting from the pivotal quantity 2aT, we obtain interval estimations for a. Because of the distribution for T, pivotal quantity 2aT has a Chi-square distribution with 2n degrees of freedom. According to this a 100(1 − α)% confidence interval for a is

[\frac{κ_{(2n,α/2)}^{2}}{2T} , \frac{κ_{(2n,1−α/2)}^{2}}{2T} ] .

With high probability, left limit of this interval can be less than one, when data come from a population with a small dependence structure, and by random is induced a negative dependence in sample, wich correspond a value for a less than one in proposed model.

3.1 Inference based on Likelihood

In the likelihood setting, sampling information w_{1}, w_{2}, ...w_{n} about proposed model is given by likelihood function

L(a; w_{1}, ..., w_{n}) = a^{n} ∏_{i=1}^{n} [2 min (w_{i}, 1 − w_{i})]^{a−1} , a > 0,  (19)

or equivalently in log-likelihood function

l(a; w_{1}, ..., w_{n}) = n log(a) + (a − 1) ∑_{i=1}^{n} log [2 min (w_{i}, 1 − w_{i})] .  (20)

Maximum estimator for a is obtained maximizing L(a; w_{1}, ..., w_{n}), and in our model, this estimator for a, is

â = −n {∑_{i=1}^{n} ln [2 min (w_{i}, 1 − w_{i})]}^{−1} .  (21)

Observed information is

I(â) = − \frac{d^{2}l(a)}{da^{2}} |_{a=â} = −\frac{n}{â^{2}} ,

while asymptotic variance is

V (â) = a^{2}/n.  (22)

The relative likelihood function R(a) = L(a)/L(â) as well as the log-relative likelihood function r(a) = l(a) − l(â) give a support measurement of data to parameter values a. The distribution of the function −2 log R(a) = −2r(a) converges to a chi squared distribution with one degree of freedom. This property can be used to obtain likelihood-confidence intervals. A likelihood-confidence (1 − α)100%, is defined as

{a > 0 : −2r(a) ≤ χ_{(1−α;1)}^{2}} ,

where χ_{(1−α;1)}^{2} is the chi squared (1 − α/2) quantile with one degree of freedom.

The left edge of the likelihood-confidence interval can be located to left one when the dependence patern of data is negative or weak positive.

3.2 Estimation of a by Method of Moments.

Note that

V (W) = E(W^{2}) − (EW)^{2} = \frac{1}{2(a + 1)(a + 2)}.

From equality V (W) = s^{2} we obtain de moments estimator

â = \frac{−3 + (1 + 2/s^{2})^{1/2}}{2} ,  (23)

where s^{2} is the second sample moment centered at the population mean 1/2,

s^{2} = \frac{1}{n} ∑_{i=1}^{n}(w_{i}^{2} − 1/4).

3.3 Bayesian Approach

The parameter space given at (10), is defined in extended form, including values a < 1. If we are sure that dependence structure is not negative, we must estimate a constrained to (a ≥ 1). Interval estimates also exclude less than one values. It is posible include this constrain using the Bayesian approach. We explore the Bayesian approach about the parameter of the proposed model. The information useful to make inferences about the parameter is in the a posteriori density function,

f(a|w_{1}, w_{2}, ..., w_{n}) = \frac{f(w_{1}, w_{2}, ..., w_{n}|a)f_{0}(a)}{∫_{A} f(w_{1}, w_{2}, ..., w_{n}|a)f_{0}(a)da},

where f(w_{1}, w_{2}, ..., w_{n}|a) is the joint density of data w_{1}, w_{2}, ..., w_{n} given the value a of the parameter, f_{0}(a) is the a priori density function.

The information about a is related to the dependence pattern between extremes X and Y, wich are positively associated. After the translation of this information to an a priori density function, we discuss the conjugated a priori density function, because the mathematical convenience to obtain the a posteriori density. Such family is the conjugated a priori.

As we pointed in (17), the density function proposed f(w; a) = a[min(w, 1 − w)]^{a−1}, 0 ≤ w ≤ 1, where the parameter a is in (0, ∞). it is a regular exponential family because the support set (0,1) is independent of the parameter a. According to Bernardo and Smith (1994), when the density function has the prvious form, the a priori density function is,

f_{0}(a) = [K(τ_{0}, τ_{1})]^{−1} [g(a)]^{τ^{0}} exp {ϕ(a)τ_{1}} ,

where g(a) = a, ϕ(a) = −(a − 1) and vector (τ_{0}, τ_{1}) is such that

K(τ_{0}, τ_{1}) = ∫_{1}^{∞}[g(a)]^{τ^{0}} exp {ϕ(a)τ_{1}} da < ∞.

In our case, the a priori distribution is Gamma(τ_{0} + 1, τ_{1}) truncated to left at 1, with density function ,

f_{0}(a|τ_{0}, τ_{1}) = [\frac{(τ_{1})^{τ_{0}+1}}{Γ(τ_{0} + 1)S(1, τ_{0} + 1, τ_{1})}] a^{τ^{0}} exp \{−τ_{1}a\} ,  (24)

where S(1, τ_{0} + 1, τ_{1}) is the survival Gamma(τ_{0} + 1, τ_{1}) evaluated in one. The a posteriori density is,

f(a|w, τ) = \frac{a^{n} [ ∏_{i=1}^{n} 2 min(w_{i}, 1 − w_{i})]^{a−1} [\frac{(τ_{1})^{τ_{0}+1}a^{^{τ}^{0}} exp\{−τ_{1}a\}}{Γ(τ_{0}+1)S(1,τ_{0}+1,τ_{1})}]}{∫_{1}^{∞}a^{n} [ ∏_{i=1}^{n} 2 min(w_{i}, 1 − w_{i})]^{a−1} [\frac{(τ_{1})^{τ_{0}+1}a^{^{τ}0} exp\{−τ_{1}a\}}{Γ(τ_{0}+1)S(1,τ_{0}+1,τ_{1})}] da} (25)

= \frac{a^{τ_{0}+n} exp \{−(t(w) + τ_{1})a\}}{∫_{1}^{∞}a^{τ_{0}+n} exp \{−(t(w) + τ_{1})a\} da}

= \frac{[(t(w) + τ_{1})]^{τ_{0}+n+1} a^{τ_{0}+n} exp \{−(t(w) + τ_{1})a\}}{Γ(τ_{0} + n + 1)S(1, τ_{0} + n + 1, τ_{1} + t(w))},

where t(w) is the sufficient statistic (18). This is the density function of the Gamma(τ_{0} +n+1, τ_{1} +t(w)) distribution truncated to left at one. Expectation and variance for a are

E(a|w, τ) = (\frac{τ_{0} + n + 1}{(t(w) + τ_{1})}) C_{1}, y  (26)

V (a|w, τ) = (\frac{τ_{0} + n + 1}{(t(w) + τ_{1})^{2}} ) C_{2},

where,

C_{1} = (\frac{S(1; τ_{0} + n + 1, t(w) + τ_{1})}{S(1; τ_{0} + n, t(w) + τ_{1})}) ,

C_{2} = \{n (\frac{S(1; τ_{0} + n + 2, t(w) + τ_{1})}{S(1; τ_{0} + n, t(w) + τ_{1})}) − (n − 1) (\frac{S(1; τ_{0} + n + 1, t(w) + τ_{1})}{S(1; τ_{0} + n, t(w) + τ_{1})})^{2}\} .

If we consider E(a|w, τ) as a Bayesian estimator for the parameter a, such estimator is a corrected version of the maximum likelihood estimator n/t(w), wich is also obtained considering the non informative of the Gamma density, where the parameters are, τ_{0} = −1, τ_{1} = 0.

Acording to the previous information we have about the dependence parameter a, extremes X y Y are positively related, with Kendall’s association coefficient 0 ≤ τ ≤ 1.

Using relationship (14) between parameter a of the model and the association coefficient τ, we obtain the approximated relationship, τ = 1 − 1/a. Taking a as well as τ as random variables, we model the information about positive dependence between extremes through a Uniform (0, 1) for τ firstly, then the priori density f_{0}(a) for a, is obtained making a variable change,

f(a) = a^{−2}1_{[1,∞)}(a).  (27)

Considering this a priori density, the correspondent a posteriori density is,

f(a|(w_{1}, ..., w_{n})) ∝ a^{n−2} exp{−at(w)}1_{[1,∞)}(a).

The exact form of this truncated distribution obtained normalizing the density with the survival function S(1; n−1, t(w)) evaluated at one for the Gamma(n− 1, t(w)) distribution. Thus the truncated a posteriori density resulting is,

f(a|(w_{1}, ..., w_{n})) = [1/S(1; n − 1, t(w))]a^{n−2} exp{−at(w)}1_{[1,∞)}(a). (28)

It is interesting to note that this is a member of the conjugated family prviously founded, with parameters τ_{0} = −2 and τ_{1} = 0.

Now the expectation and variance for a are,

E(a|w) = \frac{n − 1}{t(w)} (\frac{S(1; n, t(w))}{S(1; n − 1, t(w))})  (29)

V (a|w) = \frac{n − 1}{t^{2}(w)} \{n (\frac{S(1; n + 1, t(w))}{S(1; n − 1, t(w))}) − (n − 1) (\frac{S(1; n, t(w))}{S(1; n − 1, t(w))})^{2}\} .

A Bayesian point estimator for parameter a is the expectation E(a|w), while a minimum longitude Bayesian 100(1 − α)% confidence estimator of a is the maximum posterior density interval (a_{1}, a_{2}), discussed in Box and Tiao (1973), wich is obtained using the posterior density f(a|w), defining the limits a_{1} as well as a_{2} through constraints,

∫_{a_{1}}^{a_{2}}f(a|w)da = 1 − α,

f(a_{1}|w) = f(a_{2}|w).

Using the the a posteriori distribution for a, we are able to do forecasts about different interes quantites or about future observations, starting with the actually observations.

The excess posteriori joint probability is

P (X > x, Y > y|w_{1}, ..., w_{n}) = ∫_{0}^{∞}P (X > x, Y > y|a)f(a|w_{1}, ..., w_{n})da. (30)

This estimation considers variability of a as well as the residual variability of X, Y, when a is known.

The numerical calculation previous forecast is obtained generating the realizations sequence a_{1}, a_{2}, ..., a_{m} using the a posteriori distribution for a, then we calculate the mean of excess probability evaluated for the different a realizations,

\frac{1}{m} ∑_{i=1}^{m} P (X > x, Y > y|a_{i}).  (31)

3.4 Comparison of estimators.

A simulation study was conducted to compare, in terms of bias and meansquared error, the maximum likelihood, moments and Bayesian estimators.

While the maximum likelihood â_{L} and Bayesian â_{B} estimators are defined in terms of the complete sufficient statistic T (W_{1}, ..., W_{n}), the moments estimator â_{M} is a function of the sampling variance s^{2}, because of that, the last estimator has the greatest variance. It’s an easy task to obtain the expression for the variance of â_{L}, but, for â_{M} there is not a closed expression for the variance. Jackknife method can be used to estimate the variance of moments estimator â_{M}. In the simulation study we consider four levels of dependence, very low, low, mild and high, corresponding to the values 1.1, 1.5, 2.0 and 10 for the a parameter of the proposed model. In terms of the association coeficient Kendal’s τ, the a values correspond to the association levels, .078, .291, .445 and .871, respectively. The bias and mean squared error are shown in Tables 3.1 and 3.2 respectively, for the estimators under comparison. Table 3.1 shows a positive bias for all the compared estimators, for the considered four dependence patterns. While for the maximum likelihood and moments estimators diference is small, Bayesian estimator has a less bias than other two. This bias reduction is explained by the relationship between Bayesian and maximum likelihood estimators, because the former is less than the later by a correction factor less than one.

â_{B} = (\frac{n − 1}{1} \frac{S(1, n, t(w))}{S(1, n − 1, t(w))}) â_{L}.

Table 3.1. Bias for maximum likelihood, moments, and Bayesian estimators, for the a parameter of the proposed model.

\begin{matrix}Bias(â_{L}) \\ Bias(â_{M}) \\ Bias(â_{B}) \end{matrix}

Very low dep. a = 1.1

\begin{matrix}Low \\ dep. \\ a = 1.5 \end{matrix}

Medium dep. a = 2.0

\begin{matrix}High \\ dep. \\ a = 10.0 \end{matrix}

n = 20  \begin{matrix}.058561 \\ .048968 \\ .173216 \end{matrix}  \begin{matrix}.076316 \\ .072701 \\ .060026 \end{matrix}  \begin{matrix}.110961 \\ .111438 \\ .020554 \end{matrix}  \begin{matrix}.500462 \\ .758639 \\ −.024561 \end{matrix}

n = 50  \begin{matrix}.023631 \\ .018934 \\ .087814 \end{matrix}  \begin{matrix}.030218 \\ .025470 \\ .009088 \end{matrix}  \begin{matrix}.048219 \\ .049051 \\ −.0022741 \end{matrix}  \begin{matrix}.193516 \\ .283505 \\ −.010355 \end{matrix}

Table 3.2. Mean squared error for the estimators for the various cases of dependence.

\begin{matrix}ECM(â_{L}) \\ ECM(â_{M}) \\ ECM(â_{B}) \end{matrix}

Very low dep. a = 1.1

\begin{matrix}Low \\ dep. \\ a = 1.5 \end{matrix}

Medium dep. a = 2.0

\begin{matrix}High \\ dep. \\ a = 10.0 \end{matrix}

n = 20  \begin{matrix}.079089 \\ .096321 \\ .059434 \end{matrix}  \begin{matrix}.141872 \\ .166391 \\ .095653 \end{matrix}  \begin{matrix}.260854 \\ .290876 \\ .210330 \end{matrix}  \begin{matrix}6.368880 \\ 7.873536 \\ 5.522475 \end{matrix}

n = 50  \begin{matrix}.027084 \\ .031595 \\ .018528 \end{matrix}  \begin{matrix}.050199 \\ .054570 \\ .042802 \end{matrix}  \begin{matrix}.135848 \\ .144031 \\ .128236 \end{matrix}  \begin{matrix}2.216371 \\ 2.569900 \\ 2.092745 \end{matrix}

The estimators mean squared error, shown in Table 3.2, exhibit less variability for the Bayesian estimator. The â_{B} stability is greater when the dependence decreases, that is when a is closer to one. The explanation for such behavior is that the information about positive dependence, incorporated in the Bayesian approach, cancels the posibility of values less than one for that estimator, wich corresponds to negative dependence. This left boundary is not included in â_{L} and â_{M} estimators, wich show greater variability in the case of low dependence pattern.

4 Examples

In this section we present two examples corresponding to high dependence and independence of two componentwise maxima. We plot in each case the empirical distribution function of W.

4.1 Example 1. Flood Data

Maximum annual discharge of the Ocmulgee river. The data consist of maximum annual discharge for the years 1910-1948 for the Ocmulgee river in Georgia measured at two different stations: Hawkinsville and Macon, according to the Geological Survey (Carter, 1951).

Fig. 4.1. Maximum annual discharge at two stations on the Ocmulgee river. Transformed data to exponential standard.

Fig. 4.2. Empirical distribution function for the quotients of the maximum annual discharges at two stations of Ocmulgee river.

Gumbel distributions were fitted with location and scale parameters given, respectively, by (18.01, 26.76) for Macon and (15.05, 23.85) for Hawkinsville.

Subsecuently, the data were transformed to standard exponentials. These transformed data are presented in the scatter plot Fig. 4.1.

The high level of dependence observed in the scatter plot is also noted in the empirical distribution function of the quotients W of maximum annual discharges at the two stations, shown in Fig. 4.2, wich display high frecuency of values near to .5, corresponding to high association between both extreme discharges. An additional indication of the high level of dependence is the high value of the parameter a of the proposed model (10), shown in Table 4.1.

Table 4.1. Estimates for the parameter a of the proposed model (10).

â  V̂ ar(â)
Maximum likelihood  4.7934  .5744
Moments  4.7850

The Bayesian analysis also concludes that the value of the parameter a is high, as is shown in the Table 4.2 wich display confidence interval of minimum length for different probability levels. The point estimation for a is according to (26) is a = 4.5736.

Table 4.2. Bayesian confidence intervals for three probability levels.

Confidence 90 95 99
Confidence interval (3.442, 5.882) (3.253, 6.163) (2.903, 6.735)

We calculate the conditional probability of exceeding the thresholds in both locations, given in (30). In Table 4.3 are shown the probabilities of maximum annual discharges do not exceedes thresholds u y v, P (U < u, V < v|w) = P (X > x, Y > y|w) = S(x, y|w), for several pairs of threshold values u, v, in the two locations. In whath follows U and V represent, respectively, the maximum annual discharges at hawkinsville and Macon. In Table 4.3 also are shown the standard exponential values x and y corresponding to Gumbel values u and v.

Table 4.3. Probabilities of not exceeding different thresholds u y v in both locations.

\begin{matrix}u & v & x & y & S(x, y|w) \\ 50 & 50 & .2750 & .1760 & .7448 \\ 80 & 80 & .0519 & .0239 & .9479 \\ 100 & 100 & .0171 & .0063 & .9827 \\ 150 & 150 & .0010 & 2.4 × 10^{−4} & .9989 \\ 200 & 200 & 6.6 × 10^{−5} & 8.2 × 10^{−6} & .9999 \end{matrix}

In the context of the present example, and following to Tawn (1988), it is of interest to estimate also the following probabilities for given threshold values, u and v :

p_{1} = P (U > u), p_{2} = P (V > v)
p_{1o2} = 1 − P (U < u, V < v)
= 1 − [(1 − p_{1})(1 − p_{2})]^{A( \frac{log(1−p_{2})}{log(1−p_{1})+log(1−p_{2})} ),}
p_{1&2} = P (U > u, V > v) = p_{1} + p_{2} − p_{1o2},
p_{12} = P (V > v|U > u) = \frac{p_{1&2}}{p_{1}} ,
p_{−12} = \frac{p_{2} − p_{1&2}}{1 − p_{1}} .

Table 4.4 shows the above probabilities for the threshold values of u = v = 120 for the estimators of α(w), using (15) with substitution of parameter a for the maximum likelihood, moments and Bayesian estimates, as wel as the nonparametric estimator proposed by Rojo, Villa-Diharce, and Flores (2001). This table includes also the independence case.

Table 4.4. Marginal, joint and conditional probabilities of exceeding thresholdad levels u = v = 120 at Hawkinsville and Macon.

Max.Lik. Mom.  Bayes  Nonparam. Independence
A(w) 0.6622  0.6623  0.6635  0.6749  1.0000
p_{1}  0.02188 0.02188 0.02188 0.02188  0.02188
p_{2}  0.01221 0.01221 0.01221 0.01221  0.01221
p_{∗}  0.022533 0.022537 0.02257 0.022947  0.033811
p^{∗}  0.011565 0.011562 0.011522 0.011132  0.00256
p_{12}  0.528381 0.528227 0.526383 0.509040  0.012210
p_{−12} 0.000658 0.000662 0.000703 0.001102  0.012210

In case of wrongly assuming independence between maximum annual discharges, when in fact there is a strong dependence we obtain substantial errors in the estimation of probabilities related to excess in both locations. It is interesting to note that probabilities in Table 4.4 depend on the value of dependence function calculated in the value

w = \frac{log(1 − p_{2})}{log(1 − p_{1}) + log(1 − p_{2})} = 0.3570,

observing the importance of the value of the dependence function α(w) in any value, not just in the value w = 1/2, wich corresponds to te case of equal values of marginals, p_{1} = p_{2}.

4.2 Example 2. Environmental Data

The data set is constructed from extreme values in the ozone concentration and wind velocity within episodes that were determined in studing (Castro, 1994) the joint distribution of extreme values of ozone and several meteorological variables registered in the monitoring Pedregal station, at the Environmental Monitoring Automatic Network of Mexico City.

While Weibull distribution for minima was fitted to wind velocity, Gumbel distribution was fitted to maximum ozone concentration. In consequence velocity data were transformed to standard exponentials using an increasing transformation, while the ozone data were transformed to standard exponentials using a decreasing transformation. Because of the transformations used to obtain standard exponential data, a positive dependence between transformed data, corresponds to negative dependence between original data.The transformed data are presented in the scatter plot Fig. 4.3, in this graph we note a weak positive dependence pattern .

The low level of dependence observed in the scatter plot is also present in the empirical distribution function of the quotients W of wind and ozone transformed to standard exponentials, shown in Fig. 4.4, wich display an almost uniform distribution for quotients, but at the begining we note that empirical distribution is below the diagonal, while at the end the empirical distribution above the diagonal, this behavior correspond to weak positive dependence.

The low level of positive dependence between transformed data is noted in the value of estimates of the parameter α greater than one but very closed to the value a = 1,wich corresponds to independence.

Fig. 4.3. Maximum ozone concentration and minimum wind velocity

transformed data to exponential standard.

Fig. 4.4. Empirical distribution function for the quotients of the minimum wind velocity and maximum ozone concentration, transformed to exponentials.

In the setting of quotients used in this paper, the independence hipotesis corresponds to the hipostesis a = 1. Hence we investigate the dependence between extremes, through the behavior of estimate â.

Table 4.5. Parameter estimates for the proposed model.

â  V̂ ar(â)
Maximum likelihood  1.3996  .0408
Moments  1.3208

Table 4.5 shows parameter estimates, obtained by maximum likelihood and moments. The values of â greather than one but small, correspond to weak positive dependence between extremes. According to the confidence intervals for a, shown in Table 4.6 we can not reject the hipotesis of positive dependence between extremes transformed to standard exponentials, because the value a = 1 is located to left of interval.

Table 4.6. Confidence intervals for a, based on normal approximation likelihood.

Confidence 80 90 95
Confidence interval (1.1407, 1.6585) (1.0673, 1.7319 (1.0036, 1.7955)

Hence, the original extremes of wind velocity and ozone concentration, are negatively related. This is a logical relationship, because the low wind velocity is associated to unfavourable conditios for ozone generation, and high values of wind velocity are associated to ozone spread.

5 Conclusions

The dependence structure in bivariate extreme value distributions of (X, Y) can be studied through the distribution for the W = Y/(X + Y), so a bivariate dependence structure problem can be viewed in a univariate setting. This is particular useful for testing the hypotheses of independence or total dependence between X and Y.

Acknowledgements We acknowledge to CONACYT for their support given through the research

project I32824-E.

References

[1] Bernardo, J., M., and Smith, A. F. M. (1994)Bayesian Theory, J. Wiley.

[2] Box, G. E. P., and Tiao, G. C. (1973)Bayesian Inference in Statistical Analysis. Addison-Wesley.

[3] Carter, R. W. ”Floods in Georgia”, Geological Survey Circular 100, March, 1951.

[4] Castro, J. I. (1994). Distribución conjunta de valores extremos de ozono y de algunas variables meteorológicas. Revista de Estadística, Vol. VI, No. 8, 83-98.

[5] Deheuvels, P. (1983). Point processes and multivariate extreme values. J. Multivariate Anal. 13, 257-271.

[6] Galambos, J. (1978). The Asymptotic Theory of Extreme Order Statistics. New York: Wiley.

[7] Gumbel, E. J. (1960). Distributions de valeurs extremes en plusieurs dimensions. Publ. Inst. Statist. Paris 9,171-3

[8] Gumbel, E. J. (1964). Analysis of empirical bivariate extremal distributions, J. Amer. Statist. Assoc. 59, 794-816.

[9] Ghoudi, K., Khoudraji, A., and Rivest, P. (1998). Propriétés statistiques des copules de valeurs extrêmes bidimennsionelles, The Canadian Journal of Statistics, 26, 187-197.

[10] Haan, L. De and S. I. Resnick (1977). Limit theory for multivariate sample extremes. Zeit. Wahrscheinl. Theorie 40, 317-337.

[11] Joe, H. (1993). Parametric family of multivariate distributions with given margins. J. Multivariate Anal. 46, 262-282.

[12] Joe, H. (1994). Multivariate extreme-value distributions with applications to environmental data. Canad. J. Statist. 22, 47-64.

[13] Pickands, J. (1981). Multivariate extreme value distribution. In Proc. 43rd Session I.S.I. (Buenos Aires), pp. 859-878.

[14] Resnick, S. I. (1987). Extreme Values, Regular Variation and Point Processes. New York: Springer Verlag.

[15] Rojo,R. J. Villa-Diharce, E., and Flores, M. (2001). Nonparametric Estimation of the Dependence Functionin Bivariate Extreme Value Distributions. Journal of Multivariate Analysis, 76, 159-191.

[16] Tawn, J. A. (1988). Bivariate extreme value theory-models and estimation. Biometrika 75, 397-415.

[17] Villa E. (1999). Estimación de la Función de Dependencia de la Distruibución Bivariada de Valores Extremos. Ph. D. Thesis, CIMAT.