Extremes 5:1, 33±44, 2002 # 2002 Kluwer Academic Publishers. Manufactured in The Netherlands.

Models for Stationary Max-Stable Random Fields MARTIN SCHLATHER Soil Physics Group, University of Bayreuth, 95440 Bayreuth, Germany E-mail: [email protected] [Received May 9, 2001; Revised July 2, 2002; Accepted August 13, 2002] Abstract. Models for stationary max-stable random ®elds are revisited and illustrated by two-dimensional simulations. We introduce a new class of models, which are based on stationary Gaussian random ®elds, and whose realizations are not necessarily semi-continuous functions. The bivariate marginal distributions of these random ®elds can be calculated, and they form a new class of bivariate extreme value distributions. Key words. bivariate extreme value distribution, dependence function, Gaussian random ®eld, max-stable random ®eld, rainfall modeling, simulation of max-stable processes AMS 1991 Subject Classi®cation.

1.

PrimaryÐ60G70 SecondaryÐ60G15, 60G55

Introduction

This work was motivated by the study of annual maxima of daily spatial rainfall, for which stationary max-stable processes might be good approximations. Modeling rainfall is a very complex task and a vast amount of literature on this topic exists both in statistics and applied sciences, see, for example, Barndorff-Nielsen et al. (1998) and the references therein. Let us mention that there are two main types of precipitation: (i) convective precipitation with usually a local area of high intensity and minor or no rainfall elsewhere, and (ii) cyclonic precipitation with variable rainfall all over the region. Clearly, these two types of processes should be modeled differently, and they will also affect the behavior of the annual maxima. A max-stable process Z is the limit process of maxima of i.i.d. random ®elds Yi …x†, x [ Rd . Namely, for suitable an …x† > 0 and bn …x† [ R, maxniˆ 1 Yi …x† n?? an …x†

Z…x† ˆ lim

bn …x†

;

x [ Rd ;

provided the limits exist. We assume that Z is stationary and has non-trivial univariate marginal distributions. Interestingly, there are two canonical construction principles for stationary max-stable random ®elds that re¯ect the two types of precipitation. The max-stable random ®eld Z that corresponds to the ®rst type has been ®rst considered by Smith (1990), see also Coles and Tawn (1996). Here,

34

SCHLATHER

Z…x† ˆ sup sf …x

y†;

…y;s† [ P

x [ Rd ;

…1†

for a Poisson point process P on Rd 6…0; ?† with intensity measure L and R dL…y; s† ˆ dy s 2 ds. The function f on Rd is non-negative and f …x† dx ˆ 1, i.e. f is greater than e on a domain of ®nite size only, for any e > 0. Theorem 1 gives a generalization of (1). Theorem 2 picks up the idea for the second type of precipitation. There, the max-stable process is based on independent, stationary random ®elds that have the same correlation structure. Along with the two theorems, Section 2 provides several examples, including formulae for the bivariate distributions. In Section 3, a simulation algorithm is described.

2.

Stationary max-stable random ®elds

In this section, three classes of stationary, max-stable random ®elds Z are presented and illustrated by explicit numerical examples and two-dimensional simulations. Theorem 1 gives a generalization of Smith's approach (Smith, 1990), replacing the deterministic shape function f in (1) by a random shape Y, cf. Deheuvels (1983) and de Haan and Pickands (1986). Theorem 2 introduces a new class that is based on stationary random ®elds with ®nite expectation. For sake of completeness the section ends with a rather theoretical class based on periodic functions. Denote by Fa the FreÂchet distribution with parameter a > 0,  Fa …s† ˆ exp

1 sa

 ;

s > 0:

Due to Proposition 5.10 in Resnick (1987) we may assume without loss of generality that the univariate margins of the max-stable process Z are unit FreÂchet F1 . However, the simulation results have been ®nally transformed to F2 margins to the advantage of a convenient choice of the grey scales in Figures 1 to 4. Theorem 1 deals with arbitrary random functions Y whose positive part, maxf0; Y… ? †g, is integrable. More precisely, Y ˆ Y …o† is a random variable on a probability space d …O; a;R P† with values in a Polish subspace U of RR , Y … ? † … ? † is measurable on O6Rd , and E Rd maxf0; Y…x†g dx [ …0; ?†. Denote by o the origin of the coordinate system. R Theorem 1: Let Y be a measurable random function and m ˆ E Rd maxf0; Y…x†gdx [ …0; ?†. Let P be a Poisson process on Rd 6…0; ?† with intensity measure dL…y; s† ˆ m 1 dy s 2 ds, and Yy;s i.i.d. copies of Y; then Z…x† ˆ sup sYy;s …x …y;s† [ P

y† ˆ sup s maxf0; Yy;s …x …y;s† [ P

y†g;

is a stationary max-stable process with unit FreÂchet margins.

…2†

35

MODELS FOR STATIONARY MAX-STABLE RANDOM FIELDS

Proof: Due to the in®nite divisibility of the Poisson processes de Haan's (1984) de®nition for max-stable processes yields directly that Z is max-stable. The margins are unit FreÂchet since …y; s; Yy;s † equals in distribution a Poisson process on Rd 6…0; ?†6U with intensity measure L…B6I6V† ˆ m

1

Z

Y

jBj P …V†

I

2

s

ds;

for measurable subsets B  Rd , I  …0; ?†, and V  U. Then P…Z…o†  t†, t > 0, equals the void probability (Stoyan et al., 1995) of the set 

…x; s; x† [ Rd 6…0; ?†6U : sx… x†  t ;

i.e. log P…Z…o†  t† ˆ m

1

Z Z U

Rd

Z

?

t= maxf0;Y…

x†g

s

2

ds dx dPY ˆ t

1

;

see also Choquet's theorem in Matheron (1975).

…3† &

Theorem 1 is closely related to work by several other authors. Deheuvels (1983) deals with processes on the one-dimensional lattice; from Deheuvels' (1983) point of view Z is a moving maximum process. De Haan and Pickands (1986) characterize stationary maxstable processes on R1 that are continuous in probability. Finally Gine et al. (1990), see also Heinrich and Molchanov (1994), investigate not necessarily stationary max-stable random ®elds based on upper semi-continuous functions, hence on Matheron's Boolean functions, see Jeulin and Jeulin (1981), Matheron (1975), Norberg (1986). Statistical inference for and modeling of spatial extremes have been discussed by Coles (1993), Coles and Tawn (1996), Davis and Resnick (1993), for example. We give two examples of spatial models and their bivariate marginal distribution. The latter is completely characterized by the dependence function A, which is de®ned as (Pickands, 1981; Tawn, 1988) A…a† ˆ Ax …a† ˆ

log P…maxf…1

a†Z…0†; aZ…t†g  1† for

a [ ‰0; 1Š and

x [ Rd :

Smith (1990), see also Coles (1993), calls Z a Gaussian extreme value process if Y is deterministic and equals a multivariate Gaussian density with covariance matrix S. The bivariate marginal distribution equals (Smith, 1990) log P…Z…o†  s; Z…x†  t† ˆ s

1

  a 1 t F ‡ log ‡t 2 a s

1

  a 1 s F ‡ log ; 2 a t

…4†

36

SCHLATHER

where 2

>

a ˆx S

1

x

and

F…u† ˆ …2p†

1=2

Z

u ?

exp… v2 =2† dv:

Figure 1 shows a simulation for the two-dimensional isotropic Gaussian density with variance 9/8. Denote by b…x; r† the ball with radius r centered at x. If Y is the indicator function for a ball of random radius R, cf. Ancona-Navarrete and Tawn (2001), then Z has the bivariate distribution log P…Z…o†  s; Z…x†  t† ˆ

1 Ejb…o; R† n b…x; R†j ‡ ; minfs; tg maxfs; tgEjb…o; R†j

i.e. Ax …a† ˆ

2a…1 a† 2a…1 a† Ejb…o; R† n b…x; R†j ‡ ? : 1 j2a 1j 1 ‡ j2a 1j Ejb…o; R†j

If d ˆ 2 and P…R  r† ˆ 1

…1

4r2 †

log P…Z…o†  s; Z…x†  t† ˆ

1=2

for r [ ‰0; 1=2Š then

1 g…x† ‡ ; minfs; tg maxfs; tg

Figure 1. Smith's Gaussian extreme value process with variance 9/8; the random ®eld has been transformed to F2 -margins.

MODELS FOR STATIONARY MAX-STABLE RANDOM FIELDS

37

where 3 g…x† ˆ minf1; jxjg 2

3 1 minf1; jxjg : 2

…5†

Figure 2 presents a simulation of this process. Formula (5) can be obtained without 1=2 calculation since P…R  r† ˆ 1 …1 4r 2 † , r [ ‰0; 1=2Š, is the radius distribution of the circles obtained by a two-dimensional cross section through a stationary and isotropic set of non-overlapping spheres of radius 1=2, hence g is the spherical variogram (Wackernagel, 1998) that has form (5). In the above construction two facts have been important: (i) the distribution function of P is invariant with respect to vectors …x; 0† for x [ Rd ; (ii) the expectation of the integral of Y over Rd is ®nite. In some sense the two properties of P and Y can be exchanged, and we obtain a new class of models. Theorem 2: Let Y be a stationary process on Rd with m ˆ E maxf0; Y…o†g [ …0; ?†, and let P be a Poisson process on …0; ?† with intensity measure dL…s† ˆ m 1 s 2 ds. Then Z…x† ˆ max sYs …x† ˆ max s maxf0; Ys …x†g; s[P

s[P

…6†

is a stationary max-stable process with unit FreÂchet margins, where the Ys are i.i.d. copies of Y for all s in …0; ?†. The proof is similar to that of Theorem 1.

Figure 2. Max-stable random ®eld based on random disks; the random ®eld has been transformed to F2 margins.

38

SCHLATHER

The model given by (6) allows for practical interpretation. The single events sYs , e.g. the daily spatial rainfall, have all the same dependence structure and differ only in the magnitude s. Further, any single event, conditioned on the magnitude s, is accessible to standard geostatistical analysis, provided the variance of Y is ®nite. By way of contrast, Z does not have ®nite expectation, thus cannot be analyzed by such methods. We specify to stationary Gaussian processes. Let Y be a stationary standard normal random ®eld with correlationp function r…t†. Let P be a Poisson process on …0; ?† with  intensity measure dL…s† ˆ 2ps 2 ds. Then we call the process Z given by (6) an extremal Gaussian process. Two-dimensional simulations are shown in Figures 3 and 4, which are based on the same random seed. Extremal Gaussian processes do not belong to the class given in Gine et al. (1990) for two reasons: (i) for some covariance functions the realizations of the corresponding Gaussian random ®elds cannot be made upper semicontinuous; (ii) the L? -norm of a realization of a Gaussian random ®eld equals ?. The extremal Gaussian process differs from Smith's Gaussian extreme value process (Smith, 1990), and the processes of Casson and Coles (1999, 2000) who modeled the parameters of the generalized extreme value distribution essentially by a Gaussian random ®eld. A new class of bivariate marginal distributions is de®ned by the extremal Gaussian process, which is given by s!   1 1 1 st ‡ log P…Z…o†  s; Z…x†  t† ˆ 1 ‡ 1 2…r…x† ‡ 1† : 2 t s …s ‡ t†2

…7†

Figure 3. Extremal Gaussian process with exponential correlation model, r…h† ˆ exp… h†; the random ®eld has been transformed to F2 -margins.

MODELS FOR STATIONARY MAX-STABLE RANDOM FIELDS

39

Figure 4. Extremal Gaussian process with Gaussian correlation model, r…h† ˆ exp… h2 †; the random ®eld has been transformed to F2 -margins.

See the appendix for the proof. The correlation coef®cient r may take any value in ‰ 1; 1Š: if r 6ˆ 1 the two FreÂchet variables, Z…o† and Z…x†, are positively correlated and for r ˆ 1 we get complete dependence. The positive correlation of the tails for r 6ˆ 1 has two sources. Since the intensity measure of the Poisson process P is heavy tailed, the tail behavior of the extremal Gaussian process is essentially solely determined by the point of P with the largest value. On the other hand, if r 6ˆ 1 the probability of two jointly Gaussian random variables being positive is greater than zero. Hence, both random variables take, with positive probability, simultaneously large values, i.e. the FreÂchet variables are positively correlated. If r ˆ 1, the symmetry of the bivariate Gaussian distribution and the in®nitely divisibility of the Poisson process imply that the two corresponding FreÂchet variables can be simulated by two independent Gaussian random variables and two independent copies of the Poisson process, i.e. the FreÂchet variables are independent, as given by formula (7). Random indicator functions 1B of compact sets B and stationary standard Gaussian random ®elds Y might be combined to Z…x† ˆ max s1By;s …x …y;s† [ P

y†Yy;s …x†;

under the technical conditions of Theorem 1. Here, P is the Poisson process on p  Rd 6…0; ?† with intensity measure dL…y; s† ˆ 2p…EjBj† 1 dy s 2 ds and By;s and Yy;s are i.i.d. copies of B and Y, respectively. For the bivariate distribution of Z we get

40

SCHLATHER

log P…Z…o†  s; Z…x†  t†   1 1 EjB \ …x ‡ B†j ‡ 1 1 ˆ s t 2EjBj

s!! st 1 2…r…x† ‡ 1† : …s ‡ t†2

Finally, periodic functions, which seem to be of theoretical interest only, are mentioned for completeness. These functions are neither integrable nor ergodic. A simple example is givenR for illustration. Let Y be a deterministic, periodic function on R with period 1, 1 m ˆ 0 maxf0; Y…x†g dx > 0, and P the Poisson process on ‰0; 1†6…0; ?† with intensity measure dL…y; s† ˆ m 1 1‰0;1† …y† dy s 2 ds. Then Z…x† ˆ max sY…x

y† ˆ max s maxf0; Y…x

…y;s† [ P

…y;s† [ P

y†g;

is a max-stable stationary process on the real axis with unit FreÂchet margins.

3.

Simulation

In both, Theorems 1 and 2, the construction of a max-stable random ®eld Z involves the maximum over an in®nitely number of copies of the random function Y. In simulation studies, the number of realizations for Y are necessarily ®nite, but under certain conditions, see Theorem 4, we may get nonetheless exact simulations for Z on a ®nite sampling window B. If the conditions are not met, it is plausible that the approximations are still good. To this end, the simulation starts with the largest value of the magni®cation s, and determines on the basis of a stopping rule for s small enough. Lemma 3: Let c be a positive constant and B a subset of Rd with ®nite Lebesgue measure jBj. Let Ui be i.i.d. uniformly distributed on B and xi i.i.d. standard exponential variables, i.e. P…xi  s† ˆ 1 exp… s† for s  0. Then the random set (

, Un ; cjBj

n X kˆ1

! xk

) : n ˆ 1; 2; . . .

is a simple Poisson process on B6…0; ?† with intensity measure dL…y; s† ˆ c dy s

2

ds.

 Pn Proof: It is well known that P ˆ k ˆ 1 xk : n ˆ 1; 2; . . . is a Poisson process on the positive real axis with intensity 1. Application of the transformation t ° cjBjt 1 to the points of P yields a Poisson process whose density measure L satis®es L…‰a; ?†† ˆ cjBja 1 . As the Ui are i.i.d. uniformly distributed on B, we get the assertion of the theorem. & Denote by Br the in¯ation of B by the ball b…o; r†, i.e. Br ˆ |x [ B b…x; r†.

MODELS FOR STATIONARY MAX-STABLE RANDOM FIELDS

41

Theorem 4: Let Y, P, and Z be given by Theorem 1. Assume that Y is uniformly bounded by C [ …0; ?† and has support in b…o; r† for some r [ …0; ?†. Let B be a compact set, Yi be i.i.d. copies of Y, Ui be i.i.d. uniformly distributed on Br , xi be i.i.d. standard exponential variables, and P, Yi , Ui , and xi be all independent. Then, on B, Z * …x† ˆ

jBr j Yn …x Un † sup P ; n m n k ˆ 1 xk

x [ B;

equals the max-stable random ®eld Z in distribution, and Z * …x† ˆ

 jBr j Yn …x Un † sup P : n ˆ 1; . . . ; m; n m k ˆ 1 xk  C Yn …x Un † and m is such that Pm  max P n 1nm k ˆ 1 xk k ˆ 1 xk

…8†

almost surely: Proof: Lemma 3 yields that the process Z* P equals Z on B in distribution. Equality (8) is an immediate consequence of Yn  C and nk ˆ 1 xk being a non-decreasing sequence in n. & In case Y is a stationary, uniformly bounded random ®eld, an analogous result to Theorem 4 holds. For other random functions whose support is not included in a ball b…o; r† or which are not uniformly bounded by a constant C, approximations for r and C should be used. For Smith's Gaussian extreme value process we choose r such that the j…r† ˆ 0:001 for the standard normal density j; for the extremal Gaussian process we choose C ˆ 3. With this speci®cation, the algorithm needs about 3 and 0.2 seconds on a 1.8 GHz PC to generate Figures 1 and 2, respectively, which consist of 5126512 points; the generation of Figures 3 and 4 needs about 3 minutes each. The simulation algorithm is available as part of a contributed package on random ®eld simulation (Schlather, 2001) to R (Ihaka and Gentleman, 1996), see http://cran.r-project.org/.

Appendix: Bivariate marginal distribution of the extremal Gaussian process The extremal Gaussian process is de®ned as Z…x† ˆ max sYs …x†; s[P

where Ys are i.i.d. stationary Gaussian processes p with correlation function r and P is a Poisson process on …0; ?† with intensity 2pr 2 dr. We have

42

SCHLATHER

p Z Z log P…Z…o†  s; Z…x†  t† ˆ 2p

?

minfs= maxf0;Y…o†g;t= maxf0;Y…x†gg

z

2

dz dPY…o†;Y…x†

  p Z Y…o† Y…x† ˆ 2p max 0; ; dPY…o†;Y…x† : s t

We write r instead of r…x† for short, and de®ne r2 †;

Z ˆ 2…1 and Cˆ

p 2p…1

r2 †1=2 ˆ

p pZ:

Then,   p Z Y…o† Y…x† ; dPY…o†;Y…x† 2p max 0; s t Z ? Z ? ˆC 1 maxf0; v=s; u=tge Z ??Z ?? 2 ˆC 1 maxfv=s; u=tge …v 0

0

1

‡ …sC†

Z

Z

0 ?

?

0

ve

…v2

…v2

2rvu ‡ u2 †=Z

2rvu ‡ u2 †=Z

2rvu ‡ u2 †=Z

dv du

dv du

dv du ‡ …tC†

1

Z

?

Z

0

0

?

…v ‡ u†e

…v2

ue

…v2

2rvu ‡ u2 †=Z

dv du:

De®ning Z A…r† ˆ

?

Z

0

?

0

…v2

ve

2rvu ‡ u2 †=Z

dv du ˆ

1 2

Z 0

?

Z

?

0

and Z B…x† ˆ

0

?

Z 0

?

1v>ux v e

…v2

2rvu ‡ u2 †=Z

dv du;

x  0;

2rvu ‡ u2 †=Z

dv du; …9†

43

MODELS FOR STATIONARY MAX-STABLE RANDOM FIELDS

we get   p Z Y…o† Y…x† 2p max 0; ; dPY…o†;Y…x† s t     1 B…st 1 † B…ts 1 † 1 1 1 ‡ ‡ A… r†: ˆ ‡ C s t C t s

…10†

Let Z a…x† ˆ

Z

? 0

?

ux

…v

ru†e

…v2

2rvu ‡ u2 †=Z

ZC dv du ˆ q ; 4 x2 2rx ‡ 1

then Z B…x† ˆ

?

Z

?

…v ru†e 0 ux Z ?Z ? 2 ‡r ue …v 0

ˆ a…x†

0

rB…x

1

…v2

2rvu ‡ u2 †=Z

2rvu ‡ u2 †=Z

Z dv du

r

?

0

Z u

ux

0

e

…v2

2rvu ‡ u2 †=Z

dv du

dv du

† ‡ rA…r†:

Hence, B…x† ˆ

a…x†

ra…x

1

† ‡ r…1 1 r2

r†A…r†

:

Since (9) can be calculated by the coordinate transformation …w Z A…r† ˆ

? ?

Z

?

jzj

2we

2…w2 …1

r† ‡ z2 …1 ‡ r††=Z

dw dz ˆ

we get by (10) that   p Z Y…o† Y…x† 2p max 0; ; dPY…o†;Y…x† s t 0 !1=2 1   1 1 1 @ st A: ‡ ˆ 1 ‡ 1 2…r ‡ 1† 2 2 t s …s ‡ t†

z; w ‡ z† ° …v; u† as

ZC ; 4…1 r†

44

SCHLATHER

Acknowledgment The author thanks Jonathan Tawn for many hints and discussions. The research work has been supported by the EU TMR network ERB-FMRX-CT96-0095, and the German Federal Ministry of Research and Technology grant PT BEO 51-0339476C.

References Ancona-Navarrete, M.A. and Tawn, J.A., ``Diagnostics for extremal dependence in spatial processes,'' Submitted to Extremes, 2002. Barndorff-Nielsen, O.E., Gupta, V.K., PeÂrez-Abreu, V., and Waymire, E., (eds), Stochastic Methods in Hydrology, World Scienti®c, Singapore, New Jersey, London, Hong Kong, 1998. Casson, E. and Coles, S.G., ``Spatial regression models for extremes,'' Extremes 1, 449±468, (1999). Casson, E. and Coles, S.G., ``Simulation and extremal analysis of hurricane events,'' J. R. Statist. Soc., Ser. C 49, 227±245, (2000). Coles, S.G., ``Regional modelling of extreme storms via max-stable processes,'' J. R. Statist. Soc., Ser. B. 55, 797±816, (1993). Coles, S.G. and Tawn, J.A., ``Modelling extremes of the areal rainfall process,'' J. R. Statist. Soc., Ser. B 58, 329± 347, (1996). Davis, R.A. and Resnick, S.I., ``Prediction of stationary max-stable processes,'' Ann. Appl. Probab. 3, 497±525, (1993). de Haan, L., ``A spectral representation for max-stable processes,'' Ann. Probab. 12, 1194±1204, (1984). de Haan, L. and Pickands, J., ``Stationary min-stable stochastic processes,'' Probab. Th. Rel. Fields. 72, 477±492, (1986). Deheuvels, P., ``Point processes and multivariate extreme values,'' J. Multivar. Anal. 13, 257±272, (1983). GineÂ, E., Hahn, M.G., and Vatan, P., ``Max-in®nitely divisible and max-stable sample continuous processes,'' Probab. Th. Rel. Fields 87, 139±165, (1990). Heinrich, L. and Molchanov, I.S., ``Some limit theorems for extremal and union shot-noise processes,'' Math. Nachr. 168, 139±159, (1994). Ihaka, R. and Gentleman, R., ``R: A language for data analysis and graphics,'' J. Comput. Graph. Stat. 5(3), 299± 314, (1996). Jeulin, D. and Jeulin, P., ``Synthesis of rough surfaces of random morphological functions,'' Stereo. Iugosl. 3, 239±246, (1981). Matheron, G., Random Sets and Integral Geometry, John Wiley & Sons, New York, 1975. Norberg, T., On the existence and convergence of probabiity measures on continuous semi-lattices. Technical Report 148, Center for Stochastic Processes, University of North Carolina, 1986. Pickands, J., ``Multivariate extreme value distributions,'' Bull. Int. Stat. Inst. 49, 859±878, (1981). Resnick, S.I., Extreme Values, Regular Variation, and Point Processes, Volume 4, Springer, New York, Berlin, 1987. Schlather, M., ``Simulation and analysis of random ®elds,'' R News 1(2), 18±20, (2001). Smith, R.L., ``Max-stable processes and spatial extremes,'' Unpublished manuscript, 1990. Stoyan, D., Kendall, W.S., and Mecke, J., Stochastic Geometry and its Applications (second edition), John Wiley & Sons, Chichester, 1995. Tawn, J.A., ``An extreme-value theory model for dependent observations,'' J. Hydrology 101, 227±250, (1988). Wackernagel, H., Multivariate Geostatistics (second edition), Springer, Berlin, Heidelberg, 1998.