### Asymptotics for the Spatial Extremogram

```Asymptotics for the Spatial Extremogram
Richard A. Davis*
Columbia University
Collaborators:
Yongbum Cho, Columbia University
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
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
```