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)/(2d) + d) + y-1F(log(x/y)/(2d) + 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(1l)V(l, 1l; d)), where
V(l,1-l; d)
= l-1F(log ((1-l)/l) /(2d)+d) + (1-l)-1F(log (l/(1-l))/(2d)+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