Submitted to Bernoulli
arXiv: 1609.04967

arXiv:1609.04967v3 [stat.ME] 5 Nov 2017

Semiparametric estimation for isotropic
max-stable space-time processes
SVEN BUHL1,* RICHARD A. DAVIS2,** CLAUDIA
KLÜPPELBERG1,† and CHRISTINA STEINKOHL1,‡
1

Center for Mathematical Sciences and TUM Institute of Advanced Study, Technische
Universität München, Boltzmannstr. 3, 85748 Garching, Germany.
E-mail: * sven.buhl@tum.de; † cklu@tum.de; ‡ christina.steinkohl@gmail.com, url:
http://www.statistics.ma.tum.de
2

Department of Statistics, Columbia University, 1255 Amsterdam Avenue, New York, NY
10027, USA. E-mail: ** rdavis@stat.columbia.edu, url:
http://www.stat.columbia.edu/~rdavis
Abstract
Max-stable space-time processes have been developed to study extremal dependence in spacetime data. We propose a semiparametric estimation procedure based on a closed form expression
of the extremogram to estimate the parameters in a max-stable space-time process. We establish the asymptotic properties of the resulting parameter estimates and propose subsampling
procedures to obtain asymptotically correct confidence intervals. A simulation study shows that
the proposed procedure works well for moderate sample sizes. Finally, we apply this estimation
procedure to fitting a max-stable model to radar rainfall measurements in a region in Florida.
Keywords: Brown-Resnick process, extremogram, max-stable process, semiparametric estimation, space-time process, subsampling, mixing.

1. Introduction
Max-stable processes are a natural extension of the generalized extreme value distributions to infinite dimensions and provide a useful framework for modeling extremal
dependence in continuous time or space. In this paper we focus on the max-stable BrownResnick process, which was introduced in a time series framework in Brown and Resnick
[3], in a spatial setting in Kabluchko et al. [24], and extended to a space-time setting in
Davis et al. [10].
In the literature, various max-stable models and estimation procedures have been proposed for extremal data. For the Brown-Resnick process with parametrized dependence
structure, inference has been based on composite likelihood methods. In particular, pairwise likelihood estimation has been found useful to estimate parameters in a max-stable
process. A description of this method can be found in Padoan et al. [26] for the spatial
setting, and Huser and Davison [23] in a space-time setting. Asymptotic results for pairwise likelihood estimates and detailed analyses in the space-time setting for the model
1
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

2

S. Buhl et al.

analysed in this paper are given in Davis et al. [11]. Unfortunately, parameter estimation
using composite likelihood methods can be laborious, since the computation and subsequent optimization of the objective function is time-consuming. Also the choice of good
initial values for the optimization of the composite likelihood is essential.
In this paper we introduce a new semiparametric estimation procedure which is based
on the extremogram as a natural extremal analog of the correlation function for stationary processes. The extremogram was introduced in Davis and Mikosch [9] for time
series (also in Fasen et al. [19]), and they show consistency and asymptotic normality
of an empirical extremogram estimate under weak mixing conditions. The empirical extremogram and its asymptotic properties in a spatial setting have been investigated in
Buhl and Klüppelberg [6] and Cho et al. [8]. It can serve as a useful graphical tool for assessing extremal dependence structures in spatial and space-time processes that provides
clues about potential parametric models to consider, a critical step in the model building
paradigm. For example, compatibility with various process assumptions such as isotropy
and stationarity (see Buhl and Klüppelberg [5] and Davis et al. [11] for some examples),
can be assessed by examining invariance of the empirical extremogram when computed
over specially chosen subsets of the data. Ultimately, a number of families of proposed
parametric models are often fitted before deciding on a particular class of models. Therefore it is of interest to be able to not only have a procedure that can compute estimates
rapidly, but also to serve as a check on the efficacy of model choice. Using the estimation
procedure described in this paper allows one to provide initial parameter estimates that
can be used as initial values in more refined procedures, such as composite likelihood.
Our semiparametric estimation method assumes the same dependence structure for the
Brown-Resnick space-time process as in [10, 11], yielding a closed form expression of the
extremogram containing the parameters of interest. We first estimate the extremogram
nonparametrically by its empirical version, where we separate space and time. Weighted
linear regression is then applied in order to produce parameter estimates. Asymptotic
normality of these semiparametric estimates requires asymptotic normality of the empirical extremogram, and we apply the CLT with mixing conditions as provided in [6].
Since the dependence structure given by spatial and temporal parameters is isotropic in
space and additively separable, the proofs of the asymptotic properties of semiparametric spatial and temporal parameter estimates are analogous, and we present the details
on the spatial parameters only, referring to Buhl [4], Chapter 3, for details about the
asymptotic properties of the semiparametric temporal parameter estimates. The rate of
convergence can be improved by a bias correction term, which we explain in detail.
In a second step we prove then asymptotic normality of the weighted least squares
parameter estimates, where constrained optimization has to be applied, since one of
the parameters has bounded support. Also the limit law differs depending whether the
parameter lies on the boundary or not. Since the asymptotic covariance matrix in the
normal limit is difficult to access, we apply subsampling procedures to obtain pointwise
confidence intervals for the parameters, also taking care of the different limit laws.
The rate of convergence of the semiparametric estimates is not as fast as that of a
fully parametric procedure such as pairwise likelihood estimation whose rate equals the
square root of the effective sample size. However, it is known that such estimates may
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

3

be inefficient and even not consistent if the model is slightly misspecified. On the other
hand, the semiparametric estimates are often unaffected by slight deviations in the model.
This is illustrated in Section 5, where data are generated from a Brown-Resnick process,
but with observational noise. The semiparametric estimates outperform pairwise likelihood estimates in this case. On the other hand, the semiparametric estimates perform
admirably well relative to the pairwise likelihood estimates even when the underlying
process is in fact a Brown-Resnick process.
Our paper is organized as follows. Section 2 defines the isotropic Brown-Resnick process with its choice of dependence function used throughout for modelling extremes observed in space and time. The extremogram is introduced and its parametric form for
our model is given. Based on gridded data, the nonparametric extremogram estimation
is derived. Asymptotic normality of the parameter estimates is established in Section 3.
Section 3.1 is dedicated to the asymptotic normality of the empirical extremogram and
its bias correction; Section 3.2 deals with the asymptotic properties of the parameter
estimates. In Section 4 we explain the subsampling procedure. We test our new semiparametric estimation procedure in a simulation study in Section 5 and compare with
pairwise likelihood estimation, both when applied to data generated by a Brown-Resnick
process and when the data are affected by observational noise. The paper concludes
with an analysis of daily rainfall maxima in a region in Florida in Section 6, where we
also compare the semiparametric estimates with previously obtained pairwise likelihood
estimates. Some auxiliary results are summarized in an appendix.

2. Model description and semiparametric estimates
Throughout the paper we consider a strictly stationary Brown-Resnick process in space
and time with representation
η(s, t) =

∞ n
o
_
ξj eWj (s,t)−δ(ksk,t) ,

(s, t) ∈ R2 × [0, ∞),

(2.1)

j=1

where k · k is the Euclidean norm on R2 , {ξj : j ∈ N} are points of a Poisson process
on [0, ∞) with intensity ξ −2 dξ and the dependence function δ is nonnegative and conditionally negative definite; i.e., for every m ∈ N and every (s(1) , t(1) ), . . . , (s(m) , t(m) ) ∈
R2 × [0, ∞), it holds that
m X
m
X

ai aj δ(ks(i) − s(j) k, |t(i) − t(j) |) ≤ 0

i=1 j=1

for all a1 , . . . , am ∈ R summing up to 0. The processes {Wj (s, t) : s ∈ R2 , t ∈ [0, ∞)}
are independent replicates of a Gaussian process {W (s, t) : s ∈ R2 , t ∈ [0, ∞)} with
stationary increments, W (0, 0) = 0, E[W (s, t)] = 0 and covariance function
Cov[W (s(1) , t(1) ),W (s(2) , t(2) )]
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

4

S. Buhl et al.
= δ(ks(1) k, t(1) ) + δ(ks(2) k, t(2) ) − δ(ks(1) − s(2) k, |t(1) − t(2) |).

Representation (2.1) goes back to de Haan [13], Giné et al. [21] and Kabluchko et al.
[24]. All finite-dimensional distributions are multivariate extreme value distributions
with standard unit Fréchet margins, hence they are in particular multivariate regularly varying. Furthermore, they are perfectly characterized by the dependence function δ, which is termed the semivariogram of the process {W (s, t)} in geostatistics: For
(s(1) , t(1) ), (s(2) , t(2) ) ∈ R2 × [0, ∞), it is given by
Var[W (s(1) , t(1) ) − W (s(2) , t(2) )] = 2δ(ks(1) − s(2) k, |t(1) − t(2) |).
Since we assume δ to depend only on the norm of s(1) − s(2) , the associated process is
(spatially) isotropic.
In this paper we assume the dependence function δ to be given for v, u ≥ 0 by
δ(v, u) = 2θ1 v α1 + 2θ2 uα2 ,

(2.2)

where 0 < α1 , α2 ≤ 2 and θ1 , θ2 > 0. This is the fractional class frequently used for
dependence modelling, and here defined with respect to space and time.
The bivariate distribution function of (η(0, 0), η(h, u)) is given for x1 , x2 > 0 by
!
(
r
δ(khk, |u|)
log(x2 /x1 )
1
+
F (x1 , x2 ) = exp − Φ p
x1
2
2δ(khk, |u|)
!)
r
1
log(x1 /x2 )
δ(khk, |u|)
− Φ p
+
,
(2.3)
x2
2
2δ(khk, |u|)
where Φ denotes the standard normal distribution function (cf. Davis et al. [10]).
The parameters of interest are contained in the dependence function δ. We refer to
(θ1 , α1 ) as the spatial parameter and to (θ2 , α2 ) as the temporal parameter. From the
bivariate distribution function in (2.3), the pairwise density can be derived and pairwise
likelihood methods can be used to estimate the parameters; cf. Davis et al. [11], Huser
and Davison [23] and Padoan et al. [26]. Full likelihood inference is typically hardly
tractable in a general multidimensional setting, as the number of terms occurring in
the likelihood explode. More recently, however, parametric inference methods based on
higher-dimensional margins have been proposed that work in specific scenarios, see for
instance Genton et al. [20], who use triplewise instead of pairwise likelihood, Engelke
et al. [18], who propose a threshold-based approach, or Thibaud and Opitz [31] and
Wadsworth and Tawn [32], who use a censoring scheme for bias reduction.
In the following we introduce an alternative estimation approach, which is based on
a closed form expression of the extremogram. The latter was introduced for time series
by Davis and Mikosch [9] and for spatial and space-time processes by Cho et al. [8] and
Steinkohl [30], respectively, and can be regarded as a correlogram for extreme events.
In this paper we consider an isotropic Brown-Resnick process as a regularly varying
stochastic processes {η(s) : s ∈ Rd } where d = 1 corresponds to a time series and d = 2
to a spatial process, such that d = 3 holds for the space-time process.
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

5

More precisely, we consider strictly stationary regularly varying processes {η(s) : s ∈
Rd } for d ∈ N, where all finite-dimensional distributions are regularly varying (cf. Hult
and Lindskog [22] for definitions and results in a general framework and Resnick [28]
for details about multivariate regular variation). Throughout, f (n) ∼ g(n) means that
(n)
limn→∞ fg(n)
= 1. As a prerequisite, we define for every finite set I ⊂ Rd with cardinality
|I| the vector
|
ηI := (η(s) : s ∈ I) .
Definition 2.1 (Regularly varying stochastic process). A strictly stationary stochastic
process {η(s) : s ∈ Rd } is called regularly varying, if there exists some normalizing
sequence 0 < an → ∞ such that P(|η(0)| > an ) ∼ n−d as n → ∞, and if for every finite
set I ⊂ Rd ,

η
v
I
∈ · → µI (·), n → ∞,
(2.4)
nd P
an
for some non-null Radon measure µI on the Borel sets in R
µI (xC) = x−β µI (C),
|I|

|I|

\{0}. In that case,

x > 0,
v

for every Borel set C in R \{0}. The notation → stands for vague convergence, and
β > 0 is called the index of regular variation.
For every s ∈ Rd and I = {s} we set µ{s} (·) = µ{0} (·) =: µ(·), which is justified by
stationarity.
Assuming strict stationarity and spatial isotropy of a regularly varying space-time
process {η(s, t) : s ∈ R2 , t ∈ [0, ∞)} we can define its extremogram at two points (s1 , t1 )
and (s2 , t2 ) only in terms of the spatial and temporal lags v := ks1 −s2 k and u := |t1 −t2 |.
Definition 2.2 (The extremogram). For a regularly varying strictly stationary isotropic
space-time process {η(s, t) : (s, t) ∈ R2 × [0, ∞)} we define the space-time extremogram
for two µ-continuous Borel sets A and B in R\{0} (i.e. µ(∂A) = µ(∂B) = 0) such that
µ(A) > 0 by
ρAB (v, u) = lim

n→∞

P (η(s1 , t1 )/an ∈ A, η(s2 , t2 )/an ∈ B)
,
P (η(s1 , t1 )/an ∈ A)

(2.5)

where v = ks1 − s2 k and u = |t1 − t2 |.
Setting A = B = (1, ∞), we rediscover the tail dependence coefficient χ(v, u) =
ρ(1,∞)(1,∞) (v, u). For the isotropic Brown-Resnick process there is a closed form expression for χ(v, u), which is the basis for our estimation procedure.

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

6

S. Buhl et al.

Lemma 2.3 (Davis et al. [10], equation (3.1)). Let η be the strictly stationary isotropic
Brown-Resnick process in R2 × [0, ∞) as defined in (2.1) with dependence function given
in (2.2). Then for A = B = (1, ∞) the extremogram of η is given by

r 1

p

(2.6)
χ(v, u) = 2 1 − Φ
δ(v, u) = 2 1 − Φ( θ1 v α1 + θ2 uα2 ) , v, u ≥ 0.
2
Solving equation (2.6) for δ(v, u) leads to

 2
δ(v, u)
1
.
= θ1 v α1 + θ2 uα2 = Φ−1 1 − χ(v, u)
2
2

(2.7)

For temporal lag 0 and taking the logarithm on both sides we have


1
2 log Φ−1 1 − χ(v, 0) = log(θ1 ) + α1 log v.
2
In the same way, we obtain


1
2 log Φ−1 1 − χ(0, u) = log(θ2 ) + α2 log u.
2
These equations are the basis for parameter estimates. We replace the extremogram on
the left hand side in both of these equations by nonparametric estimates at different
lags. Then we use constrained weighted least squares estimation in a linear regression
framework to obtain parameter estimates.
The estimation procedure is based on the following observation scheme for the spacetime data.
Condition 2.4. (1) The locations lie on a regular 2-dimensional grid


Sn = (i1 , i2 ) : i1 , i2 ∈ {1, . . . , n} = si : i = 1, . . . , n2 .
(2) The time points are equidistant, given by the set {t1 , . . . , tT }.
Remark 2.1. The assumption of a regular grid can be relaxed in various ways. A
simple, but notationally more involved extension is the generalization to rectangular
grids, cf. Buhl and Klüppelberg [6], Section 3. Furthermore, it is possible to assume that
the observation area consists of random locations given by points of a Poisson process, see
for instance Cho et al. [8], Section 2.3, or Steinkohl [30], Section 4.5.2. Also deterministic,
but irregularly spaced locations, could be considered as treated in [30] in Section 4.5.1 in
the context of pairwise likelihood estimation. In order to make our method transparent
we focus on observations on a regular grid.

The following scheme provides the semiparametric estimation procedure in detail.
Denote by V ⊂ [0, ∞) and U ⊂ [0, ∞) finite sets of spatial and temporal lags, on which

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

7

the estimation is based. We denote by “lag” the norm or absolute value of the difference of two spatial locations or two time points, respectively. Concerning the choice of V
and U, we generally include those lags which show clear extremal dependence between
locations or time points. Larger lags should not be considered, since they may introduce
a bias in the least squares estimates, similarly as in pairwise likelihood estimation; cf.
Buhl and Klüppelberg [5], Section 5.3. One way to determine the range of clear extremal
dependence are permutation tests, which we describe at the end of Section 6.
(1) Nonparametric estimates for the extremogram:
Summarize all pairs of Sn which give rise to the same spatial lag v ∈ V into
N (v) = {(i, j) ∈ {1, . . . , n2 }2 : ksi − sj k = v}.
For all t ∈ {t1 , . . . , tT } estimate the spatial extremogram by
2

n
P
1
|N (v)| i=1

χ
b(t) (v, 0) =

2

n
P

1{η(si ,t)>q,η(sj ,t)>q}

j=1
ksi −sj k=v

,

n2
1 P
1{η(si ,t)>q}
n2 i=1

v ∈ V,

(2.8)

where q is a large quantile (to be specified) of the standard unit Frechét distribution.
For all s ∈ Sn estimate the temporal extremogram by
1
T −u
(s)

χ
b

(0, u) =

TP
−u

1{η(s,tk )>q,η(s,tk +u)>q}

k=1
1
T

T
P

,

u∈U,

(2.9)

1{η(s,tk )>q}

k=1

where again q is a large (possibly different) quantile of the standard unit Frechét distribution
(2) The overall “spatial” and “temporal” extremogram estimates are defined as averages
over the temporal and spatial locations, respectively; i.e.,
χ
b(v, 0)

=

T
1 X (tk )
χ
b (v, 0),
T

v ∈ V,

(2.10)

k=1
2

χ
b(0, u)

=

n
1 X (si )
χ
b (0, u),
n2 i=1

u ∈ U.

(2.11)

(3) Parameter estimates for θ1 , α1 , θ2 and α2 are found by using constrained weighted
least squares estimation:
 

X 

 2
1
θb1
= arg min
wv 2 log Φ−1 1 − χ
b(v, 0) − log(θ1 ) + α1 log(v)
, (2.12)
α
b1
2
θ1 ,α1 >0
α1 ∈(0,2]

v∈V

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

8

S. Buhl et al.


θb2
α
b2


= arg min
θ2 ,α2 >0

α2 ∈(0,2]

X
u∈U




 2
1
wu 2 log Φ−1 (1 − χ
b(0, u)) − log(θ2 ) + α2 log(u)
, (2.13)
2

with weights wu > 0 and wv > 0.
We call the estimates (θb1 , α
b1 ) and (θb2 , α
b2 ) weighted least squares estimates (WLSE).
This approach bears similarity with that proposed by Einmahl et al. [17], who suggest
semiparametric weighted least squares estimation of the parameters of parametric models
of the stable tail dependence function based on iid random vector observations.

3. Estimation of the isotropic Brown-Resnick process
In this section we investigate aymptotic properties of the WLSE (θb1 , α
b1 ) and (θb2 , α
b2 ).
For a central limit theorem of the extremogram we need a sufficiently precise estimate
for the extremogram (2.6), which we give now.
Lemma 3.1. Let s, h ∈ R2 and t ∈ [0, ∞). For every sequence an → ∞ we have
P(η(s, t) > an , η(s + h, t) > an )
P(η(s, t) > an )
h 1

i
χ(khk, 0) − 2 χ(khk, 0) − 1 (1 + o(1)).
=χ(khk, 0) +
2an
Lemma 3.1 is a direct application of Lemma A.1(b) of Buhl and Klüppelberg [6] for
A = B = (1, ∞) and their equation (A.4). This applies since {η(s, t) : s ∈ R2 } has
finite-dimensional standard unit Fréchet marginal distributions. In order to satisfy the
condition P(|η(0, 0)| > an ) ∼ n−2 as n → ∞ from Definition 2.1, we choose an ∼ n2 .
Since the WLSE are functions of the spatial and temporal extremograms, we first derive the asymptotic properties of χ
b(t) and χ
b(s) for a fixed time point t and a fixed location
s, respectively. Sections 3.1 and 3.2 focus on the spatial parameters. The√corresponding
results for the temporal case can be derived similarly by replacing n with T and can be
found with full details in Buhl [4], Chapter 3. We use several results for the extremogram
provided in Appendix A and in Buhl and Klüppelberg [6].

3.1. Asymptotics of the empirical spatial extremogram
We prove a central limit theorem for the empirical spatial extremogram of the BrownResnick process (2.1) based on a finite set of spatial lags
V = {v1 , . . . , vp },
which show clear extremal dependence as explained in Section 2. First we show that the
empirical extremogram centred by the pre-asymptotic version is asymptotically normal.
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

9

Theorem 3.2.  For a fixed time point t ∈ {t1 , . . . , tT }, consider the spatial BrownResnick process η(s, t) : s ∈ R2 as defined in (2.1) with dependence function given in
(2.2). Set mn = nβ1 for β1 ∈ (0, 1/2). Then the empirical spatial extremogram χ
b(t) (v, 0)
2
defined in (2.8) with the quantile q = mn satisfies

n
d
(iso)
χ
b(t) (v, 0) − χn (v, 0) v∈V → N (0, Π1 ),
mn

n → ∞,

(iso)

where the covariance matrix Π1
is specified in equation (3.6) below, and χn is the
pre-asymptotic spatial extremogram,
χn (v, 0) =

P(η(0, 0) > m2n , η(h, 0) > m2n )
,
P(η(0, 0) > m2n )

v = khk ∈ V .

(3.1)

Proof. As η(0, t) has standard unit Fréchet marginal distributions, we can choose amn =
m2n by Definition 2.1 of regular variation.
We apply Theorem 4.2 of Buhl and Klüppelberg [6] by verifying conditions (M1)-(M4)
of that theorem for {η(s, t) : s ∈ R2 }, d = 2, and A = B = (1, ∞). Condition (M1) is
satisfied by equation (A.2).
To show conditions (M2)-(M4) we choose sequences mn = nβ1 and rn = nβ2 for 0 <
β1 < 1/2 and 0 < β2 < β1 . For this choice mn and rn increase to infinity with mn = o(n)
and rn = o(mn ) as required.
Condition (M2); i.e., m2n rn2 /n = n2(β1 +β2 )−1 → 0 holds if and only if β2 ∈ (0, min{β1 , (1/2−
β1 )}).
We now show condition (M3). Choose γ > 0, such that all lags in V lie in B(0, γ) := {s ∈
Z2 : ksk ≤ γ}. Denote by B(h, γ) := {s ∈ Z2 : ks − hk ≤ γ} = h + B(0, γ) for h ∈ R2 .
For ε > 0, like in Example 4.6 of Buhl and Klüppelberg [6], we have for s, s0 ∈ R2 by a
Taylor expansion,
P(η(s, t) > εm2n , η(s0 , t) > εm2n )
= 1 − 2P(η(0, 0) ≤ εm2n ) + P(η(s, t) ≤ εm2n , η(s0 , t) ≤ εm2n )
n 1o
n 2 − χ(ks − s0 k, 0) o
+ exp −
= 1 − 2 exp −
x
εm2n


1
1
=
χ(ks − s0 k, 0) + O
, n → ∞.
2
εmn
m4n
Therefore, for khk ≥ 2γ,
P( max η(s, t) > εm2n , 0 max η(s0 , t) > εm2n )
s∈B(0,γ)
s ∈B(h,γ)
X
X
≤
P(η(s, t) > εm2n , η(s0 , t) > εm2n )
s∈B(0,γ) s0 ∈B(h,γ)

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

10

S. Buhl et al.
=

X

X

s∈B(0,γ) s0 ∈B(h,γ)

≤

 1 o
n 1
0
χ(ks
−
s
k,
0)
+
O
εm2n
m4n

 1 
p

2|B(0, γ)|2
α1 + O
,
1
−
Φ(
θ
(khk
−
2γ)
1
εm2n
m4n

(3.2)

as n → ∞, where we have used (2.6). Summarize V := {v = khk : h ∈ Z2 } and note
that |{h ∈ Z2 : khk = v}| = O(v). Therefore, for k ≥ 2γ,


X
Lmn := lim sup m2n
P
max η(s, t) > εm2n , 0 max η(s0 , t) > εm2n
n→∞

h∈Z2
k<khk≤rn



2

≤ 2|B(0, γ)| lim sup
n→∞

≤ C1 lim sup
n→∞

s ∈B(h,γ)

s∈B(0,γ)

o
 r 2 
p
n
α
1
(1 − Φ( θ1 (khk − 2γ) )) + O
ε
mn

n1

X
h∈Z2
k<khk≤rn

X nv
p
o
2 1 − Φ( θ1 (v − 2γ)α1 ) ,
ε
v∈V :
k<v≤rn

for some constant C1 > 0. For the term O((rn /mn )2 ) we use that rn /mn → 0. From
Lemma A.3 and the fact that 1 − Φ(x) ≤ exp{−x2 /2} for x > 0, we find for C2 > 0,
 1
Lmn ≤C2 k 2 exp − θ1 (k − 2γ)α1 .
2
Since α1 > 0, the right hand side converges to 0 as k → ∞ ensuring condition (M3).
Now we turn to the mixing conditions (M4).
We start with (M4i). With V as before, and with equation (A.2), we estimate, recalling
from above that the number of lags khk = v is of oder O(v),
X
X
X
α1
m2n
α1,1 (khk) ≤ C1 m2n
v α1,1 (v) ≤ 4C1 m2n
v e−θ1 v /2 .
h∈Z2 :khk>rn

v∈V :v>rn

v∈V :v>rn

By Lemma A.3 we find
X
α1
α1
α1 β2
/2
m2n
v e−θ1 v /2 ≤ cm2n rn2 e−θ1 rn /2 = cm2n rn2 e−θ1 n
→ 0,

n → ∞.

v∈V :v>rn

By the same arguments condition (M4ii) is satisfied.
Condition (M4iii) holds by equation (A.2), since
α1

mn n α1,n2 (rn ) ≤ 4n3 mn e−θ1 rn

/2

→ 0,

n → ∞.

For the specification of the asymptotic covariance matrix we apply Theorem 4.2 of [6] for
the isotropic case, where each spatial lag vi arises from a set of different vectors h, all
with same Euclidean norm vi . For i ∈ {1, . . . , p} such that vi ∈ V, we summarize these
into
(i)
(i)
L(vi ) := {h ∈ Z2 : khk = vi } = {h1 , . . . , h`i },
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

11

where `i := |L(vi )|. Based on the preceding steps of the proof, we conclude from that
theorem that
|
n
d
(i)
(i)
(i)
(i)
(space)
χ
b(t) (h1 , 0) − χn (h1 , 0), . . . , χ
b(t) (h`i , 0) − χn (h`i , 0) i=1,...,p → N (0, Π1
),
mn
(space)

where Π1
is specified in equation (4.3)-(4.6) of [6]. Furthermore, χ
b(t) (h, 0) is the
empirical extremogram for each vector h as specified above.
Define N (h) := {(i, j) ∈ {1, . . . , n2 } : si −sj = h}, then the numerator in (2.8) normalizes
by |N (h)| (instead of |N (v)|) and the sum
P runs over si − sj = h vector-wise (instead
of equality in norm). Hence, |N (vi )| = h∈L(vi ) |N (h)|. Isotropy implies for the preasymptotic extremogram that χn (vi , 0) = χn (h, 0) for all h ∈ L(vi ), such that
χn (vi , 0) =

X
h∈L(vi )

|N (h)|
χn (vi , 0) =
|N (vi )|

X
h∈L(vi )

|N (h)|
χn (h, 0)
|N (vi )|

(3.3)

|N (h)| (t)
χ
b (h, 0).
|N (vi )|

(3.4)

as well as, by the definition of the estimator in (2.8),
χ
b(t) (vi , 0) =

X
h∈L(vi )

|N (h)| (t)
χ
b (vi , 0) =
|N (vi )|

X
h∈L(vi )

We conclude by (3.3) and (3.4) that
χ
b(t) (vi , 0) − χn (vi , 0) =

X
h∈L(vi )


|N (h)| (t)
χ
b (h, 0) − χn (h, 0) .
|N (vi )|

To obtain a concise representation of the asymptotic normal law for the isotropic extremogram,
we define row vectors (|N (h)|/|N (vi )| : h ∈ L(vi )) for i = 1, . . . , p. Set
Pp
L := i=1 `i and define the p × L−matrix



|N (h)|
:
h
∈
L(v
)
0
0
0
1

 |N (v1 )|




|N (h)|
0
:
h
∈
L(v
)
0
0


2
|N (v2 )|
.

N := 

..
..
.
.


. 
0
.
.


|N (h)|
0
0
0
|N (vp )| : h ∈ L(vp )
(3.5)
Then we find

=
d

→

|
n
χ
b(t) (vi , 0) − χn (vi , 0) i=1,...,p
mn
|
n
(i)
(i)
(i)
(i)
b(t) (h`i , 0) − χn (h`i , 0) i=1,...,p
N χ
b(t) (h1 , 0) − χn (h1 , 0), . . . , χ
mn
(space)

N (0, N Π1

|

N ),

n → ∞,

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

12

S. Buhl et al.

such that
(iso)

Π1

(space)

:= N Π1

|

N .

(3.6)

Corollary 3.3. Under the conditions of Theorem 3.2 the averaged spatial extremogram
in (2.10) satisfies
T

n  1 X (tk )
d
(iso)
χ
b (v, 0) − χn (v, 0)
→ N (0, Π2 ),
mn T
v∈V

n → ∞,

k=1

(iso)

with covariance matrix Π2

specified in (3.10) below.

Proof. For the first part of the proof, we neglect spatial isotropy. This part is similar
to the proof of Theorem 4.2 in Buhl and Klüppelberg [6] and Corollary 3.4 of Davis and
Mikosch [9]. We use the notation of the proof of Theorem 3.2. Enumerate the set of spatial
(i)
(i)
lag vectors inherent in the estimation of the extremogram as {h1 , . . . , h`i : i = 1, . . . , p}
and let γ ≥ max{v1 , . . . , vp }. Define the vector process
|

{Y (s) : s ∈ R2 } = {(η(s + h, tk ) : h ∈ B(0, γ))k=1,...,T : s ∈ R2 }.
Let A = B = (1, ∞). Consider i = 1, . . . , p, j = 1, . . . , `i , and k = 1, . . . , T . Define sets
(i)
Dj,k by
(i)

(i)

{Y (s) ∈ Dj,k } = {η(s, tk ) ∈ A, η(s0 , tk ) ∈ B : s − s0 = hj },
and the sets Dk by
{Y (s) ∈ Dk } = {η(s, tk ) ∈ A}.
2

For h ∈ R let BT (h, γ) := B(h, γ) × {t1 , . . . , tT }. For µBT (0,γ) -continuous Borel sets C
and D in R

T |B(0,γ)|

\{0}, regular variation yields the existence of the limit measures
µBT (0,γ) (C) := lim m2n P
n→∞

τBT (0,γ)×BT (h,γ) (C × D) := lim m2n P

 Y (0)
m2n
 Y (0)

n→∞

m2n

∈C



∈ C,


Y (h)
∈
D
.
m2n

By time stationarity we have µBT (0,γ) (Dk ) = µ(A),
(i)
(i)
bm (D(i) , Dk ) := µ
χ
b(tk ) (hj , 0) ∼ R
bBT (0,γ),mn (Dj,k )/b
µBT (0,γ),mn (Dk ),
n
j,k

n → ∞, (3.7)

where the µ
bBT (0,γ),mn (·) are empirical estimators of µBT (0,γ) (·) defined as
µ
bBT (0,γ),mn (·) :=

 m 2 X
n
n

s∈Sn

1{ Y (s)
.
2 ∈·}

(3.8)

mn

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

13

Likewise we have for the pre-asymptotic quantities
(i)

(i)

(i)

χn (hj , 0) = Rmn (Dj,k , Dk ) :=

P(Y (0)/m2n ∈ Dj,k )
P(Y (0)/m2n ∈ Dk )

(i)

=:

µBT (0,γ),mn (Dj,k )
µBT (0,γ),mn (Dk )

,

(3.9)

which are independent of time tk by stationarity. For notational ease we abbreviate in
the following
µBT (0,γ) (·) = µγ (·),

µBT (0,γ),mn (·) = µγ,mn (·),

and µ
bBT (0,γ),mn (·) = µ
bγ,mn (·)

For each k ∈ {1, . . . , T } we now define the matrices
(k)

F (k) = [F1 , F2 ]
(k)

with F1 ∈ RL×L and F2
F1 = diag(µ(A))

∈ RL given by
(k)

and F2

(1)

(1)

(p)

:= (−µγ (D1,k ), . . . , −µγ (D`1 ,k ), . . . , −µγ (D`p ,k ))> .

(k)

Although F2 is constant over k ∈ {1, . . . , T } by time stationarity, we keep the index to
b − χn
clarify the notation. Define the T L × T (L + 1)-matrix F and the column vector χ
with T L components as


(1)
(1)
χ
b(t1 ) (h1 , 0) − χn (h1 , 0)


..


.



 (1)


(1)
F
0
0
0
(t1 ) (1)


χ
b
(h
,
0)
−
χ
(h
,
0)
n
`1
`1
(2)


 0

F
0
0 
..



b
and
χ
−
χ
:=
F :=  .


.
..
.
n
..



 ..
.
0
.
(p)
(t1 ) (p)
χ
b (h`p , 0) − χn (h`p , 0) 
(T )


0
0
0 F


..


.


(p)
(p)
χ
b(tT ) (h`p , 0) − χn (h`p , 0)
b m − Rm ) with the quantities from (3.7) and the corresponding
Define the vector (R
n
n
pre-asymptotic quantities from (3.9) exactly in the same way. Furthermore, define for
k = 1, . . . , T the vectors in RL+1
µ(k)
γ,mn =
|
(1)
(1)
(p)
(p)
µγ,mn (D1,k ), . . . , µγ,mn (D`1 ,k ), . . . . . . , µγ,mn (D1,k ), . . . , µγ,mn (D`p ,k ), µγ,mn (Dk ) ,
b γ,mn
which we stack one on top of the other giving a vector µγ,mn in RT (L+1) , and µ
analogously. Then we obtain
b m − Rm ) = 1 + op (1) F (µ
b − χn = (1 + o(1))(R
b γ,mn − µγ,mn ),
χ
n
n
µ(A)2
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

14

S. Buhl et al.

where the last step follows as in the proof of Theorem 4.2 of [6] and involves Slutzky’s
theorem. Using ideas of the proof of their Lemma 5.1, we observe that as n → ∞,
h
i
Cov µ
bBT (0,γ),mn (C), µ
bBT (0,γ),mn (D)
∼

 m 2 
n

n

µBT (0,γ) (C ∩ D) +

X
06=h∈Z


 m 2
n
τBT (0,γ)×BT (h,γ) (C × D) =:
cC,D .
n
2

With Σ ∈ RT (L+1)×T (L+1) defined as

cD(1) ,D(1) · · · cD(1) ,D1
1,1
 1,1. 1,1 .
..
.
.
Σ=
.
.
.

cDT ,D(1) · · · cDT ,D1

···

cD(1) ,D(p)
1,1
1,T
..
.
cDT ,D(p)

..

.
···

1,1

···
..

.
···

1,T


cD(1) ,DT
1,1

..
,
.

cDT ,DT

we thus conclude that


(1)
(1)
χ
b(t1 ) (h1 , 0) − χn (h1 , 0)
 d
n 
..

 → N (0, µ(A)−4 F Σ(F )> ).
.

mn 
(p)
(p)
χ
b(tT ) (h`p , 0) − χn (h`p , 0)
To obtain the asymptotic covariance matrix in the spatially isotropic case, we proceed
as in the proof of Theorem 3.2. We define the T p × T L-matrix


N 0 0 0
0 N 0 0


N :=  .

.. . .
 ..
. 0
.
0

0

0

N

with N given in equation (3.5). Then we have



(1)
(1)
χ
b(t1 ) (h1 , 0) − χn (h1 , 0)
χ
b(t1 ) (v1 , 0) − χn (v1 , 0)


n 
n

..
..

N

=
.
.


mn
mn
(p)
(p)
χ
b(tT ) (vp , 0) − χn (vp , 0)
χ
b(tT ) (h`p , 0) − χn (h`p , 0)


d

→

N (0, µ(A)−4 N F Σ(N F )> ),

n → ∞,

and we conclude that for the averaged spatial extremogram the statement holds with
1
0

=µ(A)−4 T −2 


0···0
1···0

0

0···1



(iso)

Π2

1
0

0···0
1···0

···
···

1
0


0···0
1 · · · 0

>
 N F Σ(N F )


0···1

···

0

0···1

..

.
0

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes
1
0




0···0
1···0

0

0···1



1
0

0···0
1···0

···
···

1
0

..

.
0

0···1

···

0

>
0···0
1 · · · 0

 .


15

(3.10)

0···1

Remark 3.1. In the central limit theorem the pre-asymptotic extremogram (3.1) can
be replaced by the theoretical one (2.6), provided that
n
(χn (v, 0) − χ(v, 0)) → 0, n → ∞,
(3.11)
mn
is satisfied for all spatial lags v ∈ V. For the Brown-Resnick process (2.1) we obtain from
Lemma 3.1,
n
(χn (v, 0) − χ(v, 0))
mn


P (η(s, t) > m2n , η(s + h, t) > m2n )
n
− χ(v, 0)
=
mn
P(η(s, t) > m2n )


n
∼
χ(v, 0) − 2 χ(v, 0) − 1
2m3n


1
→ 0 if and only if β1 > 1/3;
= n1−3β1 χ(v, 0) − 2 χ(v, 0) − 1
2
cf. Theorem 4.4 of Buhl and Klüppelberg [6]. Thus we have to distinguish two cases:
(I) For β1 ≤ 1/3 we cannot replace the pre-asymptotic extremogram by the theoretical
version, but can resort to a bias correction, which is decribed in (3.14) below.
(II) For 1/3 < β1 < 1/2 we obtain indeed

d
(iso)
n1−β1 χ
b(t) (v, 0) − χ(v, 0) v∈V → N (0, Π1 ), n → ∞.
(3.12)

We now turn to the bias correction needed in case (I). By Lemma 3.1 the preasymptotic extremogram has representation
h 1

i
χn (v, 0) = χ(v, 0) +
χ(v,
0)
−
2
χ(v,
0)
−
1
(1 + o(1))
2m2n
1
= χ(v, 0) +
ν(v, 0)(1 + o(1)), n → ∞,
(3.13)
2m2n


where ν(v, 0) := χ(v, 0) − 2 χ(v, 0) − 1 .
Consequently, we propose for fixed t ∈ {t1 , . . . , tT } and all v ∈ V the bias corrected
empirical spatial extremogram
χ
b(t) (v, 0) −

 (t)

1
1 (t)
χ
b(t) (v, 0) − 2 χ
b (v, 0) − 1 =: χ
b(t) (v, 0) −
νb (v, 0),
2
2mn
2m2n

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

16

S. Buhl et al.

and set

1 (t)
χ
νb (v, 0)
b(t) (v, 0) −
(t)
2m2n
χ
e (v, 0) :=
 (t)
χ
b (v, 0)

if mn = nβ1 with β1 ∈ ( 51 , 13 ],

(3.14)

if mn = nβ1 with β1 ∈ ( 13 , 12 ).

Theorem 3.4 below shows asymptotic normality of the bias corrected extremogram centred by the true one and, in particular, why β1 has to be larger than 1/5.
Theorem 3.4.  For a fixed time point t ∈ {t1 , . . . , tT } consider the spatial BrownResnick process η(s, t), s ∈ R2 defined in (2.1) with dependence function given in (2.2).
Set mn = nβ1 for β1 ∈ 51 , 31 . Then the bias corrected empirical spatial extremogram
(3.14) satisfies

n
d
(iso)
χ
e(t) (v, 0) − χ(v, 0) v∈V → N (0, Π1 ),
mn

n → ∞,

(3.15)

(iso)

where Π1
is the covariance matrix as given in equation (3.6). Furthermore, the corresponding bias corrected version of (2.10) satisfies
T

n  1 X (tk )
d
(iso)
χ
e (v, 0) − χ(v, 0)
→ N (0, Π2 ),
mn T
v∈V

n → ∞,

k=1

(iso)

with covariance matrix Π2

specified in (3.10).

Proof. For simplicity we suppress the time point t in the notation. By (3.13) and (3.14)
we have as n → ∞,
n
n
n
(b
ν (v, 0) − ν(v, 0)).
(e
χ(v, 0) − χ(v, 0)) ∼
(b
χ(v, 0) − χn (v, 0)) −
mn
mn
2m3n
P

ν (v, 0) − ν(v, 0)) → 0. Setting
By Theorem 3.2 it suffices
to show that (n/(2m3n ))(b

νn (v, 0) := χn (v, 0) − 2 χn (v, 0) − 1 we have
n
n
n
(b
ν (v, 0) − ν(v, 0)) =
(b
ν (v, 0) − νn (v, 0)) +
(νn (v, 0) − ν(v, 0)) =: A1 + A2 .
3
3
2mn
2mn
2m3n
We calculate


n
νb(v, 0) − νn (v, 0)
mn (2χ(v, 0) − 3)


n
=
χ
b2 (v, 0) − 3b
χ(v, 0) − (χ2n (v, 0) − 3χn (v, 0))
mn (2χ(v, 0) − 3)


n
=
(b
χ(v, 0) − χn (v, 0))(b
χ(v, 0) + χn (v, 0)) − 3(b
χ(v, 0) − χn (v, 0))
mn (2χ(v, 0) − 3)
χ
n 
b(v, 0) + χn (v, 0) − 3
=
χ
b(v, 0) − χn (v, 0)
.
mn
2χ(v, 0) − 3
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

17

The first term converges by Theorem 3.2 weakly to a normal distribution, and the second
P
P
term, together with the fact that χ
b(v, 0) → χ(v, 0) and χn (v, 0) → χ(v, 0), converges to
P
1 in probability. Hence, it follows from Slutzky’s theorem that A1 → 0. Now we turn to
A2 and calculate
νn (v, 0) = χ2n (v, 0) − 3χn (v, 0) + 2

2


1
1
∼ χ(v, 0) +
ν(v, 0) − 3 χ(v, 0) +
ν(v, 0) + 2
2
2
2mn
2mn
1
1
3
= χ2 (v, 0) − 3χ(v, 0) + 2 + 2 χ(v, 0)ν(v, 0) +
ν(v, 0)2 −
ν(v, 0)
4
mn
4mn
2m2n


1
1
3
= χ(v, 0) − 2 χ(v, 0) − 1 + 2 χ(v, 0)ν(v, 0) +
ν(v, 0)2 −
ν(v, 0)
mn
4m4n
2m2n
ν(v, 0) 
1
3
= ν(v, 0) +
χ(v,
0)
+
ν(v,
0)
−
,
m2n
4m2n
2
where we have used (3.13). Therefore, A2 converges to 0, if n/m5n → 0 as n → ∞. With
mn = nβ1 it follows that β1 > 51 . Finally, the last statement follows from Corollary 3.3.
Remark 3.2. Note that in (3.12) the rate of convergence is of the order na for a ∈
(1/2, 2/3). On the other hand, after bias correction in (3.15) we obtain convergence of
the order na for a ∈ [2/3, 4/5); i.e. a better rate.

Example 3.5. We generate 100 realizations of the Brown-Resnick process in (2.1) using
the R-package RandomFields [29] and the exact method via extremal functions proposed
in Dombry et al. [15], Section 2. We then compare the empirical estimates of the spatial
extremogram χ
b(v, 0) in (2.8) and the bias corrected
e√
(v, 0) in
√ ones
√ χ
√ (3.14)
√ with
√the true
theoretical extremogram χ(v, 0) for lags v ∈ {1, 2, 2, 5, 8, 3, 10, 13, 4, 17}. We
choose the parameters θ1 = 0.4 and α1 = 1.5. The grid size and the number of time
points are given by n = 70 and T = 10. The results are summarized in Figure 3.1. We
see that the bias corrected extremogram is closer to the true one.


3.2. Asymptotic properties of spatial parameter estimates
In this section we prove asymptotic normality of the WLSE (θb1 , α
b1 ) of Section 2. We use
the following notation:


1
yv := 2 log Φ−1 1 − χ
e(v, 0)
and xv := log(v),
2

v ∈ V,

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

18

S. Buhl et al.

Figure 3.1: Empirical spatial extremogram (left) and its bias corrected version (right) for one
hundred simulated max-stable random fields in (2.1) with δ(v, 0) = 2 · 0.4v 1.5 . The dotted line
represents the theoretical spatial extremogram and the solid line is the mean over all estimates.

PT
e(tk ) (v, 0) as in (2.10), possibly after a bias correction, which
with χ
e(v, 0) = T1 k=1 χ
depends on the two cases described in Remark 3.1. Then (2.12) reads as
 
X
2
θb1
= arg min
wv yv − log(θ1 ) + α1 xv
(3.16)
α
b1
θ1 ,α1 >0
α1 ∈(0,2]

v∈V

and we are in the setting of weighted linear regression. To show asymptotic normality of
the WLSE as in (3.16), we define the design matrix X and weight matrix W as
|

X = [1, (xv )v∈V ] ∈ Rp×2

and

|

W = diag{wv : v ∈ V} ∈ Rp×p ,
|

respectively, where 1 = (1, . . . , 1) ∈ Rp . Let ψ 1 = (log(θ1 ), α1 ) be the parameter vector
with parameter space Ψ = R × (0, 2]. Then the WLSE; i.e., the solution to (3.16) is given
by


|
|
|
log(θb1 )
b
= (X W X)−1 X W (yv )v∈V .
ψ 1 :=
α
b1
b may produce estimates of α1 outside its parameter space
Without any constraints ψ
1
(0, 2]. In such cases we set the parameter estimate equal to 2, and we denote the resulting
c |
bc = (log(θbc ), α
estimate by ψ
1
1 b1 ) .
c

|

c
b = (log(θbc ), α
Theorem 3.6. Let ψ
resulting from the constrained
1
1 b1 ) denote the WLSE
|
∗
minimization problem (3.16) and ψ 1 = (log(θ1∗ ), α1∗ ) ∈ Ψ the true parameter vector. Set
mn = nβ1 for β1 ∈ (1/5, 1/2). Then as n → ∞,
(

Z1
if α1∗ < 2,
n  bc
d
∗
ψ1 − ψ1 →
(3.17)
mn
Z2
if α1∗ = 2,

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes
(iso)

where Z 1 ∼ N (0, Π3
Z
P (Z 2 ∈ B) =

), and the distribution of Z 2 is given by

B∩{(b1 ,b2 )∈R2 :b2 <0}
∞Z

Z
+

{b1 ∈R:(b1 ,0)∈B}

0

19

ϕ0,Π(iso) (z1 , z2 )dz1 dz2
3


1
ϕ0,Π(iso) z1 − P
3

v∈V

X
wv


(wv xv ) z2 , z2 dz1 dz2

v∈V

(3.18)
for every Borel set B in R2 , and ϕ0,Σ denotes the bivariate normal density with mean
vector 0 and covariance matrix Σ. In particular, the joint distribution function of Z 2 is
|
given for (p1 , p2 ) ∈ R2 by


|

P Z 2 ≤ (p1 , p2 )



min{0,p
Z 2 } Zp1

ϕ0,Π(iso) (z1 , z2 )dz1 dz2

=

3

+ 1{p2 ≥0}

−∞
Z∞ Zp1

−∞


1
ϕ0,Π(iso) z1 − P
3

v∈V

0 −∞

X
wv


(wv xv ) z2 , z2 dz1 dz2 .

(3.19)

v∈V

The covariance matrix of Z 1 has representation
(iso)

Π3
(iso)

where Π2

(iso)

= Q(w)
x GΠ2

|

GQ(w)
,
x

(3.20)

is the covariance matrix given in (3.10),
|

|

Q(w)
= (X W X)−1 X W
and
x

s
n
2π
1 ∗ α∗1 o
θ
v
:
v
∈
V
.
G = diag
exp
∗
2 1
θ1∗ v α1

(3.21)
(3.22)

Proof. For v ∈ V we have yv = g(e
χ(v, 0)) with g(x) = 2 log(Φ−1 (1 − x/2)). The derivative of g is given by

x
x −1
g 0 (x) = − Φ−1 (1 − ) ϕ(Φ−1 (1 − ))
,
2
2
where ϕ is the univariate standard normal density. Thus,
s
q
q
−1
n1
o

2π
∗
∗
∗ α∗
1
θ1∗ v α1 ϕ θ1∗ v α1
=−
exp
θ
v
.
g 0 (χ(v, 0)) = −
∗
1
2
θ1∗ v α1
Using the multivariate delta method together with Theorems 3.2 and 3.4 it follows that

n
d
(iso)
yv − g(χ(v, 0)) v∈V → N (0, GΠ2 G),
mn

n → ∞,

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

20

S. Buhl et al.

where G is defined in (3.22). Since
X
2 X
2
min
wv g(χ(v, 0)) − log(θ1 ) + α1 xv
=
wv g(χ(v, 0)) − log(θ1∗ ) + α1∗ xv ,
θ1 ,α1

v∈V

α1 ∈(0,2] v∈V

we find the well-known property of unbiasedness of the WLSE,
X
2
|
Q(w)
wv g(χ(v, 0)) − log(θ1 ) + α1 xv
= ψ ∗1 .
x (g(χ(v, 0)))v∈V = arg min
θ1 ,α1

α1 ∈(0,2]

v∈V

It follows that, as n → ∞,




n (w)
n b
d
(iso)
(w) |
ψ 1 − ψ ∗1 =
Qx yv − g(χ(v, 0)) v∈V → N 0, Q(w)
GΠ
GQ
.
x
x
2
mn
mn
We now turn to the constraints on α1 . Since the objective function is quadratic, if the
unconstrained estimate exceeds two, the constraint α1 ∈ (0, 2] results in an estimate
α
b1c = 2. We consider separately the cases α1∗ < 2 and α1∗ = 2; i.e., the true parameter
lies either in the interior or on the boundary of the parameter space. The constrained
bc can be written as
estimator ψ
1
|
bc = ψ
b 1{b
b
ψ
α1 ≤2} + (θ1 , 2) 1{b
α1 >2} .
1
1

We calculate the asymptotic probabilities for the events {b
α1 ≤ 2} and {b
α1 > 2},
 n

n
P(b
α1 ≤ 2) = P
(b
α1 − α1∗ ) ≤
(2 − α1∗ ) .
mn
mn
Since for α1∗ < 2 as n → ∞


 d
|
n
(iso)
α
b1 − α1∗ → N 0, (0, 1)Π3 (0, 1)
mn

n
(2 − α1∗ ) → ∞,
mn

and

it follows that
P(b
α1 ≤ 2) → 1
Therefore, for

α1∗

and P(b
α1 > 2) → 0,

n → ∞.

(3.23)

< 2,
 d
n bc
(iso)
ψ 1 − ψ ∗1 → N (0, Π3 ),
mn

n → ∞.

We now consider the case α1∗ = 2 and α
b1 > 2 (the unconstrained estimate exceeds 2). In
this case (3.16) leads to the constrained optimization problem
|

|

|

min{[W 1/2 ((yv )v∈V − Xψ 1 )] [W 1/2 ((yv )v∈V − Xψ 1 ]},
ψ1

s.t.

(0, 1)ψ 1 = 2.

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

21

bc − ψ ∗ , the vector ψ
b − ψ ∗ is projected onto the line
To obtain asymptotic results for ψ
1
1
1
1
Λ = {ψ ∈ R2 , (0, 1)ψ = 0}, i.e., denoting by I2 the 2 × 2-identity matrix, the projection
|
|
matrix with respect to the induced norm ψ 7→ (ψ X W Xψ)1/2 is given by (cf. Andrews
[1], page 1365)
|

|

|

|

PΛ = I2 − (X W X)−1 (0, 1) ((0, 1)(X W X)−1 (0, 1) )−1 (0, 1).
P
P
For simplicity we use the abbreviation pwx = v∈V wv xv / v∈V wv . We calculate
∗
bc − ψ ∗ )1{b
b
(ψ
α1 >2} = PΛ (ψ 1 − ψ 1 )1{b
α1 >2}
1
1



|
|
|
| −1
−1
−1
b − ψ ∗ )1{b
= (ψ
−
(X
W
X)
(0,
1)
(0,
1)(X
W
X)
(0,
1)
(b
α1 − 2)1{b
α
>2}
α1 >2}
1
1
1


pwx
b − ψ ∗ )1{b
(b
α1 − 2)1{b
= (ψ
α1 >2} +
α1 >2} .
1
1
−1
For the joint constrained estimator ψ c1 we obtain
∗
bc − ψ ∗ = (ψ
bc − ψ ∗ )1{b
bc
ψ
α1 ≤2} + (ψ 1 − ψ 1 )1{b
α1 >2}
1
1
1
1



pwx
∗
∗
b
b
(b
α1 − 2)1{b
= (ψ 1 − ψ 1 )1{b
α1 >2} +
α1 ≤2} + (ψ 1 − ψ 1 )1{b
α1 >2}
−1


b − ψ ∗ ) + pwx (b
= (ψ
α1 − 2)1{b
α1 >2} .
1
1
−1
This implies
n
n bc
(ψ − ψ ∗1 ) =
mn 1
mn



(log(θb1 ) − log(θ1∗ )) + pwx (b
α1 − 2)1{b
α1 >2}
.
(b
α1 − 2) − (b
α1 − 2)1{b
α1 >2}
|

Let f (x1 , x2 ) = (x1 + pwx x2 1{x2 >0} , x2 − x2 1{x2 >0} ) and observe that f (c(x1 , x2 )) =
cf (x1 , x2 ) for c ≥ 0. For the asymptotic distribution we calculate, denoting by f −1 the
inverse image of f ,
 n

bc − ψ ∗ ) ∈ B
P
(ψ
1
mn 1
 n




b − ψ ∗ ) ∈ B = P f n (ψ
b − ψ∗ ) ∈ B
=P
f (ψ
1
1
1
m
mn 1
 nn

b − ψ ∗ ) ∈ f −1 (B ∩ {(b1 , b2 ) ∈ R2 : b2 < 0}) ∪ f −1 (B ∩ {(b1 , 0) : b1 ∈ R})
=P
(ψ
1
1
m
 nn
b − ψ ∗ ) ∈ [B ∩ {(b1 , b2 ) ∈ R2 : b2 < 0}]
=P
(ψ
1
mn 1

∪ [{(b1 − pwx b2 , b2 ), b2 ≥ 0, (b1 , 0) ∈ B}]
Z
→
ϕ0,Π(iso) (z1 , z2 )dz1 dz2
B∩{(b1 ,b2 )∈R2 ,b2 <0}

3

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

22

S. Buhl et al.
Z

∞

Z

+
0

{b1 ∈R,(b1 ,0)∈B}

ϕ0,Π(iso) (z1 − pwx z2 , z2 )dz1 dz2 ,

n → ∞.

3

Plugging in B = (−∞, p1 ] × (−∞, p2 ] and using the Fubini-Tonelli theorem yields (3.19).

Remark 3.3. The derivation of the asymptotic properties for the constrained estimate
is in fact a special case of Corollary 1 in Andrews [1], who shows asymptotic properties
of parameter estimates in a very general setting, when the true parameter is on the
boundary of the parameter space. The asymptotic distribution of the estimates in the
case α1∗ = 2 results from the fact that approximately half of the estimates lie above the
true value and are therefore equal to two, which is reflected by the second term in the
asymptotic distribution of the estimates.


4. Subsampling for confidence regions
Let {η(s, t) : s ∈ R2 , t ∈ [0, ∞)} be the isotropic Brown-Resnick process as in (2.1) with
dependence function δ given in (2.2); i.e.,
δ(v, u) = 2θ1 v α1 + 2θ2 uα2 ,

v, u ≥ 0,

for θ1 , θ2 > 0 and α1 , α2 ∈ (0, 2]. We assume to observe the process on a regular grid Sn =
{(i1 , i2 ) : i1 , i2 ∈ {1, . . . , n}} and at time points {t1 , . . . , tT } = {1, . . . , T } as specified in
Condition 2.4. The estimation method based on the empirical extremogram described in
c |
b c = (log(θbc ), α
bc
bc bc )|
Sections 2 and 3 yields consistent estimators ψ
1
1 b1 )| and ψ 2 = (log(θ2 ), α
2|
?
?
?
?
of the true spatial and temporal parameters ψ 1 = (log(θ1 ), α1 ) and ψ 2 = (log(θ2? ), α2? ) ,
respectively. Furthermore, we have the limit theorems


 c
 c
d
d
?
(1)
(2)
b − ψ? →
b
τn ψ
Z
,
n
→
∞,
and
τ
−
ψ
, T → ∞,
ψ
T
1
2 →Z
1
2
where the bivariate distributions of Z (1) and Z (2) depend on the true parameter q
values
α1? and α2? , respectively. The rates of convergence are given by τn := mnn and τT := mTT ,
where mn and mT are appropriately chosen scaling sequences.
Due to the complicated forms of the covariance matrices of Z (1) and Z (2) (cf. Theorem 3.6 and Theorem 3.19 of [4]) we use resampling methods to construct asymptotic
confidence regions of the parameter vectors ψ ?1 and ψ ?2 . One appealing method is subsampling (see Politis et al. [27], Chapter 5), since it works under weak regularity conditions
asymptotically correct. The central assumption is the existence of weak limit laws, which
is guaranteed by Theorems 3.6 and Theorem 3.19 of [4]. Below we consider the spatial
case, the temporal case is described in Section 3.4.2 of [4].
We have applied subsampling successfully already for confidence bounds of pairwise
likelihood estimates of the space-time Brown-Resnick process in Buhl and Klüppelberg
[5], Section 4. The procedure is as follows: understanding inequalities between vectors
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

23

componentwise, we choose block lengths (1, 1, 1) ≤ b = (bs , bs , bt ) ≤ (n, n, T ) and the
degree of overlap (1, 1, 1) ≤ e = (es , es , et ) ≤ (bs , bs , bt ), where e = (1, 1, 1) corresponds to
maximum overlap and e = b to no overlap. The blocks are indexed by i = (i1 , i2 , i3 ) ∈ N3
T −bt
s
with ij ≤ qs for qs := b n−b
es c + 1 and j = 1, 2 and i3 ≤ qt := b et c + 1. This results in
a total number of q = qs2 qt blocks, which we summarize in the sets

Ei,b,e = (s1 , s2 , t) ∈ Sn × {1, . . . , T } :(ij − 1)es + 1 ≤ sj ≤ (ij − 1)es + bs for j = 1, 2;
(i3 − 1)et + 1 ≤ t ≤ (i3 − 1)et + bt .
We estimate θ1 and α1 based on the observations in each block as described in the
previous sections. This yields different estimates for the spatial parameters, which we
bc .
denote by ψ
1,i
Theorem 4.1 below provides a basis for constructing asymptotically valid confidence
intervals for the true spatial parameters θ1? and α1? . We define τbs as the analogue of τn =
n/mn = n1−β where β ∈ (1/5, 1/2); i.e., τbs := b1−β
(cf. Remark 3.1 and Theorem 3.4).
s
Theorem 4.1.

Assume that the conditions of Theorem 3.6 hold and, as n → ∞,

(i) bs → ∞ such that bs = o(n) (hence, τbs /τn = (bs /n)1−β → 0),
(ii) e does not depend on n.
Define the empirical distribution function Lbs ,s
Lbs ,s (x) :=

qs
qs
qt
1 X X X
1
,
bc
bc
q i =1 i =1 i =1 {τbs kψ1,i −ψ1 k≤x}
1

2

x ∈ R,

3

and the empirical quantile function
cbs ,s (1 − α) := inf {x ∈ R : Lbs ,s (x) ≥ 1 − α} ,

α ∈ (0, 1).

Then


b c − ψ ? k ≤ cb ,s (1 − α) → 1 − α,
P τn kψ
1
1
s

n → ∞.

(4.1)

Proof. We apply Corollary 5.3.4 of Politis et al. [27]. Their main Assumption 5.3.4 is the
b c − ψ ? k. By Theorem 3.6, the continuous
existence of a weak limit distribution of τn kψ
1
1
mapping theorem and the Fubini-Tonelli theorem we have for γ ≥ 0, as n → ∞,
√ 2 2
γ −r
Zγ Z
c
?
b − ψ k ≤ γ) → P(kZ 1 k ≤ γ) = P(Z 1 ∈ B(0, γ)) = 2
P(τn kψ
ϕ0,Π(iso) (r, s)dsdr
1
1
3

−γ

0

if α1? < 2. For α1? = 2 we obtain
b c − ψ ? k ≤ γ) → P(kZ 2 k ≤ γ) = P(Z 2 ∈ B(0, γ))
P(τn kψ
1
1
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

24

S. Buhl et al.
Zγ

Z0
ϕ0,Π(iso) (r, s)dsdr

=
−γ −

√

3

γ 2 −r 2

Zγ Z∞
ϕ0,Π(iso) (r − P

+

3

v∈V

−γ 0

Zγ (
−

X
wv

√

(wv xv )s, s)dsdr

v∈V

Z∞

Z0



ϕ0,Π(iso) r − P

ϕ0,Π(iso) (r, s)ds +

=
−γ

1

3

3

0

γ 2 −r 2

1
v∈V

X
wv

)

(wv xv )s, s ds dr.

v∈V

b c − ψ ? k is conIn particular, the limiting distribution function of the scaled norm τn kψ
1
1
tinuous in γ both for α1? < 2 and α1? = 2. Assumptions (i) and (ii) are also presumed in
Politis et al. [27]. The required condition on the α-mixing coefficients is satisfied, similarly as in the proof of Theorem 3.2, by equation (A.2) and Lemma A.3, and the result
follows.
As a consequence of equation (4.1), for n large enough, an approximate (1 − α)confidence region for the true parameter vector ψ ?1 = (log(θ1? ), α1? ) is given by
b c k ≤ cb ,s (1 − α)/τn }.
{ψ ∈ R × (0, 2] : kψ − ψ
1
s
The one-dimensional approximate (1 − α)-confidence intervals for the parameters θ1? and
α1? can be read off from this as
n c (1 − α) oi
h
n c (1 − α) o
bs ,s
b ,s
, θb1c exp
and
θb1c exp − s
τn
τn
h
cb ,s (1 − α) c cbs ,s (1 − α) i
α
b1c − s
,α
b1 +
∩ (0, 2].
τn
τn

5. Simulation study
We examine the performance of the WLSEs by a simulation study. The estimation of
the spatial parameters relies on a rather large number of spatial observations and the
estimation of the temporal parameters on a rather large number of observed time points.
However, simulation of Brown-Resnick space-time processes based on the exact method
proposed by Dombry et al. [15] can be time consuming, if both a large number of spatial
locations and of time points is taken. For a time-saving method we generate the process
on two different space-time observation areas, one for examining the performance of the
spatial estimates and one for the temporal ones, which we call S (1) ×T (1) and S (2) ×T (2) ,
respectively. The design for the simulation experiment is given in more details as follows:
1. We choose two space-time observation areas
S (1) × T (1) = {(i1 , i2 ) : i1 , i2 ∈ {1, . . . , 70}} × {1, . . . , 10}
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

25

S (2) × T (2) = {(i1 , i2 ) : i1 , i2 ∈ {1, . . . , 5}} × {1, . . . , 300}
√
√ √
√ √
√
and the sets V = {1, 2, 2, 5, 8, 3, 10, 13, 4, 17} and U = {1, . . . , 10}.
2. We simulate the Brown-Resnick space-time process (2.1) based on the exact method
proposed in Dombry et al. [15], using the R-package RandomFields [29]. The dependence function δ is modelled as in (2.2); i.e.,
δ(v, u) = 2θ1 v α1 + 2θ2 uα2 ,

v, u ≥ 0,

with parameters
θ1 = 0.4, α1 = 1.5,

θ2 = 0.2, α2 = 1.

3. The parameters θ1 , α1 , θ2 and α2 are estimated.
• For the estimation of the empirical extremograms (cf. equations 2.8-2.11) we have
to choose high empirical quantiles q. In practice, q is chosen from an interval of
high quantiles for which the empirical extremogram is robust, see the remarks of
Davis et al. [12] after Theorem 2.1. We choose the 90%−empirical quantile for the
estimation of the spatial parameters and the 70%−quantile for the temporal part.
The quantile for the temporal part is lower to ensure reliable estimation of the
extremogram, because the number of time points (300) used for the estimation
of the temporal parameters is much smaller than the number of spatial locations
(70 · 70 = 4900) used for the estimation of the spatial parameters.
• The weights in the constrained weighted linear regression problem (see 2.12 and
2.13) are chosen such that locations and time points which are further apart of each
other have less influence on the estimation. More precisely, we choose
wu = exp{−u2 } for u ∈ U

and

wv = exp{−v 2 } for v ∈ V.

This choice of weights reflects the exponential decay of χ(v, 0) and χ(0, u) defined in
(2.6), which are tail probabilities of the standard normal distribution Φ.
4. Pointwise confidence bounds are computed by subsampling as described in Section 4
for the spatial parameters and in Section 3.4.2 of [4] for the temporal parameters. We
choose block lengths b = (50, 50, 5) and overlap e = (3, 3, 3) for the space-time process
with observation area S (1) × T (1) and b = (5, 5, 200), e = (1, 1, 1) for the process with
observation area S (2) × T (2) .
5. Steps 1 - 5 are repeated 100 times.
Figure 5.1 shows the WLSEs of the spatial parameters θ1 and α1 for each of the 100
realizations of the Brown-Resnick space-time process. The dotted lines above and below
the dots are pointwise confidence intervals based on subsampling. Panel (a) of Table 1
shows the mean, root mean squared error (RMSE) and mean absolute error (MAE) of
both the spatial and the temporal WLSEs based on the 100 simulations. Altogether, we
observe that the estimates are close to the true values. The spatial estimates are slightly
superior to the temporal ones, which is due to the larger number of observations in space
than in time.
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

26

S. Buhl et al.

θ1
α1
θ2
α2

MEAN
0.4033
1.4984
0.2249
0.9563

RMSE
0.0678
0.0521
0.0649
0.0939
(a)

MAE
0.0559
0.0400
0.0526
0.0767

θ1
α1
θ2
α2

MEAN
0.4008
1.4946
0.2188
0.9275

RMSE
0.0668
0.0525
0.0597
0.0976

MAE
0.0552
0.0400
0.0489
0.0799

(b)

Table 1.: Mean, RMSE and MAE of the WLSEs when applied to exact realizations from the
Brown-Resnick process (a) and to realizations with observational noise (b).

Figure 5.1: WLSEs of θ1 (left) and α1 (right) for 100 simulated Brown-Resnick space-time
processes together with pointwise 95%−subsampling confidence intervals (dotted). The middle
solid line is the true value and the middle dotted line represents the mean over all estimates.

In comparison with pairwise likelihood estimation in finite samples (cf. Davis et al.
[11]), a big advantage of the semiparametric method is the substantial reduction of computation time by about a factor 15. Moreover, the semiparametric estimation method is
much more robust when applied to observations that reveal slight deviations from the
model assumptions. To illustrate this point, we repeated the simulation study described
above with data obtained from the original ones by adding to each measurement an
independent N (0, 0.2)-distributed error. The results of the semiparametric estimation
remain practically unaffected; see the summary measures in the panel (b) of Table 1 and
Figure 5.2. In contrast, when applying pairwise likelihood estimation to the same simulated data we observe that the estimates are much more sensitive to small perturbations
than the semiparametric estimation. Whereas for the unperturbed data the estimates are
slightly biased, for the perturbed data the bias increases considerably; their variances,
however, remain nearly unaffected and are (not surprisingly) smaller than the corresponding variances of the semiparametric estimates. Figure 5.3 illustrates the pairwise
likelihood estimates of the spatial parameter θ1 for the simulated Brown-Resnick spaceimsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

27

Figure 5.2: WLSEs of θ1 (left) and α1 (right) for 100 perturbed simulated Brown-Resnick
space-time processes together with pointwise 95%−subsampling confidence intervals (dotted).
The middle solid line is the true value and the middle dotted line represents the mean over all
estimates.

time processes together with 95%-subsampling confidence intervals and for the perturbed
data.

6. Analysis of radar rainfall measurements
Finally, we apply the Brown-Resnick space-time model in (2.1) and the WLSE to radar
rainfall data. The data were collected by the Southwest Florida Water Management
District (SWFWMD)1 . Our objective is to quantify the extremal behaviour of radar
rainfall data in a region in Florida by using spatial and temporal block maxima and
fitting a Brown-Resnick space-time model to the block maxima.
The data base consists of radar values in inches measured on a 120 × 120km region
containing 3600 grid locations. We calculate the spatial and temporal maxima over subregions of size 10 × 10km and over 24 subsequent measurements of the corresponding
hourly accumulated time series in the wet season (June to September) from the years
1999-2004 for further analysis. In this way we obtain 12 × 12 locations during 732 days
containing space-time block maxima of rainfall observations. Taking block maxima allows
one to consider the process to be more consistent with the assumption of a max-stable
process, or at least to lie in the domain of attraction of a max-stable process. Taking
daily maxima, we can furthermore ignore diurnal patterns.
We denote the set of locations by S = {(i1 , i2 ), i1 , i2 ∈ {1, . . . , 12}} and the spacetime observations by {η(s, t), s ∈ S, t ∈ {t1 , . . . , t732 }}. This setup is also considered in
Buhl and Klüppelberg [5], Section 5, and Steinkohl [30], Chapter 7. To make the results
obtained there comparable with the results here, we use the the same preprocessing steps;
for precise description cf. [5], Section 5.1.
1 http://www.swfwmd.state.fl.us/

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

28

S. Buhl et al.

Figure 5.3: Pairwise likelihood estimates of θ1 for 100 simulated Brown-Resnick space-time
processes together with pointwise 95%−subsampling confidence intervals (dotted) in the lefthand plot; the corresponding estimates and pointwise 95%−subsampling confidence intervals
(dotted) for the perturbed simulated processes are presented in the right-hand plot. The middle
solid line is the true value and the middle dotted line represents the mean over all estimates.

The data do not fail the max-stability check described in Section 5.2 of [5], such that
we assume that {η(s, t), s ∈ S, t ∈ {t1 , . . . , t732 }} are realizations of a max-stable spacetime process with standard unit Fréchet margins. Nevertheless, the assumption that the
data are in fact an exact realization from a max-stable process is only approximate. Hence
there is no guarantee that composite likelihood estimation applied to these transformed
data outperforms the semi-parametric estimation introduced in Section 2, cf. the results
obtained in Section 5 when dealing with data with observational noise. Here we use this
data example to illustrate our new semiparametric methodology. As stated in Section 1,
one can then use these estimates to check model adequacy and as initial estimates in
more refined procedures such as composite likelihood.
We fit the Brown-Resnick model (2.1) by estimating the extremal dependence structure
(2.2) as follows:
(1) We estimate the parameters
by WLSE
√ θ1 ,√α1 ,√θ2 and
√ α2 √
√ as described in Section 2
based on the sets V = {1, 2, 2, 5, 8, 3, 10, 13, 4, 17} and U = {1, . . . , 10} for
the linear regression. Permutation tests as described below and visualized in Figure 6.3
indicate that these lags are sufficient to cover the relevant extremal dependence structure. Since the true extremogram χ is unknown, we choose as weights for the different
spatial and temporal lags v ∈ V and u ∈ U the corresponding estimated averaged
PT
Pn2 (si )
extremogram values; i.e., wv = T −1 k=1 χ
e(tk ) (v, 0) and wu = n−2 i=1 χ
e (0, u),
respectively. Since the so defined weights are random, what follows is conditional on
the realizations of these weights.
As the number of spatial points in the analysis is rather small, we cannot choose a
very high empirical quantile q, since this would in turn result in a too small number of
exceedances to get a reliable estimate of the extremogram. Hence, we choose q as the
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

29

empirical 60%−quantile, relying on the fact that the block maxima generate at least
approximately a max-stable process and on the reported robustness of the estimates.
For the temporal estimation, we choose the empirical 90%−quantile for q.
(2) We perform subsampling (see Section 4) to construct 95%-confidence intervals for each
parameter estimate. As subsample block sizes we choose bs = 12 (due to the small
number of spatial locations) and bt = 300. We further choose es = et = 1, which
corresponds to the maximum degree of overlap.
The results are shown in Figures 6.2, 6.3 and Table 2. Figure 6.1 visualizes the daily
rainfall maxima for the two grid locations (1, 1) and (5, 6). The semiparametric estimates
together with subsampling confidence intervals are given in Table 2.
For comparison we present the parameter estimates from the pairwise likelihood estimation (for details for the isotropic Brown-Resnick model see Davis et al. [10] and [30],
Chapter 7), where we obtained θe1 = 0.3485, α
e1 = 0.8858, θe2 = 2.4190 and α
e2 = 0.1973.
From Table 2 we recognize that these estimates are close to the semiparametric estimates
and even lie in most cases in the 95%-subsampling confidence intervals.
Figure 6.2 shows the spatial and temporal mean of empirical temporal (left) and spatial
(right) extremograms as described in (2.10) and (2.11) together with 95% subsampling
confidence intervals. We perform a permutation test to test the presence of extremal independence. To this end we randomly permute the space-time data and calculate empirical
extremograms as before. More precisely, we compute the empirical temporal extremogram
as before and repeat the procedure 1000 times. From the resulting temporal extremogram
sample we determine nonparametric 97.5% and 2.5% empirical quantiles, which gives a
95%−confidence region for temporal extremal independence. The analogue procedure is
performed for the spatial extremogram.
The results are shown in Figure 6.3 together with the extremogram fit based on the
WLSE. The plots indicate that for time lags larger than 3 there is no temporal extremal
dependence, and for spatial lags larger than 4 no spatial extremal dependence.
Estimate
Subsampling-CI
Estimate
Subsampling-CI

θb1
θb2

0.3611
[0.3472,0.3755]
2.3650
[1.9110,2.7381]

α
b1
α
b2

0.9876
[0.9482,1.0267]
0.0818
[0.0000,0.2680]

Table 2.: Semiparametric estimates for the spatial parameters θ1 and α1 and the temporal
parameters θ2 and α2 of the Brown-Resnick process in (2.1) together with 95% subsampling
confidence intervals.

7. Conclusions and Outlook
For the isotropic Brown-Resnick space-time process with flexible dependence structure
we have suggested a new semiparametric estimation method, which works remarkably
well in an extreme value setting. The method results in quite reliable estimates, much
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

30

S. Buhl et al.

Figure 6.1: Daily rainfall maxima over hourly accumulated measurements from 1999-2004 in
inches for two grid locations.

faster than the composite likelihood methods used so far. These estimates can also be
used as initial values for a composite likelihood optimization routine to obtain more
accurate estimates.
Meanwhile, we have generalized the semiparametric method based on extremogram
estimation. The paper Buhl and Klüppelberg [7] is dedicated to the three topics:
1. Generalize the dependence function (2.2) to anisotropic and appropriate mixed
models.
2. Generalize the sampling scheme to allow for a fixed (small) number of spatial
observations and consider limit results for the number of temporal observations to
tend to infinity.
3. Generalize the least squares estimation to estimate spatial and temporal parameters
simultaneously, also in the situation described in topic 2.
Another interesting question concerns the optimal choice of the weight matrix W ,
such that the asymptotic variance of the WLSE is minimal. Some ideas can be found
in the geostatistics literature in the context of least squares estimation of the variogram
parameters; see e.g. Lahiri et al. [25], Section 4. They describe the situation, where the
optimal choice of the weight matrix is given by the inverse
asymptotic
covariance
|
PT of(tthe
matrix of the nonparametric estimates; i.e., of T −1 k=1 χ
e k ) (v, 0) v∈V in the spatial
|
Pn2 (si )
case and of n−2 i=1 χ
e (0, u)−χ(0, u) u∈U in the temporal case. In our case, however,
(iso)
(time)
(given in equations 4.3-4.6 of Buhl and
this involves the matrices Π2
and Π2
Klüppelberg [6]), whose components are infinite sums.

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

31

Figure 6.2: Empirical spatial (left) and temporal (right) extremogram based on spatial and
temporal means for the space-time observations as given in (2.10) and (2.11) together with
95%−subsampling confidence intervals.

Appendix A: α-mixing of the Brown-Resnick
space-time process
In the following we define α-mixing for spatial processes; see e.g. Doukhan [16] or
Bolthausen [2].

Definition A.1. For d ∈ N, consider a strictly stationary random field X(s) : s ∈ Rd
and let d(·, ·) be some metric induced by a norm on Rd . For Λ1 , Λ2 ⊂ Zd set
d(Λ1 , Λ2 ) := inf {d(s1 , s2 ) : s1 ∈ Λ1 , s2 ∈ Λ2 } .
Further, for i = 1, 2 denote by FΛi = σ {X(s), s ∈ Λi } the σ-algebra generated by
{X(s) : s ∈ Λi }.
(i) The α-mixing coefficients are defined for k, l ∈ N ∪ {∞} and r ≥ 0 by
αk,l (r) = sup {|P(A1 ∩ A2 ) − P(A1 )P(A2 )| : Ai ∈ FΛi , |Λ1 | ≤ k, |Λ2 | ≤ l, d(Λ1 , Λ2 ) ≥ r} ,
(A.1)
where |Λi | is the cardinality of the set Λi for i = 1, 2.
(ii) The random field is called α-mixing, if αk,l (r) → 0 as r → ∞ for all k, l ∈ N.
For a strictly stationary max-stable processes Corollary 2.2 of Dombry and EyiMinko [14] shows that the α-mixing coefficients can be related to the extremogram of the
max-stable process. Equations (A.2) and (A.3) follow as in the proofs of Proposition 1
and 2 of Buhl and Klüppelberg [5].

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

32

S. Buhl et al.

Figure 6.3: Permutation test for extremal independence: The gray lines show the 97.5%− and
2.5%−quantiles of the extremogram estimates for 1000 random space-time permutations for the
empirical spatial (left) and the temporal (right) extremogram estimates.


Proposition A.2. For all fixed time points t ∈ N the random field η(s, t), s ∈ Z2
(2.1) is α-mixing with mixing coefficients satisfying
αk,l (r) ≤ 2kl sup χ(s, 0) ≤ 4kle−θ1 r

α1
/2

,

k, l ∈ N, r ≥ 0.

(A.2)

s≥r

For all fixed locations s ∈ R2 the time series {η(s, t) : t ∈ [0, ∞)} in (2.1) is α-mixing
with mixing coefficients satisfying for some constant c > 0
α(r) := α∞,∞ (r) ≤ c

∞
X

ue−θ2 u

α2
/2

r ≥ 0.

,

(A.3)

u=r

We will make frequent use of the following simple result.
Lemma A.3. Let z ∈ N. For (θ, α) ∈ {(θ1 , α1 ), (θ2 , α2 )} and sufficiently large r such
α
that the sequence uz e−θu /2 is monotonously decreasing for u ≥ r, we have
gz (r) =

∞
X

α

uz e−θu

/2

≤ ce−θr

α

/2 z+1

r

r ∈ N.

,

u=r

for some constant c = c(z) > 0.
Proof. An integral bound together with a change of variables yields
z −θr α /2

gz (r) = r e

+

∞
X
u=r+1

z −θuα /2

u e

z −θr α /2

≤r e

Z
+

∞

α

uz e−θu

/2

du

r

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes
 (z+1)/α Z ∞
2
1
t(z+1)/α−1 e−t dt
θ
α θrα /2

= rz e−θr

α

/2

+

≤ rz e−θr

α

/2

+ c1 Γ (d(z + 1)/αe , θrα /2)

= rz e−θr

α

/2

+ c1 (d(z + 1)/αe − 1)! e−θr

d(z+1)/αe−1
α

/2

X
k=0

z −θr α /2

≤r e

+ c2 e

−θr α /2 z+1

≤ ce

r

33

θk rαk
2k k!

−θr α /2 α(d(z+1)/αe−1)

r

,

R∞
Ps−1
where Γ(s, r) = r ts−1 e−t dt = (s − 1)!e−r k=0 rk /k!, s ∈ N, is the incomplete gamma
function and c1 , c > 0 are constants depending on z.

Acknowledgements
The three last authors gratefully acknowledge support by the TUM Institute for Advanced Study. Furthermore, all authors would like to thank Chin Man Bill Mok and
Daniel Straub for their help in finding the data and discussions regarding the results. We
further acknowledge the Southwest Florida Water Management District (SWFWMD) for
providing the data. We also thank Dennis Leber for comparing simulation methods of the
Brown-Resnick process, and we thank Ton Dieker and Marius Hofert for providing and
improving the simulation programme. SB and CS also thank the International Graduate
School of Science and Engineering (IGSSE) of the Technical University of Munich for
support. The research of RD was supported in part by the National Science Foundation
grant DMS-1107031, and ARO MURI grant W11NF-12-1-0385.

References
[1] D. W. K. Andrews. Estimation when a parameter is on a boundary. Econometrica,
67(6):1341–1383, 1999.
[2] E. Bolthausen. On the central limit theorem for stationary mixing random fields.
Annals of Probability, 10(4):1047–1050, 1982.
[3] B.M. Brown and S.I. Resnick. Extreme values of independent stochastic processes.
Journal of Applied Probability, 14(4):732–739, 1977.
[4] S. Buhl. Statistical Modelling and Estimation of Space-Time Extremes. Dissertation,
Technische Universität München, München, 2017.
[5] S. Buhl and C. Klüppelberg. Anisotropic Brown-Resnick space-time processes: estimation and model assessment. Extremes, 19(4):627–660, 2016.
[6] S. Buhl and C. Klüppelberg. Limit theory for the empirical extremogram of random
fields. Submitted. arXiv 1609.04961, 2016.

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

34

S. Buhl et al.

[7] S. Buhl and C. Klüppelberg.
Generalised least squares estimation of regularly varying space-time processes based on flexible observation schemes.
arXiv:1862319[], 2017.
[8] Y. Cho, R.A. Davis, and S. Ghosh. Asymptotic properties of the spatial empirical
extremogram. Scandinavian Journal of Statistics, 43(3):757–773, 2016.
[9] R.A. Davis and T. Mikosch. The extremogram: A correlogram for extreme events.
Bernoulli, 15(4):977–1009, 2009.
[10] R.A. Davis, C. Klüppelberg, and C. Steinkohl. Max-stable processes for extremes
of processes observed in space and time. Journal of the Korean Statistical Society,
42(3):399–414, 2013.
[11] R.A. Davis, C. Klüppelberg, and C. Steinkohl. Statistical inference for max-stable
processes in space and time. Journal of the Royal Statistical Society B, 75(5):791–
819, 2013.
[12] R.A. Davis, T. Mikosch, and Y. Zhao. Measures of serial extremal dependence and
their estimation. Stochastic Processes and Their Applications, 123(7):2575–2602,
2013.
[13] L. de Haan. A spectral representation for max-stable processes. Annals of Probability, 12(4):1194–1204, 1984.
[14] C. Dombry and F. Eyi-Minko. Strong mixing properties of max-infinitely divisible
random fields. Stochastic Processes and their Applications, 122(11):3790–3811, 2012.
[15] C. Dombry, S. Engelke, and M. Oesting. Exact simulation of max-stable processes.
Biometrika, 103:303–317, 2016.
[16] P. Doukhan. Mixing: Properties and Examples. Springer, New York, 1994.
[17] J.H.J. Einmahl, A. Kiriliouk, and J. Segers.
A continuous updating
weighted least squares estimator of tail dependence in high dimension.
arXiv:1601.04826vl[stat.ME], 2016.
[18] S. Engelke, A. Malinowski, Z. Kabluchko, and M. Schlather. Estimation of HüslerReiss distributions and Brown-Resnick processes. Journal of the Royal Statistical
Society B, 77(1):239–265, 2015.
[19] V. Fasen, C. Klüppelberg, and M. Schlather. High-level dependence in time series
models. Extremes, 13(1):1–33, 2010.
[20] M.G. Genton, Y. Ma, and H. Sang. On the likelihood function of Gaussian maxstable processes. Biometrika, 98(2):481–488, 2011.
[21] E. Giné, M. G. Hahn, and P. Vatan. Max-infinitely divisible and max-stable sample
continuous processes. Probab. Theory Rel., 87:139–165, 1990.
[22] H. Hult and F. Lindskog. Regular variation for measures on metric spaces. Publications de l’Institut Mathématique (Beograd), 80:121–140, 2006.
[23] R. Huser and A. Davison. Space-time modelling of extreme events. Journal of the
Royal Statistical Society B, 76(2):439–461, 2014.
[24] Z. Kabluchko, M. Schlather, and L. de Haan. Stationary max-stable fields associated
to negative definite functions. Annals of Probability, 37(5):2042–2065, 2009.
[25] S. N. Lahiri, Y. Lee, and N. Cressie. On asymptotic distribution and asymptotic
efficiency of least squares estimators of spatial variogram parameters. J. Stat. Plan.
Inf., 103(1):65–85, 2002.
imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

Semiparametric estimation for isotropic max-stable space-time processes

35

[26] S.A. Padoan, M. Ribatet, and S.A. Sisson. Likelihood-based inference for max-stable
processes. Journal of the American Statistical Association, 105(489):263–277, 2010.
[27] D.N. Politis, J.P. Romano, and M. Wolf. Subsampling. Springer, New York, 1999.
[28] S.I. Resnick. Heavy-Tail Phenomena, Probabilistic and Statistical Modeling. Springer
Series in Operations Research and Financial Engineering, 2007.
[29] M. Schlather. RandomFields, contributed package on random field simulation for
R. http://cran.r-project.org/web/packages/RandomFields/.
[30] C. Steinkohl. Statistical Modelling of Extremes in Space and Time using Max-Stable
Processes. Dissertation, Technische Universität München, München, 2013.
[31] E. Thibaud and T. Opitz. Efficient inference and simulation for elliptical Pareto
processes. Biometrika, 102(4):855–870, 2015.
[32] J.L. Wadsworth and J.A. Tawn. Efficient inference for spatial extreme value processes associated to log-Gaussian random functions. Biometrika, 101(1):1–15, 2014.

imsart-bj ver. 2014/10/16 file: BDKS_rev2_final_arxiv.tex date: January 11, 2018

