Report

Asymptotics for the Spatial Extremogram Richard A. Davis* Columbia University Collaborators: Yongbum Cho, Columbia University Souvik Ghosh, LinkedIn Claudia Klüppelberg, Technical University of Munich Christina Steinkohl, Technical University of Munich Zagreb June 5, 2014 1 Plan Warm-up • Desiderata for spatial extreme models The Brown-Resnick Model • Limits of Local Gaussian processes • Extremogram • Estimation—pairwise likelihood • Simulation examples Estimating the extremogram • Asymptotics in random location case Composite likelihood Simulation examples Two examples Zagreb June 5, 2014 2 Extremogram in Space Setup: Let be a stationary (isotropic?) spatial process defined on ∈ ℝ2 (or on a regular lattice ∈ ℤ2 ). 5 4 3 2 1 14 12 10 10 8 8 6 6 4 4 2 Zagreb June 5, 2014 12 14 2 3 Extremal Dependence in Space and Time Zagreb June 5, 2014 4 Illustration with French Precipitation Data Data from Naveau et al. (2009). Precipitation in Bourgogne of France; 51 year maxima of daily precipitation. Data has been adjusted for seasonality and orographic effects. Zagreb June 5, 2014 5 Warm-up Desiderata for Spatial Models in Extremes: {Z(s), s D} 1. Model is specified on a continuous domain D (not on a lattice). 2. Process is max-stable: assumption may be physically meaningful, but is it necessary? (i.e., max daily rainfall at various sites.) 3. Model parameters are interpretable: desirable for any model! 4. Estimation procedures are available: don’t have to be most efficient. 5. Finite dimensional distributions (beyond second order) are computable: this allows for computation of various measures of extremal dependence as well as predictive distributions at unobserved locations. 6. Model is flexible and capable of modeling a wide range of extremal dependencies. Zagreb June 5, 2014 6 Building a Max-Stable Model in Space Building blocks: Let () be a stationary Gaussian process on ℝ2 with mean 0 and variance 1. • Transform the () processes via () = −1/log(Φ(()), where F is the standard normal cdf. Then () has unit Frechet marginals, i.e., ( () ≤ ) = exp{−1/}. Note: For any (nondegenerate) Gaussian process (), we have lim ∞ > 0 > ) = 0. and hence lim ∞ > 0 > ) = 0. In other words, • observations at distinct locations are asymptotically independent. • not good news for modeling spatial extremes! Zagreb June 5, 2014 7 Building a Max-Stable Model in Space-Time Now assume () is isotropic with covariance function Cov ℎ , 0 = ℎ = exp − ℎ , where > 0 is the range parameter and α (0,2] is the shape parameter. Note that 1 − ℎ ~ ℎ =: (ℎ) ℎ → 0, and hence log ( 1 − ℎ) → ℎ , where = log 1 − . It follows that Cov ℎ , 0 Zagreb June 5, 2014 = ℎ ~1 − (ℎ)/ log . 8 Then, Building a Max-Stable Model in Space-Time 1 ≔ − =1 1 log Φ → () on ℝ2 . Here the are IID replicates of the GP , and is a BrownResnick max-stable process. Representation for (): ∞ = exp{ − ()} =1 where • are points of a Poisson process with intensity measure −2 . • ( ) are iid Gaussian processes with mean zero and variogram , . . , , Zagreb June 5, 2014 = + − − . 9 Then, Building a Max-Stable Model in Space-Time 1 ≔ − =1 1 log Φ → () on ℝ2 . Bivariate distribution function: ℎ ≤ , 0 ≤ = exp{− , ; } where , ; = −1 Φ log / 2√ + √ + −1 Φ log / 2√ + √ , and = ℎ . Note • V(x,y; ) → x-1 + y-1 as → ∞ (asymptotic indep) • V(x,y; ) → 1/min(x,y) as → 0 (asymptotic dep) Zagreb June 5, 2014 10 Scaled and transformed Gaussian random fields (fixed time point) −1 log(Φ( ) −1 log(Φ( ) Zagreb June 5, 2014 11 Scaled and transformed Gaussian random fields (fixed time point) 1 ≔ Zagreb June 5, 2014 − =1 1 log Φ 12 Lattice vs cont space Setup: Let be a RV stationary (isotropic?) spatial process defined on ∈ ℝ2 (or on a regular lattice ∈ ℤ2 ). Consider the former—latter is more straightforward. ℎ ∈ ℝ2 , ℎ = lim + ℎ ), regular grid 10 8 6 y 4 2 0 0 2 4 y 6 8 10 random pattern 0 2 4 6 x Zagreb June 5, 2014 8 10 0 2 4 6 8 10 x 14 France Precipitation Data regular or random? K-Function 0.5 1.0 distance Zagreb June 5, 2014 19000 18500 0.5 y L(t) 1.0 19500 1.5 20000 Site Locations 1.5 5500 6000 6500 7000 7500 8000 x 15 Regular grid 0 2 4 y 6 8 10 regular grid 0 2 4 h = 1; # of pairs = 4 6 x h = sqrt(2); # of pairs = 4 8 10 h = sqrt(10); # of pairs = 8 h = sqrt(18); # of pairs = 4 h = 2; # of pairs = 4 Zagreb June 5, 2014 16 Random pattern 0 2 4 y 6 8 10 random pattern 0 2 4 h = 1; # of pairs = 0 6 8 10 x h = 1 ± .25 Zagreb June 5, 2014 17 Random pattern 8 9 Zoom-in 4 5 6 y 7 buffer 0 1 2 3 4 5 x Estimate of extremogram at lag h = 1 for red point: weight “indicators of points” in the buffer. Bandwidth: half the width of the buffer. Zagreb June 5, 2014 18 Random pattern 0 2 4 y 6 8 10 random pattern 0 2 4 6 8 10 x Note: • Expanding domain asymptotics: domain is getting bigger. • Not infill asymptotics: insert more points in fixed domain. Zagreb June 5, 2014 19 Estimating the extremogram--random pattern Setup: Suppose we have observations, 1 , … , at locations 1 , … , of some Poisson process with rate in a domain ↑ ℝ2 . Here, = = number of Poisson points in , ~ . Weight function (): Let ⋅ be a bounded pdf and set 1 = 2 , where the bandwidth → 0 and 2 → ∞. Zagreb June 5, 2014 20 Estimating extremogram--random pattern , ℎ = lim ( + ℎ , )/ ), ℎ ∈ ℝ2 Kernel estimate of : , ℎ = 2 | | ≠=1 ℎ − + ∈ ( ∈ ) | | =1 ∈ , ℎ = 2 (ℎ + 1 − 2 ) 1 ∈ ( 2 ∈ ) 2 (1 , 2 ) | | ∈ () Note: 2 (1 , 2 )= 1 2 1 ≠ 2 is product measure off the diagonal. Zagreb June 5, 2014 21 Limit Theory Theorem: Under suitable conditions on , (i.e., regularly varying, mixing, local uniform negligibility, etc.), then 1 2 2 , ℎ − ,, ℎ → 0, Σ , where ,, ℎ is the pre-asymptotic extremogram, ,, ℎ = ( + ℎ , )/ ), ℎ ∈ ℝ2 , ( is the 1 − 1/ quantile of ). Remark: The formulation of this estimate and its proof follow the ideas of Karr (1986) and Li, Genton, and Sherman (2008). Zagreb June 5, 2014 22 Limit theory Asymptotic “unbiasedness”: , ℎ is a ratio of two terms; , ,, (ℎ) ℎ = , will show that both are asymptotically unbiased. Denominator: By RV, stationarity, and independence of and , , = | | ∈ () = 0 ∈ = 0 ∈ → () Zagreb June 5, 2014 23 Limit Theory 2 Numerator: = = 2 1 (ℎ + 1 − 2 )( 1 ∈ ℎ + 1 − 2 ( 0 ∈ , 2 − 1 ∈ ) 2 1 2 1 2 ℎ+1 −2 (2 − 1 ) 1 2 where ℎ = 0 ∈ , ℎ ∈ . Making the change of variables = ℎ+1−2 and = 2 , the expected value is 1 − +ℎ = Zagreb June 5, 2014 ∩( − +ℎ) − +ℎ (ℎ − )) ℎ − | ∩ − + ℎ |/|Sn | 24 Limit Theory − +ℎ ℎ − → | ∩ − + ℎ | Sn , ℎ = , ℎ . ℝ2 Remark: We used the following in this proof. • | ∩ − +ℎ Sn | → 1 and − +ℎ → ℝ2 . • ℎ − = 0 ∈ , ℎ − ∈ → , ℎ . Need a condition for the latter. Zagreb June 5, 2014 25 Limit theory Local uniform negligibility condition (LUNC): For any , > 0, there exists a ′ such that | − 0 | lim sup sup > < . ′ n < Proposition: If is a strictly stationary regularly varying random field satisfying LUNC, then for → 0, (0) + ∈ , ∈ → , This result generalizes to space points, 0, 1 + , … , + . Zagreb June 5, 2014 26 Limit Theory Outline of argument: • Under LUNC already shown asymptotic unbiasedness of numerator and denominator. • , → • ,, ℎ → , ℎ with , ℎ = , / ℎ . Strategy: Show joint asymptotic normality of , and ,, ℎ var , → Zagreb June 5, 2014 () + ℝ2 , ⟹ , ℎ → (ℎ) 27 Limit Theory Step 1: Compute asymptotic variances and covariances. i. ii. var , → 2 () + ℝ2 , var , (ℎ) → 1 2 (ℎ) ℝ2 2 Proof of (i): Sum of variances + sum of covariances 2 , = 2 | | + 2 | | () → + Zagreb June 5, 2014 ℝ2 1 ∈ (1 ) ( 1 ∈ , 2 ∈ ) 2 (1 , 2 ) , 28 Limit Theory Step 2: Show joint CLT for , and ,, ℎ using a blocking argument. Idea: Focus on ,, ℎ . Set = 1 2 2 Zagreb June 5, 2014 1 2 (ℎ + 1 − 2 )( 1 ∈ 29 Limit Theory Subdivide = 0, =∪=1 2 into big blocks and small blocks. where has diminesions × and size = 2 0 2 4 y 6 8 10 Zagreb June 5, 2014 0 2 4 6 8 10 30 Limit Theory Insert block = inside with dimension − blocks. × ( − ): Subdivide 0, 2 into big blocks and small 0 2 4 y 6 8 10 − 2 × and size = 2 || = =∪ where has diminesions =1 Zagreb June 5, 2014 0 2 4 6 8 10 31 Limit Theory Recall that is a (mean-corrected) double integral over × , . . , = × ℎ + 1 − 2 1 , 2 2 1 , 2 = =1 × =1 × 2 2 ℎ + 1 − 2 1 , 2 1 , 2 1 , 2 10 = ℎ + 1 − 2 1 , 2 8 + =1 6 0 =1 y + − 4 = 2 0 Zagreb June 5, 2014 2 4 6 x 8 10 Limit Theory Remaining steps: = × =1 ℎ + 1 − 2 1 , 2 2 1 , 2 i. Show var − ii. Let be an iid sequence with = whose sum has characteristic function iii. → 0. . Show → exp 2 2 − 2 − → 0. Intuition. (i) The sets ∖ are small by proper choice of and . (ii) Use a Lynapounov CLT (have a triangular array). (iii) Use a Bernstein argument (see next page). Zagreb June 5, 2014 . Limit Theory Useful identity: =1 =1 − = =1 1 ⋯ −1 − +1 ⋯ | − ()| = | =1 −1 | − =1 ( − ) = | =1 =1 ≤ ≤ ≤ where , =+1 −1 |cov( , )| =1 =1 −1 |(cov( , ) =1 =1 =1 4 (∪−1 =1 ),( ) | (by indep of ) | ℎ is a strong mixing bounding function that is based on the separation h between two sets U and V with cardinality r and s. Zagreb June 5, 2014 Strong mixing coefficients Strong mixing coefficients: Let be a stationary random field on ℝ2 . Then the mixing coefficients are defined by , ℎ = sup ∩ − , E1 ,E2 where the sup is taking over all sets ∈ 1 , ∈ 2 , with (1 ) ≤ , (2 ) ≤ , and 1 , 2 ≥ ℎ. Proposition (Li, Genton, Sherman (2008), Ibragimov and Linnik (1971)): Let and be closed and connected sets such that ≤ , ≤ and , ≥ ℎ. If and are rvs measurable wrt () and , respectively, and bded by 1, then cov , ≤ 4, ℎ 16, ℎ , . Zagreb June 5, 2014 Limit Theory Mixing condition:sup ℎ / = (ℎ− ) for some > 2. s Returning to calculations: | − ()| = ≤ 16 =1 , )| |cov( =1 −1 =1 (∪−1 =1 ),( ) 16 ∪=1 − ≤ =1 ≤ 162 − ≤ 2 2 − =1 4−2− = → 0 4 − 2 − < 0 . Zagreb June 5, 2014 Simulations of spatial extremogram Extremogram for one realization of B-R process (function of level) Note: black dots = true; blue bands are permutation bounds Lattice Zagreb June 5, 2014 Non-Lattice 37 Simulations of spatial extremogram Max-moving average (MMA) process: Let be an iid sequence of Frechet rvs and set = max2 − , ∈ℤ where ≥ 0, < ∞. MMA(1): = max − . ≤1 Extremogram: 1, 2/5 ℎ = lim ℎ > 0 > = 1/5 0 Zagreb June 5, 2014 ℎ ℎ ℎ ℎ = 0, = 1, √2, = 2, > 2. 38 Simulations of spatial extremogram Box-plots based on 1000 (100) replications of MMA(1) (left) and BR (right) Lattice Non-lattice; = 1/ log (left) = 5/ log (right) Zagreb June 5, 2014 39 Space-Time Extremogram Extremogram. , (ℎ, ) = lim (( + ℎ, + ) | (, ) ) For the special case, = = (1, ∞), ℎ, = lim + ℎ, + > , > ) For the Brown-Resnick process described earlier ℎ, = 2(1 − Φ 1 ℎ1 + 2 2 , we find that 1 2 2 log(Φ−1 (1 − ℎ, 0 ) = log 1 + 1 log ℎ and 1 2 2 log(Φ−1 (1 − 0, ) = log 2 + 2 log Zagreb June 5, 2014 40 Inference for Brown-Resnick Process Semi-parametric: Use nonparametric estimates of the extremogram and then regress function of extremogram on the lag. Regress: 1 2 2 log(Φ−1 (1 − ℎ, 0 )) on 1 and log ℎ 1 2 2log(Φ−1 (1 − 0, )) on 1 and log , The intercepts and slopes become the respective estimates of log and . Asymptotic properties of spatial extremogram derived by Cho et al. (2013). Pairwise likelihood: Maximize the pairwise likelihood given by PL ( , ) log f ( ( s i , t k ), ( s j , t l ); , ) ( k , l ) :| t k t l | D i ( i , j ) :| s j s i | D i Zagreb June 5, 2014 41 Inference for Brown-Resnick Process Semi-parametric: Use nonparametric estimates of the extremogram and then regress function of extremogram on the lag. Regress: 1 2 2 log(Φ−1 (1 − ℎ, 0 ) on 1 and log ℎ 1 2 2log(Φ−1 (1 − 0, ) on 1 and log , The intercepts and slopes become the respective estimates of log and . Asymptotic properties of spatial extremogram derived by Cho et al. (2013). Pairwise likelihood: Maximize the pairwise likelihood given by PL ( , ) log f ( ( s i , t k ), ( s j , t l ); , ) ( k , l ) :| t k t l | D i ( i , j ) :| s j s i | D i Zagreb June 5, 2014 42 Data Example: extreme rainfall in Florida Zagreb June 5, 2014 43 Data Example: extreme rainfall in Florida Radar data: Rainfall in inches measured in 15-minutes intervals at points of a spatial 2x2km grid. Region: 120x120km, results in 60x60=3600 measurement points in space. Take only wet season (June-September). Block maxima in space: Subdivide in 10x10km squares, take maxima of rainfall over 25 locations in each square. This results in 12x12=144 spatial maxima. Temporal domain: Analyze daily maxima and hourly accumulated rainfall observations. Fit our extremal space-time model to daily/hourly maxima. Zagreb June 5, 2014 44 Data Example: extreme rainfall in Florida Hourly accumulated rainfall time series in the wet season 2002 at 2 locations. Zagreb June 5, 2014 45 Data Example: extreme rainfall in Florida Hourly accumulated rainfall fields for four time points. Zagreb June 5, 2014 46 Data Example: extreme rainfall in Florida Empirical extremogram in space (left) and time (right) Zagreb June 5, 2014 47 Data Example: extreme rainfall in Florida Empirical extremogram in space (left) and time (right): spatial indep for lags > 4; temporal indep for lags > 6. Zagreb June 5, 2014 48 Data Example: extreme rainfall in Florida Empirical extremogram in space (left) and time (right) Zagreb June 5, 2014 49 Data Example: extreme rainfall in Florida Computing conditional return maps. Estimate , such that , > , ∗ , ∗ > ∗ = , where ∗ satisfies ∗ , ∗ > ∗ = ∗ is pre-assigned. A straightforward calculation shows that , must solve, 1 1 = 1 − ∗ exp − , Zagreb June 5, 2014 1 + ∗ ( , , 1 − ∗ ) , 50 100-hour return maps ( = .01): ∗ = 5,6 , time lags = 0,2,4,6 hours (left to right on top and then right to left on bottom), quantiles in inches. Zagreb June 5, 2014 51 Estimation—composite likelihood approach For dependent data, it is often infeasible to compute the exact likelihood based on some model. An alternative is to combine likelihoods based on subsets of the data. To fix ideas, consider the following data/model setup: (Here we have already assumed that the data has been transformed to a stationary process with unit Frechet marginals.) Data: Y(s1), …, Y(sN) (field sampled at locations s1, …, sN ) Model: max-stable model defined via the limit process maxj=1,…,nYn(j)(s) →d X(s), • Yn(s) = Y(s /(log n)1/b)) = 1/log(F(Z(s/(log n)1/b)) • Z(s) is a GP with correlation function r(|s-t|) = exp{-|s-t|b/f} Zagreb June 5, 2014 52 Estimation—composite likelihood approach Bivariate likelihood: For two locations si and sj, denote the pairwise likelihood by f(y(si), y(sj); di,j ) = ∂2/(∂x∂y) F(Y(si) ≤ x, Y (sj) ≤ x) where F is the CDF F(Y(si) ≤ x, Y (sj) ≤ x) = exp{-(x-1F(log(y/x)/(2d) + d) + y-1F(log(x/y)/(2d) + d))}, and di,j |si – sj|b/f is a function of the parameters b and f. Pairwise log-likelihood: N PL (f , b ) log i 1 i j Zagreb June 5, 2014 N f ( y ( s i ), y ( s j ); d i , j ) j 1 53 Estimation—composite likelihood approach Potential drawbacks in using all pairs: • Still may be computationally intense with N2 terms in sum. • Lack of consistency (especially if the process has long memory) • Can experience huge loss in efficiency. Remedy? Use only nearest neighbors in computing the pairwise likelihood terms. 1. Use only neighbors within distance L. The new criterion becomes N PL (f , b ) log f ( y ( s i ), y ( s j ); d i , j ) i 1 j :| s j s i | L 2. Use only nearest K neighbors. If Di denotes the distance of kth nearest neighbor to si , then criterion is N PL (f , b ) Zagreb June 5, 2014 log i 1 j :| s j s i | D i f ( y ( s i ), y ( s j ); d i , j ) 54 Simulation Examples Simulation setup: • Simulate 1600 points of a spatial (max-stable) process Y(s) on a grid of 40x40 in the plane. • Choose a distance L or number of neighbors K (usually 9, 15, 25) • Maximize N PL (f , b ) log f ( y ( s i ), y ( s j ); d i , j ) i 1 j :| s j s i | D i with respect to f and b • Calculate summary dependence statistics: r(l; |s-t|) = l-1(1 – l(1l)V(l, 1l; d)), where V(l,1-l; d) = l-1F(log ((1-l)/l) /(2d)+d) + (1-l)-1F(log (l/(1-l))/(2d)+d) d = |s-t|b/f Zagreb June 5, 2014 55 Simulation Examples Process: Limit process d = |s-t|, 1 , .2; 9 nearest neighbors (; | − ) = lim > 1 − →∞ Zagreb June 5, 2014 () > ) 56 Simulation Examples Process: Limit process d = |s-t|, 1 , .2; 25 nearest neighbors (; | − ) = lim > 1 − →∞ Zagreb June 5, 2014 () > ) 57 Illustration with French Precipitation Data Data from Naveau et al. (2009). Precipitation in Bourgogne region of France; 51 year maxima of daily precipitation. Data has been adjusted for seasonality and orographic effects. Zagreb June 5, 2014 60 Illustration with French Precipitation Data Estimated spatial correlation function before transformation to unit Frechet. 0.4 0.0 0.2 cor 0.6 0.8 1.0 Correlation function 0 Zagreb June 5, 2014 200 400 600 lag 800 61 Illustration with French Precipitation Data After transforming the data to unit Frechet marginals, we estimated and −1 =: using pairwise likelihood (nearest neighbors with K = 25). 1.143; f 185.1 (if constrained to 1, then f 89.2) 1.0 1.0 Correlation function 0.8 0.8 Green: 1.143 ; f 185.1 Black: 1 ; f 139.6 (MLE) Red: 1.17 ; f 147.8 0.4 0.4 cor cor 0.6 0.6 Blue: 1 ; f 89.2 0.2 0.2 (Matern) 0.0 0.0 (MLE based on GP likelihood) 0 0 Zagreb June 5, 2014 200 200 400 600 400 600 lag 800 800 62 Illustration with French Precipitation Data Plot of r(l; |s-t|) = limn→∞ P(Yn(s) > n(1-l) | Yn(t) > nl) Zagreb June 5, 2014 63 Illustration with French Precipitation Data Plot of r(l; |s-t|) = limn→∞ P(Yn(s) > n(1-l) | Yn(t) > nl) Zagreb June 5, 2014 64 Illustration with French Precipitation Data Results from random permutations of data—just for fun. Since random permutations have no spatial dependence, r should be 0 for |s-t| >0. Zagreb June 5, 2014 65