#### Transcript Asymptotics for the Spatial Extremogram

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 𝑘𝑛 ≤ 16𝑖𝑛2𝛼 𝑛−𝜖𝜂 ≤ 𝐶𝑘𝑛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