ESTIMATION W I T H MULTIVARIATE E X T R E M E V A L U E DISTRIBUTIONS, WITH APPLICATIONS TO ENVIRONMENTAL DATA By YOUSHI C H E N B.S c , Beijing Institute of Technology, 1983 M.S c , Beijing Institute of Technology, 1986 A THESIS S U B M I T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES D E P A R T M E N T O F STATISTICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA September 1991 ©Youshi Chen, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may department or by his or her representatives. be granted by the head of my It is understood that copying or publication of this thesis for financial gain shall not be allowed without my permission. Department of The University of British Columbia Vancouver, Canada Date DE-6 (2/88) written Abstract Several parametric families of multivariate extreme value distributions (Hiisler and Reiss 1989, Tawn 1990, Joe 1990a, 1990b) have been proposed recently. Applications to multivariate extreme value data sets are needed to assess the adequacy of the known families in their fit to data. Different families are compared in their range of multivariate dependence and their ease of use for maximum likelihood estimation. Some useful conclusions have been made from experience with several environmental data sets. n TABLES OF CONTENTS Abstract ii Table of Contents iii List of Tables iv Acknowledgement v Chapter 1. Introduction 1 Chapter 2. Multivariate Extreme Value Distributions 3 2.1 Introduction 3 2.2 Multidimensional Distributions 3 2.3 Copulas 6 2.4 Univariate Extremes and the Generalized Extreme Value Distribution 7 2.5 Multvariate Extreme Value Distributions 11 2.6 Some Relation Between Extreme Value Theory and Central Limit Theory 18 2.7 Parametric Families of Multivariate Extreme Value Distributions 20 2.8 Estimation for Multivariate Extreme Value Distribution 23 Chapter 3. Maximum likelihood Estimation With Husler-Reiss Model 25 3.1 Bivariate Case 26 3.2 Trivariate Case 28 3.3 Four Dimensional Case 30 3.4 Computer Implementation of the Husler-Reiss model 32 3.5 Data Analysis 32 Chapter 4. Analysis with Missing Data 44 Chapter 5. Conclusion and Further Research 51 Bibliography 52 in List of tables Table 1. Bay Area-Ozone. (n=120) 38 Table 2. GVRD-Sulphur Dioxide. (n=70) 39 Table 3. GVRD-Nitrogen Dioxide. (n=84) 40 Table 4. GVRD-Ozone. (n=84) 41 Table 5. GVRD-Station 5 (Confederation Park). (n=64) 42 Table 6. GVRD-Station 7 (Anmore). (n=76) 42 Table 7. GVRD-Station 9 (Rocky Point Park). (n=63) 43 Table 8. GVRD-Sulphur Dioxide. (n=104) 48 Table 9. GVRD-Nitrogen Dioxide. (n=84) 49 Table 10. GVRD-Ozone. (n=84) ...50 iv ACKNOWLEDGMENT I would like to express my deep gratitude to my supervisor, Professor Harry Joe, for his constant encouragement, guidance and assistance in my study and in the producing of this thesis. Very special thanks are due Professor Mohan Delampady for reading the manuscripts and making constructive comments. Finally, I would like to thank the Department of Statistics, The University of British Columbia, for supporting me with financial aid and encouragement. v Chapter 1 Introduction Extreme value theory has a great number of applications. We list two cases where the largest or the smallest "measurement" are of interest. 1. A i r Pollution. A i r pollution concentration is expressed in terms of proportion of a specific pollutant in the air. Concentrations are recorded at equal time intervals, and it is required by law to keep the largest annual concentrations below given limits. 2. Natural Disasters. Floods, heavy rains, extreme temperatures, extreme atmospheric pressures, winds and other phenomena can cause extensive human and material loss. Communities can take preventive action to minimize their effects even if such disasters cannot be completely avoided. In dams, dikes, canals, and other structures the choice of building materials and methods of architecture can take some of these potential disasters into account. Engineering decisions that confront such problems should be based on a very accurate theory, because inaccuracies can be very expensive. For example, a dam built at a huge expense may not last long before collapsing. In example 1, at each time point, there may be a vector X of measurements, consisting of concentrations of a pollutant at several air quality monitoring stations in a network or consisting of concentrations of several pollutants at a given monitoring stations. In example 2, there may be a vector X of measurements at several positions in a spatial grid. In the examples, we have over time a number n of random measurements X i , ' . . . , X „ , and the behaviour of either Z = m a x ( X i , . . . , X ) or W n n n = m i n ( X i , . . . , X ) can be of n interest. Here max and min refer to componentwise maxima and minima; see Chapter 2 for details. Statistical inferences of the upper or lower tail of the multivariate distribution of X can be made from using extreme value theory to model the multivariate distribution of 1 Z or W . n n There has been much recent research in multivariate extreme value theory; see Galambos (1987, Chapter 5) and Resnick (1987, Chapter 5) for some general theory. The set of possible of multivariate extreme value distributions has been shown to be of infinite dimension. Smith, Tawn and Yuen (1990) study a density estimation based method for nonparametrically fitting a multivariate extreme value distribution. Except in the bivariate case, the method was not too promising. A n alternative is to find a finite-dimensional parametric subfamily of the set of multivariate extreme value distributions. A good subfamily would be "dense" in the complete family. Parametric families have recently been proposed by Hiisler and Reiss (1989), Tawn (1990), Joe (1990a,b), Coles and Tawn (1991). Experience is now needed to assess the adequacy of fit of these families to multivariate extreme value data sets. Joe (1990b) did some work on this but did not include the Hiisler and Reiss family in his comparisons. One goal of this thesis is to do maximum likelihood estimation with the Hiisler and Reiss family. This turns out be very difficult and details are given in Chapter 3, where comparisons are also made with results in Joe (1990b). Another goal of the thesis is to deal with missing values in computing maxima. Joe (1990b) deleted cases where there were many missing values so that maxima could not be properly computed. These maxima were assumed to be missing at random. Here, in Chapter 4, we treat these cases as right-censored maxima, that is the actual maxima is known to be above the maxima of all non-missing observations. The thesis proceeds next (in Chapter 2) with definitions and properties of general multivariate distributions and of multivariate extreme value distribution. The last chapter (Chapter 5) of this thesis makes some final conclusions and points out a direction for further research. 2 Chapter 2. Multivariate Extreme Value Distribution 2.1 Introduction In this chapter we describe some general properties of multidimensional distributions and copulas. Then we introduce some univariate extreme value theory and the generalized extreme value distribution. Following this, we concentrate on multivariate extreme value distributions, the main topic of this chapter. 2.2 Multidimensional Distributions Let us define the p-dimensional random variable X as the vector X = (Xi,... The distribution function -F(x) = F(xi,..., P x ) is defined as p F(x) = P ( X < x) = P(X where X < Y means Xi < Yi, ,X ) <x ...,X < X u x ), p p for 1 < i < p. For future reference, we also define X.+ Y = (Xi + Y ,...,X p X Y = (X\Yi,..., p 1 + Y ), p X Yp), and X/Y = (X /Y ,...,X /Y ). 1 1 P P Random variables with the same distribution as X will be denoted by X i , . . . , X , p and the components of X j by Xij, that is, Xij is the ith component of X j , or X j = {Xij, • • • ,X j). p The order statistics of the ith. component (Xn,... ,Xi ) n x<£ <xW <...<xj \ n and we denote Wi — XJ;l\ Z{ = 3 X$. are In extreme value theory, an objective is to investigate the existence of the asymptotic distribution of W = n For (W ,... W ) l 1 Z p — (Zi, • • • > Zp). n this, we need some theory of multidimensional distribution functions. A multidimensional distribution function F ( x ) of a random vector X has the following elementary properties: PI (Bounds): 0 < F ( x i , . . . , x ) < 1 Vxi,..., x. P2 (Monotonicity): F is nondecreasing in each of its arguments x,-, 1 < i < p. p P3 (Limits): linx^-oo F ( x p ...,x ) = 0. p u P4 (Marginal distributions): If Xj —• +oo, then F ( x ) tends to an (p—l)-dimensional distribution, which is the distribution of the vector obtained by removing its jth component. We can obtain the univariate marginal distribution FJ(XJ) by letting each x^, i ^ j tend to +oo. P5 (Density): If F has derivatives of order p, then f = >0, dx\ • • •v dx. and / is the probability density function. P6 (Rectangles): P(a <X <b ,...,a <X <b ) 1 l 1 p where for a subset I of { 1 , . . . p p = F(b ,b ,...,b )-J2Fi 1 2 ^F i j p p + ... + i=l + (-lYF(a ,...,a )>0, 1 p F j is the value of F ( s i , . . . , Sk) with s; = a; if i 6 / , and Si = b{ if i £ I. P7 (the Frechet Bounds): p max(0, J2 i( i) F x ~ P + !) < ( ^ F x ••• x ) < min(Fi(a;i),..., p F (x )). p p If Fi,...,F are all continuous, then the upper bound corresponds to distribution of p (Xi, F ~ F i ( X i ) , . . . , F ^ F ^ X i ) ) . In general the lower bound is not a multivariate dis1 2 tribution function for p > 3. For example, suppose F is a three-dimensional distribution function with marginal functions satisfying Fi(xi) F (x ) 3 = 1/2, ^2(^2) = 1/2 and = 1/2. Then 3 F(+oo, +00, - Fi(ari) - F2(x2) +00) + max(0, F ^ ) max(0, F2{x2) = + F {x ) 2 2 + F (x ) 3 1-3/2 =-1/2 3 - F (x ) 2 3 - 1) + max(0, F i ( x i ) + F3{x3) - 1) - max(0, F f a ) + F2(x2) - 1) + + F (x ) 3 3 - 2) <0, which contradicts Property 6. P8 (Survival function): Let Bj be the event (Xj > Xj) and j(k) = (ji,... ,jk) be a subset of { 1 , . . . ,p}, then let G j ( k ) (x ,..., h X j k ) =P(B n • • • n B ), h 5 (x) = 1, 0 jk and k(x) S For = G i<ji<h-ik<p h ( h • •• > p > 2, then F(x ,...,x ) 1 = p x p > x < <Pk 3k)> 1 J2(- ) k(xi,...,x ). l ks p k=0 In addition, for any integer 0 < q < (p — l ) / 2 , E ( - i ) ^ W < f(x) < E ( - i ) ^ ) . f c f c fc=0 k-0 The proof of most of the above results is straightforward; for the proof of properties P7 and P8, readers are referred to Galambos (1987). 5 2.3 Copulas Assume that F(x) is a p-dimensional distribution function with univariate marginals Fi(x{), for 1 < i< p. Let C(y) be a p-dimensional function over the unit cube 0 < yi < 1, 1 < ^ < P, and such that it increases in each of its variables and F(x ,...,x ) 1 = p C[F (x ),...,F (x )]. 1 1 p p Then the function C(y) is called a copula of JP(X). If C is a copula, then G(Vu•••,y ) P = C[G (y ),...,G (y )] 1 1 p p is a multivariate distribution with univariate margins G,-, for i — 1, ...,p, where Gi is an arbitrary univariate distribution function. When needed, we write CF = CF^Y) = C(y) for ^(x). Note that some people call C(y) a dependence function. If F(x\,..., x ) is a continuous p-variate cumulative distribution function with unip variate margins Fj(xj),j = 1,... ,p, then the associated copula is C( ,..., UL tip) = F(Fr (u ),..., F -\u )), 1 1 p 0 < m < 1, i = 1,... ,p. p This is a multivariate distribution with uniform (0,1) margins. Returning to extreme values, the ith. marginal of F (x) is F^Xi), n F ( ,..., xp) = CFn[F?( ,..., n Xl F?(x )]. Xl p On the other hand, we have F (x ,..., p) = n 1 X C [F(x ),...,F(xp)}. F 1 A comparison of these last two equations leads to CF4y) = c [yl F 6 / n ,...,yl / n }. and thus 2.4 Univariate Extremes and the Generalized Extreme Value Distribution Let Xi,...,X denote independent and identically distributed (iid) random vari- n ables, with cumulative distribution function F. Let = max(Xi,..., X ), Z n By W n n = min(X ,..., X ). a n the iid assumption, we have H (x) = P(Z < x ) = n F {x) n n and L (x) = P(W n We seek conditions on F(x) {d } n < x) = 1 - (1 - n F(x)) . n to guarantee the existence of sequences { a } , {&„}, {c }, n n of constants such that, as n —> oo lim H J a n n—*oo + b x) = H(x) n and lim L (c n exist for all continuity points of H ( x ) + d x) n = L(x) n and L(x) respectively, where H ( x ) nondegenerate distribution functions. Equivalently, {a }, n that (Zn — a )/b n n and (Wn — cn)/dn {&„}, and L(x) {c } and { d } n n are are such converge in distribution as n —> oo. A sufficient condition (see Galambos, 1987) is given next. Assume that there are sequences a , c , b > 0, and d > 0, of real numbers such n n n n that, for all x and y, n[l — F(an lim + b x)] = u(x), n n—>-+oo lim nF(c + d y) n n = w(y) n—*-+oo exist. Then lim P(Z n < a + b x) n n = exp[—u(x)] and Jim P(W - 1 - < c + d y) n n n exp[-w(y)]. Example 1. ( The Exponential Distribution). Let Xi,... ,X n be iid random variables with common distribution function F(x) = 1 - e~ , x > 0. x Then 1 - F{a + b x) = e - e a n n h n X n . In order to satisfy n[l — F(a lim + b x)\ — u(x), n n n—•+oo we can choose a = logn. and b = 1. Hence, u(z) n n lim On P(Z and thus, = e~ , z < logre+ z) = exp(—e~ ). z n the other hand, we can choose c = 0 and d = 1/re, such that n nF(c lim -f- c? y) = n n n—>+oo n lim nFiyln) lim n ( l — e~y^n) = n—t-+oo = y. n—t-+oo Hence, lim P ( W < c + d y) = n n n n-+oo . lim P { W n n—*+oo < y/n) = 1 - e~ . y - Example 2. Let X i , . . . , X be iid random variables with common distribution function n F(x) To = \-x~ , x > l . 1 determine the limiting distribution of Z , we note that with a = 0 and b = re, we n n obtain lim nil - F(a n + b y)\ = lim re — n = -, y > 0. n and thus lim P ( Z < ny) = e x p ( - - ) . n n-++oo y For maxima, it is well-known that there are only three types of nondegenerate distributions H(x) that can be univariate limiting extreme value distributions. These results are contained in books devoted to extreme value theory, such as Galambos (1987). The three types are: H (x) = < x > 0,7 > 0, exp(—x~) lf1 0 H (x) = < otherwise, x < 0,7 < 0, exp(—(—x) i ) 2n x > 0, 1 H3,o(x) - exp(—e ) — oo < x < x + 0 0 , 7 = 0. After adding location and scale parameters the above three types of distributions can be summarized together as the generalized extreme value distribution ( G E V ) : , H(x) = exp{-[max(l + ^ ~ 7 ^ , 0 ) ] ^ } - 00 < 7 < + 0 0 . ( 2 - 4 - 1 ) The parameter 7 can be interpreted as a tail index parameter of F. It measures the thickness of the tail of F; 7 is larger if F has a thicker tail, and vice versa. The tail thickness is related to the rate of convergence of 1 — F to 0 as x —> 00. For example, let F(x) =. [max(l + 7 s , 1 - 0)] _ 1 / 7 . For 72 < 71 < 0 then — 1 / 7 1 > — 1 / 7 2 > 0 and limi/T4 ———;— = OO. [max(l + 72X, O ) ] - ^ [max(l+7 x,0)]- ^ 1 1 1 2 For 72 < 0, 71 = 0, then li m ^-1/72 7 „, ,,— = 0 0 . T, [max(l + x 7 x,0)]- /i'2 1 2 9 For 72 = 0, 71 > 0, then [max(l-f :r,0)]~ r x lim ^oo e e 1/71 7l r e x x - hm . x^oo [max(l + 7 i x , 0 ) ] / ' -x 1 > lim r—,— = oo. (x7i) /' 71 1 yi For 71 > 72 > 0, then 0 > - I / 7 1 > - I / 7 2 , I/72 — I/71 > 0 and ,. h [max(l + 7 i ^ , 0 ) ] " (1 + 7 i * ) ~ T\ ,,— = l i m — —-;— [max(l + 72Z, O)]- /^ (1 + 7 z)- /'y2 1 / 7 1 m 1 / 7 1 r 7 1 1 2 > lim i 1 + l x X ) X ' T = lim (1 + ; r ) 1 / 7 2 7 l a;-+oo M +7lX) l^ 1 X—KX> ~ = 00. 1 / 7 1 ' K So, the tail is heavier as 7 increases. For 7 > 0, from Theorem 2.1.1 of Galambos (1987), we find a n - [ ( l / r a ) " - 1] so that 1 7 = 0 and b = n 7 = 1 - F(b x) F(b x) n = [1 + ( n - l)x}7 n ~ n- x- . ih x l h and F { b x ) = (1 - F(b x)) n n = (1 - n " ^ - ^ ) " ™ exp(-a;- ). 1 n n 1/7 For 7 < 0, from Theorem 2.1.2 of Galambos (1987), we find a = — 1 / 7 , n —7 n _ 1 7 b = n so that F(x) = (1 + 7 * ) " 1 / 7 +M , n^i-x)- ^ = (-n^x)- ^ = 1 1 and F ( a + 6 x) = [1 - n - ^ - x ) - ^ - expi-i-x)- ^), n 1 n x < 0. 1 n For 7 = 0 , from Theorem 2.1.3 of Galambos (1987), we find a — logn, n so that F(x) = e~ x F(a + b x) = n ^ e * - n n and F ( a + b x) = [1 - n - V * ) " = e x p ( - e - ) . n x n n 10 b = 1 n The G E V distribution is useful for statistical inference when the distribution is unknown, and the tail index must be estimated. Since min(Xi,..., X ) - - max(-AV,..., - X ) , n n then, similarly for minima there are also three types of nondegenerate distributions as the univariate limiting extreme value distribution which can be summarized together and can be expressed in terms of (2-4-1) L(x) = = l-H(-x) 1 - exp{-[max(l + ~ ^ , 0)]- }. 1/7 <j 2.5 Multivariate Extreme Value Distributions We call a p-dimensional distribution function -P(x) nondegenerate if all of its univariate marginals are nondegenerate. One object of the multivariate extreme theory is to seek conditions on F(x), under which there are sequences {a } and {b } of vectors such that each component of {b„} n n is positive and < a + b z ) = H n ( a n + b „ z ) -»• H(z) P(Z„ n n n (2-5-1) n for a nondegenerate p-dimensional function H(z). Suppose X i , . . . , X n are independent random vectors, with common distribution function F, then H {z) n = P(Z L (w) = P(W n n l n l < z . . . , Z < z) = u n p > ,...,W >w ) Wl np p p = where G(w) = P(Xi > wu ..., Xp > wp). 11 F {z), n G (w), n Any problem on W is equivalent to one on Z by changing the basic vector X to (—X), therefore we concentrate on the vector Z of maxima. Next we list and discuss various properties which F(x.) or H(x.) may possess. P r o p e r t y 1. If F ( x ) n is a sequence of p-dimensional distributions, let the «th uni- variate marginal of F „ ( x ) be Fjf>(xi). If F (x) converges to a nondegenerate continuous n distribution function -F(x), then for each i with 1 < i < p, F ^ ( X J ) converges to the ith marginal F^(xi) of F(x). Later we shall see that all limiting distribution functions of multivariate extremes are continuous. Hence, Property 1 tells us that we can appeal to univariate case for determining the components of a and b whenever n (2-5-1) n holds. The following several properties are important in multivariate extreme value theory. P r o p e r t y 2. Let X i , . . . , X be iid p-dimensional vectors with common distribution n function F(x.). Then there are vectors a and b > 0 such that ( Z n n to nondegenerate distribution function H(x.), n —a n )/b n converges if and only if, each marginal belongs to the G E V family, and if C [y\'\...,y J ] n l F n -> CW(vi,• • •\y ), n^oo. P Property 2 tells us that given a p-variate distribution function .F (x), we can check n if its marginals belong to the G E V family, if so, then we use the methods of the univariate case to determine the components of the vectors a and b ; furthermore we can n n determine CV(y) by its definition and check if C ^ y / " ) converges. 1 Example 3. Let (X, Y ) have a bivariate exponential distribution F ( x , y ) = l-e- -ex y + G(x,y), where G(x,y) = P ( X > x , Y > y ) . 12 (2-5-2) If ( Z - a ) / b n n converges to a distribution H(zi,z2), n (logn,logn), and b then we can choose a n = (1,1), and obtain n e~ + e~ Zl Z2 F ( l o g n + 2i,logn + z ) = 1 h G(logn + zi,lo'gn + 2 z ). 2 n For the Marshall-Olkin distribution G(x,y) = exp[-rc — y — Xmax(x,y)], The last term of (2-5-2) is exp{ — (2 + A)logn — z — z — A max(zi, 2 )} x 1 = 2 2 A > 0. ex'p{—z\ — z — \ 2 + A ri = ma,x(z ,z 2 1 1 0(n-"- ), A and lim F ( l o g n + zi,logn + 2) n 2 = lim n—>oo e 1 - 2 i _(. -^2 \ e TI n—>oo V = For / exp(—e * — e - 1 - 2 2 71 I ). \ another example, let G{x,y) = (e 6x + e - l) 6y _ 1 / e , 6> > 0. Then G(log n + log n + z2) = (nV „ 2 1 +nV I( ^i n 2 2 l)~ - 1/e )"!/«, e and lim F (logn + ,logn + * ) n 2 l 2 = lim(l-^—^ = {exp[-(e + e )]}{exp[( ^ + e* )" /*]}. 21 - K ^—L 22 22 e 13 )» 1 2 In the second example, there is dependence in the limiting extremes. Property 3. A p-variate continuous distribution function H ( x ) is a limit distribution in (2-5-1) if, and only if, its univariate marginals belong to the G E V family and if its copula CH satisfies C {y\'\ k H ..., yl ) = C (yi, lk H • • •, yP), k = 1,2, • • •. (2-5-3) Property 3 tells us that if the limit H ( x ) exists, then we check condition (2-5-3). If it holds and if C'H is a copula, then we have got the actual limit distribution. Example 4. Let H\(x),..., H ( x ) belong to the GEV p H(x) = family, then H {x)..-H {x) x p is a possible limit of (2-5-1). This limit consists of independent univariate marginals. By assumption, the condition of Property 2 on the marginals is satisfied. Furthermore, by definition, CW(y) = y\ • • • yP for which the condition of (2-5-3) is evident. A n appeal to Property 3 yields the claim. Example 5. The distribution function H(x!, ...,x ) p = exp{exp[- m i n ( x i , . . . , x )}} p is a limit in (2-5-1). The marginal distributions are Hi(x{) = exp(—e~ ) = ^ ^ ( x , ) . Therefore it remains Xi to check the validity of the condition (2-5-3): exp{— exp[— m i n ( x i , . . . , x )]} p = min{exp[— exp(—x,-)] and C (yu H • • •, y ) = m i n ( j / i , . . . , y ). P p 14 1 < i < p}, Hence for k > 1, C (yl , k yl ) ,k = [mm(yl , lk H vl )] /k /k = m i n ( , . . . , y ). k p y i Example 6. The distribution H ( x x ) = # ,o(si)#3,o(>2)[l + ^(1 - #3,of>i))(l - #3,o(z ))] u 2 3 2 does not occur as a limit in (2-5-1). (See Galambos(1987).) Even though the marginals H\(xi) =H 3 f i ( x ) and the copula i / 2 ^ 2 ) = -#3,0(22) 1 2/2) = ym[l C {yi, H - 2/2)] + ^ ( l - yi)(l fails to satisfy (2-5-3). Let X = (Xi,... be a vector with distribution function F ( x ) . Let j(k) iX ) p = (ji, • • • ,jk), 1 < k < p be a vector with components 1 < j i < J2 • • • < i t < P- The distribution function Fj(k)( h, • • • ^j*) is x is obtained from F(x) a A;-dimensional marginal distribution, which by letting X{ —> +oo for all % £ {ji, • • • ,jk}- Let Gj{k){ JIT ••i jk) X x = F (-^ji > x j n • • • •>Xj k > Xj ). k Assume that F ( x ) is such that each of its univariate marginals belongs to the G E V family. Then, there are a, and 6, > 0 such that n F ( a , „ + binx) lim where Hi(x) n n t = Hi(x), 1 < i < p, (2-5-4) belongs to the G E V family. We assume that a, and 6 n determined and put a P r o p e r t y 4. n = ( a i „ , . . ., a ) and b pn (Z — a )/b n n n n = (b ,.. tn > 0 have been ., b ). ln pn converges to a nondegenerate distribution i / ( x ) if and only if for each fixed vector j(k) and x for which Hi(xi), 1 < i < p, of (2-5-4) are positive, the limit nHS> Gj(k){ jin n a ~i~ bj Xj , in in . . . , Oj n fc "H bj n.Xj ) k 15 kn — hj (Xj k 1 , . . . , £j ) t (2 5 5) are finite, and the function H(x) = exp{£(-l) h (x ,... )} £ f c ih h (2-5-6) lXjk 1 < J 1 < J 2 — <3k<P fc=l is a nondegenerate distribution function. If the actual limit distribution of ( Z — a ) / b n n n is the one given in (2-5-6), then the following inequalities hold. Let s > 0 be an integer, then #(x;2s + l) < H ( x ) < H(x;2s), (2-5-7) where i7(x, 0) = 1 and #(x,r) = e x p { £ ( - l ) h (x ,... £ f c jk When s = 0 in (2-5-7), the inequality H ( x , ,x )}. h (2-5-8) Jk 1) < #(x) < 1 obtains. This means that iJ(x) is never exceeded by the product of its univariate marginals, and H ( x ) exceeds 1. Example 7. . Let F(x, y) = 1 — e~ — e~ + G(x, x G{x,y) = (e y 6x +e y) where -I)- ' . ey 1 0 From the univariate case we know that lim i*7(logra + x) = exv(-e~ ), x > 0, x lim F "(log n + t/) = exp(—e ), _2/ 2 Since G\(x) = e - x and G ( y ) 2 = e _ y War) = lim nG(logra + x) = lim n—»oo n—*oo Similarly, h ( y ) 2 y > 0. ne~^ n + x ) = e~ . x = e , and _J/ h (x,y) 12 = Jim nG(logn + x,logn + 16 y) never lim n(nee6x = +n e e n—*oo n—+00 — l)~ 1/,fl JI<> (e * + e )- ' . e = 9 y 0y 1 e Therefore, the limiting extreme value distribution is H(x,y) = exp{-h (x)-h (y) 1 = + 2 h (x,y)} 12 e x p { - e - - e~y + (e8x + x e )~ }. By 1/e The following definition is needed for the next property. Definition: The p-dimensional unit simplex S is the set of vectors q with nonnegative components qi such that YA=I Qi = 1Property 5 (The Pickands Representation for a min-stable exponential distribution). G(xi,..., x ), with univariate exponential margins, is a survival function satisfyp ing — \ogG(txi,... ,tx ) p = —t\ogG(x ,... 1 ,x ) p for all t > 0 if and only if G has the representation -logG?(a;i,...,a;p) = / [max. (q x )]dU(q ,... i Js p where S p — {(qi,---,q ) P i 1 ,q ), p x, > 0,i = 1,... ,p, (2-5-9) 1<»<P ' qi > 0 , i = 1,... ,p, Yli qi — 1} is the p-dimensional unit simplex and U is a finite measure on S . p Property 5 is a result of Pickands and an alternative statement of Theorem 5.4.5 of Galambos (1987). It is not an easy task to give the representation (2-5-9) for a given H(x), but one usually does not aim at giving another form of H(x) when it is already known. The value of (2-5-9) lies in its possibility of generating functions which are limits in (2-5-1). For applications of Property 5, see Joe (1990a). A n example from there is: exp{-A( , Zl . . . , z ; \ , B c ( l , . . . , p ) , 6 ) } = e x p { - £ A ( £ x ) ' },6 6 p s B t l s > 1. (2-5-10) 17 For each integer m > 2, it can be directly shown that x m ) = e x p { - ( ^ + --- + ^) / } 1 (2-5-11) 5 is a survival function over Xj > 0,j = 1,... ,m,for 1 < 8 < co. (2-5-11) satisfies the homogeneity of degree one condition of Property 5. Hence by Property 5 , if m < p, there is a measure V m on S m and a measure U on S such that m p where Um is a representation of Vm in p dimensions which puts all mass on qi,... ,q, Therefore, £ M £ ^ ) B^0 1 / f i ieB has the representation where dU has the form J2 B \BdUs- B y Property 5, (2-5-10) is a distribution. 2.6 Some relation between extreme value theory and Central Limit theory Since Central Limit theory is more familiar to statisticians, in this section we show some similarities and differences between Central Limit theory and Extreme Value theory. Consider the univariate case first. Suppose that X\,... ,Xn are iid with distribution F. Let Z n W = mm{Xi,.. = max{Xi,...,X }, n n 18 .,X }. n We would like to know whether there are any sequences { a } n S —a n — A Z- a n n ' and {b } such that n W — a 1 n n A '. ° r n converge in distribution. From Central Limit theory, if F has a finite second moment, then there exist and {bn} such that the limiting distribution (Sn — an)/bn {a } n has a normal distribution. Definition. Let H ( x ) be a nondegenerate distribution function, then F(x) domain of attraction of H ( x ) if there are sequences { a } is in the and {b } > 0 such that n n lim F ( a + b x) = H(x). n n n From extreme value theory, we know that the limiting distribution H ( x ) of ( Z — n o-n)/b n belongs to the GEV attraction of For family if the distribution function F(x) H(x). the multivariate case, assume that X i , . . . , X X; = (Xn,..., Z Are is in the domain of n n are iid with distribution F . Let S =Xi+••• + X , Xip), n = ( m a x ( X n , . . .,X ),... nl n ,max(A"i ,.. .,X )). p np there any {a } and {b } > 0 such that n n converge in distribution? If the limit of (S — a ) / b and the second order moments of n n n F exist, then from Central Limit theory the limit of (S — a ) / b has a multinormal n distribution. 19 n n If the limit of ( Z — a ) / b exists, from extreme value theory, each univariate margin n n n must be in the G E V family. A n analogy between multinormal and multivariate extreme value distribution is: (1) Suppose ( Z i , . . . , Z ) has a multinormal distribution, then linear combinations p of univariate normal, i.e. a\Z\ + • • • + a Z has a univariate normal distribution. p p (2) Suppose Xj ~ H(-,~fj, fj,j,o-j). After transformation to an exponential survival function Zj = - ^og{H(Xf^ n aj)) h h have an exponential distribution. If (Z\,..., Z ) has a min-stable multivariate exponenp tial distribution, then from Property 5, G(zi,...,z ) p = exp{-A(z ,...,z )} 1 (2-6-1) p and V =. m i n ( ^ - , . . . , ^f-) has an exponential distribution - ( i>-> p) ^ where Wj > 0. vA w w e That is, weighted minima are univariate exponential. Note also the main difference: The multinormal distributions form a finite-dimensional parametric family, the multivariate extreme value distributions form an infinite-dimensional family (as given by Pickand's representation). Since parametric inference is easier than nonparametric inference, a "goal" is to find good parametric (finite dimensional) subfamilies of the infinite-dimensional family (see the next section). 2.7 Parametric families of multivariate extreme value distributions. A multivariate extreme value distribution with margins transformed to a survival function with exponential survival functions as univariate margin is a min-stable exponential distribution. Let G be a p-dimensional min-stable exponential distribution, from extreme value theory A = — log(G) satisfies A(tzi,...,tz ) p = tA(zi,...,z ), p 20 Vt > 0. Let G be the marginal survival function of G, then G ( z ) = e x p { — A ( z ) } where A S s a s s S is obtained from A by setting Zj — 0 for j $ s and s is a subset of {1,... ,p}. Here we will list some parametric families of multivariate min-stable exponential distributions in terms of A . These families can be found in Joe (1990a,b), Tawn (1989), Coles and Tawn (1991), Husler and Reiss (1989). The interpretations of the parameters are given after all the families are listed. The families are given here for p = 2,3,4, from which the general multidimensional form can be seen. A(z z ;S) A( , s 6 1 r 2 x 2 r 2 )~ x l l 3 (2-7-1) r > 0, \ + zi*) l \ A ( z z2, z3; 8, 8 ) = ((z* + zs2)s*'s u 1 s 2 z ; r) = z + z - (z7 + z ZL 8>1, = (z + z ) ' ,- u 2 s 6 > 8 (2-7-2) (2-7-3) > 1 , 3 A(ZI,Z ,Z ;T,T ) 2 = 3 3 z + z + z - (z~ + z ) T 1 2 ( ( ^ + ^ T T ) T 2 3 1 / T 2 3 / + ^3 T A ( z z , z , z ; 8,8 ,8 ) u T 3 4 3 T 3 )" 1 / T 3 - (z7 5 1 1 2 r>r >0, , (2-7-4) 3 = ([(z + z ) *l 6 s 6 2 4 + z ^ ) - ^ - (z ^ + z ^ ) ' ^ + T3 + zi*} *' > + zi*) ' *, s s 1 8 l < 8 4 < 8 3 < 8 , (2-7-5) A ( z , z , z , z ; 8 , 8 , 8 ) = {{z[ + zs2)Si'& x 2 3 2 4 4 + (z > + z ^ * ) ' * , s 3 6 1 6 1< 6 ,8 2 4 < 8, (2-7-6) A(zi,z ,z , 2 3 z; 4 T,r ,T ) 3 4 21 = z + z + z + z - (z7 + z r 1 2 (z7 3 +z7 )- Ti T i 3/T {{ r +*3 r T T3 Z (((*r + ^ ) T {z * - 1 / T i (c*r + ^ T 3 T 4 T 4 4 +z " ) T3 T3 )" 1 / T 4 T 3 ~ +^" 1 / T 3 T 4 ) -(z T + ((*r +zrr 1/T4 3 + ^" )" T4/T3 T4 3 1/T4 , 4 r i/T3 +z7 )~ T 3 T 3 2 + 1/T4 +*r )~ 4 1/T4 r > r > r > 0. 2 - 1 / T 3 + _ 1 / T 4 + (c*r+*-2 r' +*r )~ 1/T3 +*r )~ - (z~ + z ^ ) - 1 / T 4 T3 T 3 / T )~ +2 " T 2 +^- )" /T3 2 (2-7-7) 4 A(ZI,Z ,Z ,Z ;T,T ,T ) 2 3 4 2 4 = z + z + z + z - (z7 + 2 - T T 1 2 3 (zp + z r ) 4 (c*r + * n 4 -T2 3 O r +* "T 2 - (sp + z - T T T4 1/T4 2 +*3- )" T 4 + ^r ) 2 4 / T T4/T2 1 / T 4 - +* r r + ((*r +* - ) T T 4 / T 2 )" + (*r + (* 1/T4 4 T 2 T 2 3 _T2 3 + (*" + - ) - ( ^ +^ 1/T4 3 - (^ " +*rT _ 1 / T 4 T 4 / T (*r + (* 1 / T 2 4 T 4 / T 2 2 4 ) _ 1 / T 4 1 / r i +*r )~ 4 + *r ) 2 2 4 1 / T 4 + )~ 1 < r , r < r. , )- + 1 / T 4 T4/T2 T 4 1/T4 - (2-7-8) Finally, the Hiisler and Reiss(1989) model with p = 2 is A(zx, A) = *(A + l 0 g Z l -2A l 0 g Z 2 ), + *(A + ° 1 2 g Z 2 -2A l 0 g ^), , 1 A > 0. (2 - 7 - For p > 2, let Xij € (0,oo) VI < i , j < p with i ^ Put Moreover , for 2 < k < p and m = ( m . . . , m^) with 1 < m i < m < . . . < 1 ? 2 < define rfc,m — ^ ( A ^ . ^ + x . — \ .^ .) i,j<k-\ 2 2 2 m mk m 22 m (2-7-H Furthermore , let denote the survivor function of a normal random vector with S'(-|r) mean vector 0 and covariance matrix T. The extension of (2-7-4) is A(z ,...,z ;\ j, 1 p i \ < i < j <p)= p k J2(~ ) '52 k,m(-1ogz ,...,-\ogz ) 1 k h mi k=\ mk (2-7-11) with , /•oo ^m(y) = / S -, (y- - Z + 2 A f 2 jf- |r ,m 1 n > ) m f c 6 ~^> Vk J for 2 < k < p, where J2 nieans summation over all p-vectors m = ( m i , . . . , m^) with k 1 < mi < m < ... < m 2 < p, and /ii, (y) = e~ for m = 1,... ,p. The cases p = 3,4 y m k are given explicitly in Chapter 3. For comparison, note that the bivariate models (2-7-1), (2-7-2) and (2-7-9) each have a single dependence parameter. For (2-7-1) and (2-7-2) dependence increases as the parameter increases; for (2-7-9) dependence decreases as the parameter increases. The trivariate models (2-7-3) and (2-7-4) have respectively models (2-7-1) and (2-7-2) for each of the three bivariate margins. For the trivariate model (2-7-3) (respectively (2-7-4)), 8 (respectively T) is the dependence parameter for the (1,2) bivariate margin; £3 (respectively, T3) is the dependence parameter for both the (1,3) and (2,3) bivariate margins. We have something similar for the 4-variate models (2-7-5) to (2-7-8). The 4-variate models (2-7-5) to (2-7-8) also have respectively models (2-7-1) and (2-7-2) for each of the six bivariate margins. For (2-7-5), 8 is the dependence parameters for the (1,2) bivariate margin; <$ is the dependence parameter for the (1,3) and (2,3) bivariate 3 margins; 84 is the dependence parameter for the (1,4), (2,4) and (3,4) margins. For (27-6), 8 is the dependence parameter for the (1,2) bivariate margin; 82 is the dependence parameter for the (3,4) bivariate margin; 84 is the dependence parameter for the (1,3), (1,4), (2,3) and (2,4) margins. The model (2-7-11) has a dependence parameter for each of the p(p— l ) / 2 bivariate margins. For (2-7-11), the dependence parameter of the 23 bivariate margin, with i < j, is A,j, and the margin has form (2-7-9). Models (2-7-3) to (2-7-8) and their extensions have p — 1 parameters in total; as above, some of the bivariate margins have the same dependence parameter. Hence these models do not have as much flexibilty in the dependence pattern but they have a simpler form than (2-7-11). 2.8 Estimation for a parametric family of multivariate extreme value distribution. Joe (1990b) suggests the following procedures for fitting of a parametric multivariate distribution to iid X i , . . . , X . n (1) Fit the p univariate margins separately by maximum likelihood using the G E V family. (2) Transform each margin so that each transformed variable has exponential survival function Gi(zi) = e~ . Zi (3) Suppose that the transformed data are iid p-vectors, then compare fit of different families of min-stable multivariate exponential distributions. Check for parameter estimation consistency with bivariate and higher order margins. (4) Go back to the original data and estimate both multivariate parameters and parameters of the univariate margins simultaneously, using maximum likelihood if a copula has been decided on as being acceptable. In (3), to check the fit of a model, we compare parametric estimates of A with a nonparametric estimate of A. if (Zi,..., Z) p Besides the consistency checking from we know that is a random vector with the min-stable exponential survival function G — e~ , then A E(wm[^,...,^]) Wi for = w [A{w ,;..,w )]- 1 1 p p u>j > 0, j = 1,... , p. So if Y i , . . . , Y p 24 are a transformed random sample, expo- nential probability plots of Vi = minjYij/wj can be used to check the assumption of min-stability for multivariate exponential distribution. If the plots are adequate, a nonparametric estimate of A(wi,..., w ) is given by n/ p parametric estimates of A from different models. 25 J27=i ^ which can be compared with Chapter 3 Maximum Likelihood Estimation With The Husler-Reiss Model Hiisler and Reiss (1989) obtained a new class of multivariate extreme value distributions by taking a non-standard extreme-value limit. Let X and M n be as defined in Chapter 2. For a mutivariate normal random vec- tor, having all correlation coefficients smaller than one, Sibuya (1960) proved that the marginal maxima M i , . . . , M „ n > iP are asymptotically independent. The rate of con- vergence to asymptotic independence is slower as the correlation coefficients increase. Since in practice the sample size n never goes to infinity but is a fixed positve integer, is there another asymptotic formulation which may provide a better approximation? This motivated the work by Hiisler and Reiss. For the bivariate case, Hiisler and Reiss let the correlation coefficient p — p increase n towards 1 as the sample size increases. It is shown that the marginal maxima are neither asymptotically independent nor completely dependent if (1 — p )\ogn n converges to a positive constant as n —> oo. This was extended to dimensions p > 2. Hiisler and Reiss define their limiting distribution by the c.d.f. or survival function. To get a maximum likelihood estimation, the density is needed. To obtain the density, the following result for the conditional distribution of a multidimensional normal distribution is needed. Result 1: Let U = (Ui,..., U ) and V = (Vi,..., H ) be random vectors such that a MVN 12 0, V ^ £21 ^22 \ JI where E n , £12, £21, £22 are a x a, a x b,b x a,b x b matrices respectively, and suppose £ 1 1 and £22 are nonsingular. Then U | V = v ~ MVN ( EjaE-iy; E 26 n - EiaE^i )• Let F x and F x denote the distribution function and survival function for a random vector X , and let fx and /x|Y denote the density function of X and the conditional density function of X when Y is fixed. Then the joint distribution of X and Y is /x,y( »y) = /x|Y(x|y)/v(y)x Let <^(x;E), $(x; E ) denote the multivariate normal density and survival function with zero mean vector and covariance matrix E . Then /u,v(u,v) = ^(u - S H S ^ V J S H - S12E- E21)^(v;S22) 1 and Fu,v(u,v) = $((u,v);E) J roo f "1 ' /-co roo •••/ / Jvi, Jui too •••/ 0(s-Ei S - t;Eii-S Sj E )^(t;E )dsdt 1 2 Ju 2 1 2 1 2 2 2 1 2 2 a •••/ ^ ( u - E i a E ^ t j E u - E w E ^ E a O ^ E M j d t . Therefore, 3 <&((u,v);£) _ dv\ • • • dvb b ( _ 1 ) 6 $ ( u _ S i 2 S - i v ; E i i _ E i 2 E - i S 2 i ) ^ ( v ; ^ Now we are ready to derive the density for the Hiisler and Reiss survival function when the univariate margins exponential have been transformed to distributions with mean 1. We provide details for dimensions p = 2,3,4 from which the general case will be apparent. The notation gets increasingly difficult as p increases. 3.1 B i v a r i a t e Case Frome Section 2.7, G(wi,w ,X) 2 = e x p { - A ( i 0 i , u ; , A)}, 2 27 2 where A) = $(A + ^ \og(wi/w ))w A (w ,w , 2 1 2 2 Since J^- = — f ^ - G and J^- = — dG 2 = _d_ dw\dw 2 a = 1 2 x2 -\og[cf)(\+ dA dw 2 From 4>(x) = (2Tr)~ ^ e~ ^ A > 0. 2 we obtained the density function of t7(u>i, w , dwi dw 2 + $(A + ^ log(w /u;i))u)i, 2 1 2 &4 g A 2 = dw 2 c dw\ dw 2 _ d A 2 2 dw dw 2 x 2 we have ^'(x) = a;^(x). Also 2 ^\og(w /w ))/w ] 2 1 1 = A / 2 + (logto - logwi)/2 + ^ - ( ( l o g w - l o g ^ i ) -t-logiu! + log V^TT = A / 2 + (log w = - l o g t y ( A + ^-log(tOi/u; ))/u; ]. 2 2 2 2 2 - log w )/2 + g^-((log tux - log w ) 2 x 2 2 2 + log w 2 + log V2TT 2 Hence <KA + ^ log(u; /u)i))/u;i = <^(A + ^ log^/u^))/^. 2 Therefore dA 2 d 1 1 d w ^ ($(A + 2A — l o g ^ '/ u ^") ) ^ + $(A +2A — log^/w^w-^) v = ~ ^ ( A toV + 2Xlog(w2M))(u> M) + ^ ^ ( +$(A + = v ^(X + 2A 2 ^\og( /w )) Wl 2 ^-\og( /w )). Wl 2 Similarly, 3\A. 1 Finally, 28 A + ^-log(^iM))(^i/^2) d o^ m X + 2 ^ ( A From Ait2, 1 2X l o g { w i / w 2 ) ) + ^logKM)). A ,i and A\ we can get the density function g(w ,w , 2 2 1 A). 2 3.2 T r i v a r i a t e Case Let f E =2 2 \ \2 _i_ \2 ^13 ' ^23 "12 2A^3 / 9\ 13 \2 I \2 \2 13 "t" ^23 12 2 y A — A be the matrix in (2-7-10) and let \2 _ ^13 + ^23 ~~ \2 A P12 = P21 - P13 = psi = 2A13A23 1 \2 , \2 12 + 23 23 1 2Aj A23 \2 . \2 12 + 13 23 A A A A A 2 P23 = P32 = — A 2A12A13 i.e., Pij = 2A,/;Ajfc for i,j, k distinct. Then h 3 < 1 2 3 in (2-7-11) can be written as C(wi,w ,w , 2 3 A12, A 1 3 , A 2 3 ) = / JO $ ( A + — — ( l o g g - l o g t u i ) , A 23+777—(log g-log w2); LA ^ 13 13 A 23 where with some abuse of notation 1 1 P 1 $(•;/>)=¥(•; \ P 1 From Section 2.7, we have G(w ,w ,w ) 1 2 3 - exp{-A (w ,w ,w , 3 1 2 29 3 X , X , A )} 12 13 23 p )dq 12 where with simplification of (2-7-11), A (w ,w ,w ,Xi ,Xi ,\ ) 3 1 2 3 2 3 ~w\ - w - w + A (wx,w , 23 2 3 1 Ai ) + 2 2 A13) + A (w ,w , A (w ,w , 2 2 3 2 2 X )+ 3 C(wi,w ,w ). 23 2 3 Constraints are 0 < A < 0 0 and 0 < pij < 1. tJ By symmetry, equivalent forms for C can be obtained by permutations of subscript indices: C(iui, w , w , A , A i , A ) 2 3 rw\ = / 12 3 23 1 ^(Ais + 77T-(log q - log w ), A 3 23 1 $(^12 + 7 T T ( S 9 - logioi), A - 0 2 1 ru>2 J' 1 + -r— (log q - log w ); p )dq 12 + — — (log 9 - log w );p )dq. lo 6Ai2 The density function for G(wi,w ,w ) 2 2 3 3 is: 3 13 4A23 d 3 g(w ,w ,w ) 1 2 3 = ( - l ) ^ — 7 ; — « — e x p { - A ( w , w , w , A i , A i , A )} 3 3 OWiOW OW _ dA d A dw\dw dw 2 _ From SMdMdM dwidw dw 2 2 1 2 3 2 3 23 3 _ 8A d A dw dw\dw 2 3 2 3 2 3 3 3 2 3 _ dA d A dw dw\dw dA dw dw dw 2 3 3 3 3 + 3 2 1 <j>(xi,x ; p) = (f>( )(f)(^=^) 2 Xl 00 f°° ±t i ~ P ~L z Z 1 and ^ , x ; p ) = 1 2 $ ( 7 f = ) ^ i ) . Let Cl,23 := = $(Ai3 + atvi \og(Wi/w ), A12 + J-J—log(lWi/u7 );P23)zA 3 zAi3 2 1 2 30 2 3 B y symmetry, we have C ,i3 •- dC and C ,i2 := 2 Similarly, 3 dC 2 12,3 dw\dw<i $[(( i3 + 2^1og(u7i/u;3)) A - ^ (Ai2 + ^ M t o i / w j ) ) ) / ^ 12 - P23} 2 3 x <K i2 + ^—log(u;i/u;2)); A 2Ai2io B y symmetry, we also have C 2 := g^f^ and C , i := afffe"- Finally, 23 1 3 i 2 dC 3 C i o :— 3 dwidw dwdw 2 3 3 1 <K 13 + ^T—log(w M),A 3 + A 3 £A\ ^T—^Og(w /w2);pl )jT-^ 6A23 2 3 3 <K i3 + r r 2 ^^13 Combining C i ±A\ A2 W\ 3 3 t«2 log(^3/wi))<^[((^23 + 7 7 — l o g ( w / w ) ) - A 3 ^^23 13 log(wi/w3)))/y/l + ^77 Pi2(Ai3 2 - p\ ] X 2 4A A 1 3 ^2 \ / l ~ ^ - P i 2 1 23 , C ,i3, C7,i2, Ci , , Ci , , C ,i, Ci23,and f^ -, f ^ f ^ ^ g ^ , 2 li23 2 3 9~~9^7> and QJ^Q^£Q w 2 3 3 2 23 a£rfc have an expression for the density function <7(it>i, iw> u> ). e 2 W 3 3.3 Four-Dimensional Case Let _ \2 + A2 24 12 \2 1 \2 _ \2 \2 1 \2 _ \ 2 14 + 24 12 14 + 34 13 \2 1 \2 _ \2 2A 24 + 34 23 -r1 A\2 \2 A 2A? E = 2 A y 2 4 \2 A 14 A 4 A 34 \2 —A 13 A A A A 24 24 A 1A \2 \2 —A 34 23 A A 2A| A 4 be the matrix in (2-7-10) and let Pijk A \2 1 \2 _ \2 t j + ifc jk ij ik A A ZA A where A,j = A - and, «, j , A; are distinct and between 1 and 4. Then / i i Jt 4 ) can be written as: 31 2 3 4 in (2-7-11) D(w ,W ,W ,W ) 1 2 3 = 4 rm l / $( 14 + 7T\ ^g(q/w ), l A + p-r— A 24 4 JO &A\ A34 + 7 7 7 — ZA34 log(q/w ), 4 ^^24 4 \og(qjw ), R(p i2, 4 P423))dq, 0413, 4 where f 1 a 0A R(a, b, c) — X a l e 6 c 1 By symmetry, equivalent forms for D are D(lWl,l0 ,f>3,tU4) 2 y^i = / JO 1 $(Ai2 + 7 7 7 — log(?/w ), A 2 A^Vl A14 + rw2 2A 14 Jo 1 + 7 7 7 — log(o/w ), 3 ZA13 log(a/to ), R(p 3, 4 0124, 12 1 *(Ai2 + 7 7 7 — log(o/wi),A / a 3 23 P\3A))dq 1 + 777—log(?/w ), 3 £>A\2 ^ 23 A •^—\og(q/w ),R(p i3,p2i4,p234))dq A4 + 2 2 4 ZA2 4 / ru> Jo 3 \ \ • $(Ai3 + 7 7 7 — log(g/wi), A + ^ 7 — log(?/u;2), AA\ ZA23 A34 + 77T— I°g(?/W4), #(0312, 0314, P324))dq. 2 3 3 ZA34 Constraints are: 0 < Xij < 0 0 , i,j, i ^ j, 0 < /o.-jjk < 1, are distinct and R(P412, 0413, 0423), R(P\23i 0124, 0134), #(0213, 0214, 023 ), #(0312, 0314, 0324) 4 are positive definite. 32 Taking derivatives, we have dD -Dl,234 —.. o^r* 1 , 2^ { X i 2 + i , . ' i ^g(w /w ), TT—log(«;i/u;3),Ai4 + 7 7 7 — A13 + . l o g ( W i / W 2 ) 1 R(p , 4 123 0124,0134)) ZA13 By interchange of 1 and 2, 1 and 3, 1 and 4 in 1)1,234 we can obtain As,i24 : = | £ , -04,123 := f j - B y making use of Result 1 with E 3u>3 : / 1)2,134 1 := g ~ 0134 n ^ 0134 1 j \ 0123 J12 , E22 = 1 we have \ 0124 j d D dwidw 2 D 12,34 • — = 2 dD 234 dw lt 2 — 1 $( 13 + 7 7 T — A ZA13 log(wi/w ) 3 1 ~ 0123(Al2 1 + 7 - — \ 0 g ( w 1 / w 2 ) ) , \ l 4 + 777—log(tt>l/w 4 ) ^Ai2 ^ 1 i 4 1 1 ^°&( l /l)) 7y[ 012 (Al2 + 7 7 7 — log(wi/w2)),#34,12)<KAl2 + 77\ W 4 ZA12 A . W ^Ai2 ZA12W2 where 1 R.34,12 0134 ^ 0134 ( I \ 0123 0123 0124 1 j 1 ^ 0134 \ 0123 — — 0134 01230124 — 1 01230124 \ 2 - 0124 ) By permutation of indices we can easily get D 13,24 •— Q g. W 9D 9D 2 23 9wi dun '-^23,14 dw dw i 2 3 r) 24,13 •— iy 2 dw^dvn ' Again using Result 1 with E n — 1, / S12 — \ 0123 / 1 0134 0134 1 , E22 — 0124 ) I 33 \ / dD 2 9 ^ 9 ^ 5 -^34,12 - 5 then y-i ^22 — -#123,4 = 1 — E i S 2 1 1 1 - Pl34 ^ —Pl34 E i = 1 — (p 2 2 -/>134 2 123 +p 1 2 4 1 j — 2/)i 3/)i /) 34)/(l — p 2 24 1 1 3 4 ) and d D 3 D 123,4 dw\dw dw 2 3 ${(Al4 + -^—logiwxIWi) - [(/9123 - />124Pl34)(Al2 + T"J— l o g ^ i / ^ ) ) + Z/\i4 ^^12 f>124 - /9l23/»134)(Al3 + ^ 7 log(li>i/u> ))]/(l ~ p\ ))I\J ™A) R 3 34 1 </>(Ai2 + -J—log(t«i/ti;2),Ai3 + r-J— l o g ^ / ^ a ) ; ^123)7. . ZA12 ZA13 4Ai2Ai3U;2^3 By symmetry and permution of indices we can easily get dwidufySwi D' £*i24,3 Finally, 134,2 • dw^dw^dun1 -^234,1 d D A D 1234 dwidw dwzdw 2 4 •^—log(w /w ),X <^(Ai4 + 1 4 + -^—\og(w /w ),X 24 2 ZA\ 4 + -^—\og(w /w ),) 3 34 ZA24 4 4 ^A 4 3 R{,P412, / > 4 1 , / > 4 2 ) ) g ^ ^ ^ ^ ^ ^ 3 3 From Sections 2.7 we have G(w!,w ,w , 2 w ) = exp[—A (iwi, w , w , w , X , X , A , A , A , A , A )], 3 4 2 4 3 4 12 13 i 4 2 3 23 2 4 34 where w , A , A , A , A 3, A 3, A , A ) = A (w ,w ,w , 4 1 2 3 12 4 13 14 2 2 24 34 ^i(iwi) + h (w ) + hi(w ) +fei(iw )- h {w ,w ) - h^W!, w ) - h (w ,w ) 1 2 3 4 2 1 2 3 2 x h (w ,w ) - h (w ,w ) - h (w ,w ) + h (u)!,w ,w ) + h (w!,w ,w ) + 2 2 3 2 2 4 2 3 4 3 2 h (u)!,w ,w ) + h (w ,w ,w ) + h (w!,w ,w ,w ) 3 3 4 3 2 3 4 4 2 34 3 4 3 3 2 4 4 = IWi + w + w + w - [wi +w - A (w ,w , 2 3 4 2 2 1 A )] 12 2 [W!,W - A ( l U i , to , Ai )] - [wi + w - A (wi, 3 [w + w 2 3 2 3 2 2 3 23 [w + w - A (w ,w , 4 2 3 2 4 3 4 2 2 34 2 3 24 3 4 2 1 2 13 X , A ,A 12 23 14 A , A , A ) - D(u>i, 1 ^ 2 , ^ 3 , ^ 4 ) 4 23 = -2(wi + u> + w + w ) + A (w ,w ,\ )] A (w!,w , 4 2 C(w ,w ,w , i3 2 2 A )] + C(ioi,tw ,iW3, Ai , A , A ) + C(wi,w2,w4, 4 C(wi,ti> ,t«4,A ,A14, A 3 4 ) + 2 W , A14)] - 2 A (w , w , A )] - [w + w - A (w , w , A )] - 3 3 4 24 34 + A (wi,w ,\ )] 12 2 3 + 13 Ai )] + A (w , iw , A )] + A (u> , u>, A )] + 4 4 2 2 3 23 2 2 4 24 A (u; ,u;4, A )] + C(w , w , w , A , A 1 3 , A ) + C(u>i, w , w , A i , A , A ) + 2 3 34 1 2 23 C(w ,w ,w , C(iwi,io ,u;4, A 1 3 , A i 4 , A ) + 3 12 3 2 3 4 3 2 4 2 i 4 24 A , A 4 , A 3 4 ) - D(u>i, w2, w3, w4) 4 2 23 and g(w , ,w ) uW2 W3 = 4 ( - l ) aA aA <9A d A 4 4 4 ^ 4 F A a A i 8 ^ ^ ^ 3 , ^ ) 4 aA4a^4 2 4 > 4 dw dw dw dwidw dw dw dA dA dA dA 8A dA dA 8A dA dw\dw dw dw dw\dw dw dw dw dw dw-^ dw dA 8A dA dA 8A dA dA dA dw dw dw\ dw dw dw dw\ dw dw\dw dw dw dA dA dA dA dA dA dw\dw dw dw dw\dw dw dw dw\dw dw dw dA dA dA dA dA dA dwidw dw dw dw dw dw dw dw dw dw dw dA dw\dw dw dw 2 3 4 2 2 3 4 2 4 4 3 2 4 2 4 4 4 4 2 4 2 4 3 2 2 4 2 4 4 3 2 3 4 2 4 4 2 2 4 3 2 3 4 4 4 3 3 | 4 4 3 4 3 4 2 4 3 4 4 3 4 2 | 2 4 2 4 4 2 4 3 4 4 3 2 4 4 2 4 1 3 4 4 4 2 2 3 4 4 x A 4 2 3 4 Upon differentiation, dA 4 -2 + A 1>2 + A 1<3 + A 1A + Ci,23 + Ci,24 + Ci,34 - D i ,234- By symmetry,!^ ,1^-, f^ are like f^- with indices 1 and 2, 1 and 3, 1 and 4 inter1 1 changed. Also, Pi A. 1 dw dw 1 2 A\2 + Cl2,3 + Cl2,4 35 — ^12,34 24 and o Z — W23 ~ -L/123,4 5 OW1OW2OW3 By the same reason as the above we can easily get dA dA 2 dA 2 4 dw\dw dwxdw 3 dA dwidw dw ' 3 dw dw 4 3 4 z 4 4 4 dA dw dw dw 3 4 2 4 3 dA dwidw dw ' 3 dA 2 4 dw2dw dw2dw 4 2 dA 2 4 4 4 2 3 4 .Finally, dA , 4 4 Q Q Q Q C'Zi;iC'l«2C'W3C'^'4 — ~~ ^1234- ; 3.4 Computer implementation of the Hiisler-Reiss model The parameter A,j's are estimated with a quasi-Newton routine. The multivariate normal survival function is computed with the routine of Schervish (1984). Other routines that are needed are a numerical integration routines and routines for the functions: C i 3 , C i 2 , , Ci23,-Di,234, ^12,34,-Di23,4,-Di234i 2 3 Then C 1 < 2 4 , C 1 3 i 4 etc., can be com- puted from these routines by changing the parameters of the routines. 3.5 Data Analysis Joe (1990b) constructed some multivariate extreme value environmental data sets from Bay Area and from the Great Vancouver Regional District. The first data set consists of weekly maxima of ozone concentrations ( in parts per hundred million) for several monitoring stations in the Bay Area; the weeks were for the months of April to October in 1983-1987. For comparison with Joe (1990b), the same subset of 5 stations are used. The second data set consists of weekly maxima of SO2, NO2 and ozone concentrations (in parts per hundred milion ) for several monitoring stations in the Greater Vancouver Regional District: the weeks were for the months of April to September to October for July 1984 to October 1987. Joe (1990b) eliminated all weeks with too many 36 missing values, and for comparison the same reduced data are used. In Chapter 4, weeks with missing values are used as censored data. Details of exploratory data analysis of both data sets are given in Joe (1990b). Joe (1990b) compares the fit of the models (2-7-1) to (2-7-8) and their multivariate extensions, as well as other models in Joe (1990a). The models in Joe (1990a) which have some appealing theoretical properties were found to fit worse even though the models had more parameters. A partial explanation is that the parameters in these models are not easily interpretable as in (2-7-1) to (2-7-11). For comparison with a baseline, the likelihood of the multivariate normal copula with exponential margins was computed. The comparison of the Husler-Reiss model with the best fitting of (2-7-1) to (2-7-8) and with the multivariate normal copula is given in Tables 1 to 7. Only the parameter estimates for the Hiisler-Reiss model are given. Note that the likelihoods given are for the dependence parameters, after the univariate margins have been transformed to exponential survival functions (see Section 2.8). Also all likelihoods based on margins ( or subsets of data from subsets of stations) are given. This is to check on consistency of the parameter estimates for multivariate margins of different orders. This check provides an indication of the goodness of fit of the multivariate models. A final note for the tables is that for 3 or more stations, labelled as l , . . . , p , the parameter estimates are given in lexicographical order, that is, A12, A 1 3 , . . . , A i , A 2 3 , . . . , A p 2 p , . . . , A _i . p i P For example, consider the pairs of stations P T , SJ in Table 1, where the other stations are abbreviated C C , V A , S T . In the trivariate subset, ( C C , P T , SJ), the A parameter for (PT,SJ) is the third value (0.588); in the trivariate subset, (PT, SJ, V A ) , it is the first value (0.595). in the 4-variate subset, ( C C , P T , SJ, V A ) , it is the fourth value (0.567), etc. Recall from Chapter 2, that a smaller value of A means more dependence. Some overall conclusions from the tables are: 37 (1) The parameter estimates for a particular A are similar over different likelihoods. (2) Of the 37 bivariate likelihoods, the likelihood of the Hiisler-Reiss model is higher than those of both models (2-7-1) and (2-7-2) in 13 cases, and it is higher than that of the multivariate copula in 18 cases. (3) Of the 25 trivariate likelihoods, the likelihood of the Hiisler-Reiss model is higher than those of both models (2-7-3) and (2-7-4) in 11 cases. With the adjustment for the number of parameters ((2-7-3) and (2-7-4) have 2 parameters, the Hiisler-Reiss trivariate model has 3 parameters), through a penalty of half the number of parameters for the log likelihood( the Akaike information criterion), the Hiisler-Reiss model is only better for the first data set in Table 1. (4) Of the 9 four-variable likelihoods, the likelihood of the Hiisler-Reiss model is higher than those of the models (2-7-5) to (2-7-8) in 6 cases. This is true even with adjustment for the number of parameters ( 6 for the Hiisler-Reiss model and 3 for models (2-7-5) to (2-7-8) ). Overall, it appears that the Hiisler-Reiss model provides a better fit for Tables 1 and 3, and maybe Table 5. Therefore it is a useful model in addition to those in Joe (1990b). A n indication of whether the Hiisler-Reiss model is better can be seen from the estimates for the bivariate likelihoods. If the pattern of the A's does not fit the bivariate patterns that are possible with models (2-7-3) to (2-7-8), then the Hiisler-Reiss model should be better. ( The restrictions of the bivariate patterns for models (2-7-3) to (2-7-8) can be seen from taking the bivariate margins in these models.) 38 Table 1. Bay Area - Ozone(n=120) pay A Notation: Concord(CC),Pittsburg(PT),Vallejo(VA), San Jose(SJ) Santa RosT). Tikelihood for dependence parameters negative log likeli.hood estimate r^s subset T F T T normal H-R model CC, PT — O S 3 — 154.75 154.28 154.08 153.50 175.96 179.40 176.91 0.451 176.17 cd, SJ 190.30 0.532 190.70 190.38 192.41 CC, VA 210.47 211.60 211.03 214.24 0.708 CC, ST 199.06 199.85 199.34 200.48 0.595 PT, SJ 197.15 196.07 196.05 198.52 0.581 PT, VA 213.92 214.80 214.37 214.54 0.749 PT, ST 197.15 197.68 197.40 202.49 0.579 SJ, VA 213.52 214.56 213.96 218.74 0.740 SJ, ST 199.54 199.86 199.62 203.12 0.599 VA. ST 212.23 221.32 216.41 212.94 0.354,0.453 CCJ^SJ 0.5^88 222.92 CC,PT,VA 0.353,0.532 225.26 221.36 223.87 0.&0 243.86 CC,PT,ST 0.352,0.704 244.94 245.20 245.70 0.750 242.80 CC,SJ,VA 0.451,0.536 245.38 243.37 248.67 0.576 0.450,0.716 268.11 265.32 272.20 265.69 CC,SJ,ST 0.741 268.28 CC,VA,ST 0.533.0.708 271.62 271.40 273.54 0.606 0.595,0.583 267.09 265.27 269.90 264.27 PT,SJ,VA 0.579 287.04 PT,SJ,ST 0.594,0.755 290.80 288.08 291.56 0.?38 275.45 PT,VA,ST 0.582.0.752 279.04 277.97 278.94 0.d06 0.580.0.747 282.20 278.01 284.13 275.85 SJ,VA,ST 0.607 CC,PT,SJ 0.357,0.441 288.42 280.47 280.01 275.72 VA 0.539,0.567 0.581,0.574 CC,PT,SJ 0.355,0.441 311.54 305.04 303.61 299.69 0.70. 0.567 ST 0.738,0.741 CC,PT,VA 0.353,0.534 307.17 303.67 304.20 300.66 0.706,0.581 ST 0.754,0.607 CC,SJ,VA 0.450,0.538 329.94 324.70 329.57 321.01 0.714,0.577 ST 0.750,0.609 0.594,0.585 352.88 347.34 349.94 PT,SJ,VA 342.36 0.755,0.581 ST 0.746,0.610 39 Table 2. GVRD-Sulphur Dioxide(n=70) Notation: Confederation Park(CP), Second Narrows(SN), Anmore(AN),Rocky Point Park(RP). Likelihood for dependence parameters. estimate negative log likelihood subset H-R 6's T'S normal H-R model CP, SN 0.845 129.55 129.54 129.37 129.58 CP, A N 1.086 136.11 135.96 134.82 135.94 CP, R P 0.905 131.45 131.67 129.61 132.00 SN, A N 1.068 135.36 135.54 133.97 135.62 SN, R P 1.229 137.04 137.40 136.48 137.70 AN, RP 0.756 123.39 124.54 123.63 126.08 CP,SN,AN 0.840,1.068 193.83 193.82 192.16 193.44 190.28 192.97 189.26 190.89 183.89 188.01 182.94 187.05 187.61 190.59 187.30 190.93 241.86 247.15 240.07 244.05 1.062 CP,SN,RP 0.846,0.891 1.157 CP,AN,RP 1.045,0.909 0.750 SN,AN,RP 1.064,1.172 0.750 CP,SN,AN 0.840,1.060 RP 0.900,1.032 1.155,0.748 40 Table 3. GVRD-Nitrogen Dioxide(n=84) Notation: Kensington Park(KP), Confederation Park(CP), Anmore(AN), Rocky Point Park(RP). Likelihood for dependence parameters. negative log likelihood estimate subset H-R 6's r's normal H-R model KP, CP 0.378 112.97 113.10 120.15 114.06 KP, A N 0.415 119.26 119.31 124.01 119.59 KP, RP 0.362 110.64 .110.64 119.83 111.53 CP, A N 0.466 127.84 127.30 128.81 126.55 CP, R P 0.346 108.53 108.27 107.87 107.71 AN, RP 0.356 110.48 110.22 112.64 109.91 KP,CP,AN 0.379,0.417 147.96 146.20 154.95 146.27 133.11 126.51 137.25 128.55 136.99 .131.88 143.35 132.88 139.29 137.49 135.81 132.40 161.42 151.95 160.77 149.77 0.463 KP,CP,RP 0.382,0.363 0.343 KP,AN,RP 0.416,0.363 0.353 CP,AN,RP 0.463,0.345 0.356 KP,CP,AN 0.382,0.418 RP 0.363,0.457 0.342,0.352 41 Table 4. GVRD-Ozone(n=84) Notation: Marpole(MA), Confederation Park(CP), Anmore(AN), Rocky Point Park(RP). Likelihood for dependence parameters. negative log likelihood estimate subset H-R 8's r's normal H-R model MA, CP 0.521 89.17 89.23 88.00 89.69 MA, A N 0.829 101.63 102.47 101.80 104.54 MA, RP 0.763 99.94 100.50 100.88 102.21 CP, A N . 0.747 98.78 99.60 95.96 101.93 CP, R P 0.680 97.23 97.72 93.19 98.95 AN, RP 0.284 61.55 61.73 72.02 64.66 MA,CP,AN 0.526,0.803 129.02 131.08 126.47 133.31 125.73 129.42 123.81 130.07 103.97 105.08 115.26 108.96 101.13 101.70 107.48 105.76 130.01 133.28 137.97 136.71 0.741 MA,CP,RP 0.526,0.738 0.676 MA,AN,RP 0.816,0.747 0.284 CP,AN,RP 0.734,0.670 0.285 MA,CP,AN 0.505,0.763 RP 0.705,0.700 0.642,0.291 42 Table 5. GVRD-Station 5 (Confederation Park). (n=64) Likelihood for dependence parameters. estimate subset negative log likelihood H-R <5's r's normal H - R model 1.359 127.33 126.92 125.49 126.77 so ,o 1.033 122.83 122.88 122.10 123.01 N0 ,0 0.932 120.31 120.56 118.39 120.76 1.312,1.023 181.73 180.37 176.24 179.48 S0 ,N0 2 2 2 2 3 3 S0 ,N0 ,0 2 2 3 0.930 Table 6. GVRD-Station 7 (Anmore). (n=76) Likelihood for dependence parameters. estimate subset negative log likelihood H-R <5's T'S normal H - R model 1.167 148.68 148.81 146.10 148.84 so ,o 1.134 148.21 148.22 145.17 148.25 N0 ,0 0.526 118.65 119.11 124.25 120.64 1.132,1.112 189.37 190.11 192.69 191.96 S0 ,N0 2 2 2 2 3 3 S0 ,N0 ,0 2 2 3 0.526 43 Table 7. GVRD-Station 9 (Rocky Point Park). (n=63) Likelihood for dependence parameters. estimate subset negative log likelihood H-R tf's r's normal H-R model 1.065 120.87 121.48 118.51 121.99 so ,o 1.063 120.38 121.10 118.57 122.28 N0 ,0 0.622 107.55 106.99 106.00 106.30 1.019,1.050 162.60 163.61 159.92 164.04 S0 ,N0 2 2 2 2 3 3 S0 ,N0 ,0 2 2 3 0.625 44 Chapter 4. Analysis With Missing Data In statistical analysis, it may happen by accident or some other reason that some of the observations are missing, cannot be collected, or in some other way are not obtainable. In such a situation it is obvious that the routine method is not appropriate. Now we would like to reanalyze the data from Greater Vancouver Regional Dis- trict with missing values among the hourly measurements. If there are missing hourly measurements, we treat the daily or weekly maxima as right-censored. For right censored data, let x - be either the observed value or the right censored t value, and let 1 if Xi is right-censored Si 0 if Xi is not right-censored If f(t; 0) is a parametric family of models for the iid data and if S(t; 9) is the survival function, then the likelihood function with right-censored data is L = n^J(x ,0) - S(t ,9) , 1 S i 6 i i i see, for example, L'awless(1980) for details. We will use the procedure in Section 2 proposed by Joe (1990b) for fitting of a parametric multivariate distribution to independent and identically distributed p-vectors X i , . . . ,X . n The first step is to fit generalized extreme value distribution to the uni- variate margin separatly by maximum likelihood. That is, the j t h univariate margin is: Fj{x^ Hj,(Tj) h = e x p { - ( l + 7 [(x i - fij)/<Tj])+ }, lhi where (z)+ = max(z,0). Assuming that the measurements are iid, the likelihood is 45 where ~F(x; 7, [i, a) = 1 - F ( x ; 7, (i, a) and f(x; 7, fi, a) = — been added for the jth , and subscript j has component. After fitting the univariate margins, make the transformation ^• = (i+7i[(^-^)/^]r 1/7> so that original random variables become exponential random variables. The y ^ s are not treated as left censored data. Now for multivariate parameters, we take the univariate margins to be exponential. {Y\\, Yi ),. For bivariate data, suppose we have iid random pairs observe (Yu, Y , Sn, 6 ),... 12 ,{Y , Y , S , 12 nX n2 Sij = { ^12),and S ),where nl n2 1 if Yij is left censored, 0 if Y{j is not censored, i = 1 , . . . , n and j = 1,2. The bivariate cdf is F ( y , y ,6) x F is the survival function and F\ and F • • i(Yni, 2 2 = F = 1 - F - F + F, where x 2 are univariate margins of F. 2 Therefore we can get the following contribution to the likelihood of the bivariate parameter 6: P(Y =y ,Y 2 P(Y <y ,Y 2 1 1 1 1 = y )_ dF dy dy = y) _ dF_ _ df_ _ d~F dy dy dy 2 2 _ P(Y =y ,Y <y ) 1 1 2 2 _ P(Y <y .Y <y ) = 1 1 2 2 2 1 if ( £ i , M = (0,0) 2 7 2 2 8F_ _ dyi 9F _ dyi 2 = (1,0) i2 dF, dyiif (^1,^2) = F( ,y ,0) if (6 S ) 2 y i if (Siu6 ) iu = i2 (0,1) (l,l) In the trivariate case, the generalization is straightforward and some details are given below. and observe Suppose we have iid random variables (Y ,Y ,Y ,8 ,S ,S ), n 1 2 1 3 n 1 2 1 3 8^ = < (Y\\,Yi ,Y\ ), 2 ... ,(Y i, Y , Y ), 3 n (Y ,Y ,Y ,8 ,8 ,8 ), n l n 2 n 3 n l 1 if Yij is left censored, 0 if Yij is uncensored, 46 n 2 n 3 where n2 n3 i — 1,..., n and j = 1,2,3. The trivariate cdf is F(yi, y , y , 5 ) = 1 - F i ( y i ) - F ( y i ) - F ( y i ) + F ( , 2 2 2 + ^ i ( y i , y , 0) + 3 3 ^23(2/2,2/3,0) 3 12 yi 2 1 3 i2 2 - F( , y , y ,9) where JF is survival function and F\, F , F , F , F , and F 2 y ,5) yi 13 2 3 2 are the univariate and bivariate margins of F. Let C i = (Sn,Si2,Si ), then we have 2 3 3 P(Y =y ,Y = y ,Y = y )=: P(Y <y Y = y ,Y 1 1 1 u 2 2 2 3 2 dF dyidy dy 3 3 2 F = y ) =dyddy3 2 3 1 P(Y 1 = ,Y Vl < y ,Y 2 2 2 = y ,Y = y ,Y <y ) 1 2 2 3 1 2 2 3 P(Yr < y Y u 2 3 u 2 2 P{Y < ,Y 1 Vl 2 3 = 3 < y ,Y 2 3 dF dyi d F dyidy dF dyi dy 2 2 d F dy dy 3 dF dy\dy 2 _ 2 dF* dys _ _ dF dy 3 , dF-\3 " " dy 1 , 1 3 dFi ? dy 2 dF, , dF-\ 9 dyi " " dyi t dF?3 _ 1 3 dF 3 2 dy 2 2 dF_ _ a£ dy 2 3 = (0,1,0) = (0,0,1) if C*i = (1,1,0) if C\2 3 = (1,0,1) 2 dF-] 3 _ 9 F dyi dyi 2 3 23 dy3 < y ) = F(yi, 2/2,y ,0) 3 C\2 if C i 2 dy + 2 if 3 2 n 2 3 P(Y =y Y <y ,Y <y ) 1 2 < y ) = dF dy dy3 u 3 = y ,Y 2 dF dyidy = y ) =dF P(Y <y ,Y <y ,Y 1 = 3 = (0,0,0) 23 dyidyz 2 3 2 3 if C*i = (1,0,0) 2 2 dF = y ) = dyidy3 3 dF d F,3 3 2 P(Y if C i 2 if C ' i = (0,1,1) 23 if C*i 23 = (1,1,1) Clearly, this can be generalized to higher dimensions. Computer implementation of likelihood For models (2-7-1) to (2-7-8), the different contributions to the likelihood can easily be obtained using a symbolic manipulation software, such as Maple (Char et al., 1988, 5th edition), where derivatives can be output in Fortran format. Data analysis Here we reanalysize the Data Set 2 of Chapter 3 with the missing hourly values and censored weekly maxima. The results of the multivariate extreme value analysis are summarized in Table 8 — 10. The values of the maximum likelihood estimates are given in the left hand side of each table, and the negative log likelihoods are in the right side of each table. 47 A larger value of n means that the maxima for the station tend to be larger, a larger value of a means more spread and a larger value of 7 means that the hourly concentration has a heavier tail. Compared to the results of previous analysis, where weeks with many missing hourly measurements were treated as missing at random, we see that all parameter estimates of fi,cr and 7 in three tables are mostly slightly larger (compare Tables 3 to 5 of Joe 1990b). This is not a surprising result. After fitting the univariate margins, we obtain the transformed data which is left censored. The next step is to analyze the transformed data. The models for the multivariate exponential distribution were described before. The parameter estimates for the censored data are quite close to the estimates in Table 3 to 5 of Joe (1990b), but the relative values of the estimates are different in some cases. For trivariate case, the likelihoods for the 6's models and the r's models diverge whereas they are close in the bivariate case. A n explanation of the multivariate parts of Tables 8 to 10 is explained next for one case. For example, for (CP, SN, AN) in Table 8, the estimates of the parameters in model (2-7-3) are (1.567, 1.249), and the estimates of the parameters in (2-7-4) are (0.851, 0.507). For model (2-7-3) (respectively (2-7-4)), 1.567 (0.851) measures the dependence of the pair of stations (CP, SN), and 1.249 (0.507) measures the dependence of the pairs of (CP, A N ) and (SN, A N ) . The model (2-7-3) is better in Table 10, and the model (2-7-4) is better in Table 8 and 9. Conclusions are the same as in Joe (1990b) which means that there is no substantial difference with the missing values deleted. 48 Table 8. GVRD-Sulphur Dioxide(n=104) Notation: Confederation Park(CP), Second Narrows(SN), Anmore(AN),Rocky Point Park(RP). Likelihood for dependence parameters. subset estimate estimate negative log likelihood S's T'S <5's r's CP, SN 1.560 0.843 185.581 185.664 CP, A N 1.220 0.468 205.635 205.345 CP, R P 1.462 0.729 200.587 200.556 SN, A N 1.292 0.535 191.869 192.223 SN, R P 1.447 0.720 190.815 190.669 AN, RP 1.647 0.905 195.814 197.452 Multivariate CP,SN,AN 1.567,1.249 0.851, 0.507 284.370 283.937 CP,SN,RP 1.550,1.456 0.823, 0.674 277.888 279.263 AN,RP,CP 1.650,1.348 0.875, 0.562 291.025 293.808 AN,RP,SN 1.675,1.395 0.868, 0.596 278.128 282.235 Estimate for univariate margin parameter CP 0.313 1.857 0.834 SN 0.181 2.083 0.885 AN 0.261 3.215 2.024 RP 0.080 2.058 1.094 49 Table 9. GVRD-Nitrogen Dioxide(n=84) Notation: Kensington Park(KP), Confederation Park(CP), Anmore(AN),Rocky Point Park(RP). Likelihood for dependence parameters. estimate estimate <5's T'S 8's r's KP, CP 2.816 2.118 161.095 161.041 KP, A N 2.673 1.965 161.941 161.944 KP, RP 2.947 2.246 154.375 154.356 CP, A N 2.168 1.482 176.638 176.104 CP, R P 3.256 2.557 144.079 143.930 AN, RP 2.968 2.272 152.358 152.171 subset negative log likelihood Multivariate KP,CP,AN 2.741, 2.259 2.061, 1.715 219.517 215.709 RP,KP,CP 2.945, 2.718 2.364, 2.069 194.823 189.065 RP,KP,AN 2.974, 2.597 2.226, 2.150 198.191 192.105 RP,CP,AN 2.961, 2.419 2.2.9, 1.934 204.392 199.517 Estimate for univariate margin parameter KP 0.082 5.701 2.186 CP 0.179 6.093 2.442 AN 0.086 3.750 1.760 RP 0.078 5.228 1.658 50 Table 10. GVRD-Ozone(n=84) Notation: Marpole(MA), Confederation Park(CP), Anmore(AN),Rocky Point Park(RP). Likelihood for dependence parameters. estimate estimate 6's r's 6'a r's MA, CP 2.077 1.360 143.604 143.340 MA, A N 1.726 1.012 159.720 160.279 MA, RP 1.691 0.967 159.913 160.312 CP, A N 2.174 1.467 147.824 147.753 CP, R P 2.162 1.444 144.703 144.804 AN, RP 3.616 2.917 122.432 122.622 subset negative log likelihood Multivariate MA,CP,AN 2.061, 1.835 1.322, 1.140 211.020 210.141 MA,CP,RP 1.948, 1.948 1.328, 1.073 207.476 210.365 AN,RP,MA 3.606, 1.759 2.868, 0.963 189.255 191.630 AN,RP,CP 3.586, 2.131 2.852, 1.439 176.798 177.085 Estimate for univariate margin parameter MA -0.049 4.182 0.816 CP -0.019 5.054 1.384 AN 0.047 5.250 1.595 RP 0.091 5.054 1.761 51 Chapter 5. Conclusion and Further Research This thesis has applied maximum likelihood estimation of several parametric families of multivariate extreme value distribution, including the case where some maxima are right censored. The Husler-Reiss family has more flexibility, with one parameter for each bivariate margin, than the parametric families used in Joe (1990b). However, its form is much more complicated and for diminsion p > 3, time-consuming integrations and computations of the multivariate normal cdf are needed. Also, unlike the other families, symbolic manipulation software cannot be used to obtain the density function from the survival function. The maximum likelihood estimation takes 10-20 hours on a Sun SPARCstation 330 for p =4, and is expected to take weeks for p =5 ( the programming for this case has not been completed). If both the multivariate dependence parameters and the univariate marginal parameters are simultaneously estimated, the C P U time is much more than the preceding (this also has not been done). Fortunately, it seems possible from the maximum likelihood estimates of the bivariate likelihoods to check whether the simpler models in Joe (1990b) fit the data adequately, if so, the bivariate parameter estimates must follow certain patterns. If the simpler models appear not to be adequate, then the best model, among known ones, would likely be the Hiisler-Reiss model. Further research is needed to derive families of multivariate extreme value distributions that have a dependence parameter for each bivariate margin and that have a simpler form than the Husler-Reiss model (simpler in the sense that maximum likelihood estimation becomes possible for p > 5). Some positive work in this direction has been initiated. 52 BIBLIOGRAPHY Char, Bruce W . , Geddes, K . O . , Gonnet, G . H . , Monagan, M . B . , and Watt, S.M. (1989). Maple. 5th ed., W A T C O M Publications Limited. Coles, S. G . and Tawn, J . A . (1991). Modelling extreme multivariate events. J . R. Statist. Soc. 52, 377-392. Galambos, J . (1987). The Asymptitic Theory of Extreme Order Statistics, second edition. Kreiger Publishing Co., Malabar, F L Hiisler, J . and Reiss, R . - D . (1989). Maxima of normal random vectors between independence and complete dependence. Statistics &; Probability Letters 7 283-286. Joe H . (1990a). Families of min-stable multivariate exponential and multivariate extrem value distributions. Statistics & Probability Letters 9 (1990) 75-81. Joe H . (1990b). On fitting multivariate extreme value distributions, with applications to environmental data. Technical report. Department of Statistics, University of British Columbia. Lawless J . F . (1982). Statistical Models and Methods for Lifetime data. JOHN W I L E Y & SONS, Inc., New York. Resnick, S.I. (1987). Extreme Values, Regular Variation and Point Processes. New York: Springer-verlag. 53 Schervish, J . (1984). Multivariate normal probabilities with error bound. Appl. Stat.33 81- 94. Tawn J . A . (1988). Bivariate extreme value theory, models and estimation, Biometrika 75, 397-415. Tawn J . A . (1990). Modelling multivariate extreme value distributions. Biometrika 77, 245-253. 54
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimation with multivariate extreme value distributions...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimation with multivariate extreme value distributions with applications to environmental data Chen, Youshi 1991
pdf
Page Metadata
Item Metadata
Title | Estimation with multivariate extreme value distributions with applications to environmental data |
Creator |
Chen, Youshi |
Publisher | University of British Columbia |
Date Issued | 1991 |
Description | Several parametric families of multivariate extreme value distributions (Hüsler and Reiss 1989, Tawn 1990, Joe 1990a, 1990b) have been proposed recently. Applications to multivariate extreme value data sets are needed to assess the adequacy of the known families in their fit to data. Different families are compared in their range of multivariate dependence and their ease of use for maximum likelihood estimation. Some useful conclusions have been made from experience with several environmental data sets. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0302381 |
URI | http://hdl.handle.net/2429/29795 |
Degree |
Master of Science - MSc |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1991_A6_7 C42_5.pdf [ 3.38MB ]
- Metadata
- JSON: 831-1.0302381.json
- JSON-LD: 831-1.0302381-ld.json
- RDF/XML (Pretty): 831-1.0302381-rdf.xml
- RDF/JSON: 831-1.0302381-rdf.json
- Turtle: 831-1.0302381-turtle.txt
- N-Triples: 831-1.0302381-rdf-ntriples.txt
- Original Record: 831-1.0302381-source.json
- Full Text
- 831-1.0302381-fulltext.txt
- Citation
- 831-1.0302381.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0302381/manifest