Tail waiting times and the extremes of stochastic processes

arXiv:1512.07848v2 [math.ST] 28 Jul 2016

James E. Johndrow
Duke University, Durham, NC USA
jj@stat.duke.edu

Robert L. Wolpert
Duke University, Durham, NC USA
wolpert@stat.duke.edu

March 12, 2018
Abstract
In applications where extreme dependence at different spatial locations is of interest, data
are almost always time-indexed. When extremes do not occur contemporaneously, existing
methods for inference and modeling in this setting often choose window sizes or introduce
dependence in parameters with the goal of preserving temporal information. We propose an
alternative paradigm for inference on tail dependence in stochastic processes with arbitrary
temporal dependence structure in the extremes, based on the idea that the information on
strength of tail dependence at different locations and the temporal structure in this dependence
are both encoded in waiting times between exceedances of high thresholds at different locations.
We construct a class of space/time-indexed stochastic processes with tail dependence obtained
by endowing the support points in de Haan’s spectral representation of max-stable processes with
velocities and lifetimes. We extend Smith’s model to these max-stable velocity processes and
obtain the distribution of waiting times between extreme events at multiple locations. Motivated
by this result, a new characterization of tail dependence is proposed that is a function of the
distribution of waiting times between threshold exceedances, and an inferential framework is
constructed for quantifying the strength of extremal dependence and assessing uncertainty in
this paradigm. The method is applied to climatological, financial, and electrophysiology data.

Keywords: extreme value; max-stable process; peaks-over-thresholds; tail dependence; time
series; waiting time.

1

Introduction

A max-stable process is defined by de Haan [11] as a stochastic process Y pxq on an index set X with
the property that, for all integers n P N,
D

Y p¨q “

n

1ł
Yi p¨q,
n i“1

(1)

where tYi u are iid copies of Y , where “_” denotes pointwise maximum, and where for two stochastic
D
processes Y, Z the relation Zp¨q “ Y p¨q means that all their finite-dimensional marginal distributions
agree. Slightly different definitions appear elsewhere in the literature (33; 7, §9.3; 30; 5, §8.2; 12,
§9.2). Some authors use the term “simple max-stable process” for those which satisfy (1) (and hence
1

Tail waiting times

2

have Fréchet univariate marginal distributions with shape α “ 1) and extend the class of max-stable
processes to those satisfying
n
”ł
ıL
D
Y p¨q “
Yi p¨q ´ bn p¨q an p¨q
i“1

for suitable sequences of functions an p¨q ą 0, bn p¨q. In the spatial or spatio-temporal setting,
one usually takes X “ Rd for some integer d. This model is quite general in the“ Ž
sense that if
there ‰exist X -indexed sequences an pxq and bn pxq and iid processes wn pxq such that
i“1 wi pxq ´
bn pxq {an pxq converges in distribution then the limit must be a max-stable process [12, Thm 9.2.3]
(after transforming to Fréchet marginals, for a simple max-stable process).
In applications where multivariate or spatial extremes are of interest, typically one has a collection
of observations wpx, tq “ pwpx1 , tq, . . . , wpxn , tqq of a stochastic process tW px, tqu at a collection of
locations x1 , . . . , xn and times t1 , . . . , tp in some study period rT0 , T0 ` T s. These observations could
represent hourly precipitation, maximum daily wind speed, or, if we treat the spatial index set X
as a latent coordinate in an abstract attribute space, essentially any multivariate time series, such
as daily stock prices.
There are two competing paradigms for statistical inference on max-stable processes based on
multivariate time series. The first is the max-over-windows (MoW) approach, which collapses observations wpx, tq of a space/time-indexed process tW px, tqu to observations that are (approximately)
from a max-stable process Y pxq by selecting a time window width ∆t ą 0, setting sj :“ T0 ` j∆t,
and setting yj pxq :“ suptwpx, tq : t P rsj , sj`1 qu [38, 39]. The data yj pxq, j P 0, . . . , tT {∆tu are then
used for estimation and inference.
An alternative is the peaks-over-thresholds (PoT) approach [10, 31], which keeps all those observations wpx, tq that exceed a specified threshold, then treats these observations as realizations
of a max-stable process. Various ways of thresholding have been proposed, including: keeping observations at times t for which maxni“1 wpxi , tq exceeds a pre-specified threshold [28, 6]; fixing a
specific component (say, x1 ) and keeping observations at times t where wpx1 , tq exceeds a threshold
[20, 19, 9, 1]; and, keeping all observations at times t where some vector norm }wp¨, tq} exceeds a
threshold [8, 2].
Extreme events often cluster temporally: they tend to occur nearby in time, but not necessarily
at identical times. Temporal structure can be accounted for in the context of both MoW and PoT
approaches in a variety of ways. Perhaps the most common way is by choosing the window size in
a MoW approach with the goal that dependence across consecutive time windows should be small.
Window lengths are often selected based on observed sizes of extreme observation clusters [40, 23].
A common model for temporally dependent spatial extremes that fits broadly into the MoW context
is the Multivariate Maximum of Moving Maxima (M4) model [34, 21]. However, MoW methods
have the drawback that they tend to discard a larger portion of the data than PoT. Also, artifacts
can arise from the arbitrary location of the time window boundaries. Estimation procedures for
these models also tend to be complex and multi-stage, resulting in potential loss of efficiency.
Under the PoT approach, some temporal structure must be introduced into the model to accommodate dependency for temporally lagged extremes at different locations [see, e.g., 32, 24, 22, 4, 3].
Without such structure, dependency information is lost for observations of a max-stable process
whose threshold exceedances occur at time lags exceeding the sampling rate.
This is illustrated in Figure 1. Extreme values of the variable Y1 are always followed by extreme
values at Y2 exactly two time units later (Y2 also has some additional extreme events that occur
randomly). Two settings in which this aspect of PoT methods is obviously problematic are when
2

J. E. Johndrow and R. L. Wolpert

3

the exceedance at x1 causes a subsequent exceedance at x2 ; and when x1 and x2 are arranged spatially such that extreme events tend to occur first at x1 and then later at x2 , perhaps because the
phenomenon causing the exceedances (e.g., a weather event) moves in that direction. While this can
be addressed by introducing temporal dependence in parameters, this results in more complex modeling, places more demands on usually limited data, and can complicate selection of an appropriate
threshold.
The first main contribution of this work is to propose an alternative paradigm for inference on
tail dependence between spatial locations with arbitrary temporal structure. The underlying idea is
that waiting times between threshold exceedances at different sites contain information on extremal
dependence across both space and time— specifically, that under positive extreme dependence the
time lags between thresholds will be shorter than under independence. To provide a canonical setting
in which to study properties of this approach, we propose a heavy-tailed stochastic process with
spatio-temporal dynamics, so that temporal dependence in extreme events can be accommodated
naturally. To offer intuition for our proposed model, recall that the process
Zpxq :“ supj uj kpx, ξj q is
ş
max-stable [30] if k : X ˆX Ñ R` is a nonnegative kernel satisfying X kpx, ξq mpdξq “ 1 and tuj , ξj u
are points of a Poisson random Borel measure on R` ˆ X with intensity measure proportional to
u´2 du mpdξq for some σ-finite Borel reference measure mpdxq on X . A single realization from Zpxq
for a bivariate Gaussian kernel kpx, ξq is shown in the right panel of Figure 1. Kernels associated

Figure 1: Left panel: example of two time series in which large values of one (Y1 ) are followed by large
values of the other (Y2 ) two time units later. Right panel: a realization of a Gaussian max stable process;
blue indicates small values, white large values.

with large values of uj are clearly discernible. The new stochastic process proposed here endows the
points ξj with random velocities, birth times, lifetimes, and shapes, making the process dynamic
and inducing a temporal dependence structure. As a consequence, time-indexed multivariate data in
this paradigm are viewed as a single realization from what we will call a max-stable velocity process,
rather than multiple (possibly dependent) realizations of a max-stable process. The resulting process
can be visualized as a dynamic version of the sample realization in Figure 1 in which the elliptical
light regions appear at random times with random orientations and eccentricities, move with random
3

Tail waiting times

4

velocities, and then disappear.
The other main contribution of this work is to propose a method for inference on extreme dependence that does not require modeling the process explicitly. Instead, we model waiting times between
thresholds exceedances at pairs of locations, so that inference can be done using only the exceedance
time data. Under positive dependence, waiting times between exceedances are systematically shorter
than they would be under independence.
We show that the waiting times between exceedances for the max-stable velocity process are
distributed as mixtures of exponential distributions and an atom at zero. Further, the waiting
time distributions under independence are also mixtures of exponentials and a point mass, but
with different rate and mixture parameters. By fitting these two models to data and performing
model comparison, we obtain a measure of extreme dependence that naturally captures time lags
and exposes differences in the temporal structure of dependence across space. While we focus on
inference in the setting where the data are thought to originate from a max-stable velocity process,
the use of waiting times to model extreme dependence is general and can be applied in settings
where the max-stable velocity process is considered inadequate by choosing a fully nonparametric
model for the waiting times.

2

Model

In this section we construct the max-stable velocity process and derive some of its properties, with
a focus on the distribution of waiting times between threshold exceedances.

2.1

Max-stable velocity processes

The max-stable velocity process is constructed by extending the spectral characterization of de Haan
[11] and its continuous-path extension by Resnick and Roy [27]:
Theorem 2.1 (after de Haan and Resnick and Roy). Let tpuj , ξj qujě1 be the points of a Poisson
process on R` ˆ Rd with intensity measure proportional to u´2 du dξ. Let tY pxquxPRd be a pathcontinuous max-stable process with unit Fréchet margins. Then there exist nonnegative continuous
functions tkpx, ξq : x, ξ P Rd u such that
ż
kpx, ξqdξ “ 1
@ x P Rd
Rd

for which
"
*
D
tY pxquxPRd “ sup uj kpx, ξj q
jě1

,

(2)

xPRd

D

where “ denotes equality in distribution. Moreover, any process defined by the right side of (2) is
max-stable.
A useful heuristic for de Haan’s spectral characterization is that of weather extremes: the locations tξj u of the support points are taken to be storm centers, the kernel functions kpx, ξj q describe
the shape of the storm, and the marks tuj u quantify storm severity. In this context, the process
realization is the maximum over some period of time of a climatological quantity, such as precipitation or temperature. To create a time-indexed process, we endow the points ξj with birth times,
4

J. E. Johndrow and R. L. Wolpert

5

lifetimes, velocities, and shapes. This approach has the advantage of easily extending the physical
interpretation of the points ξj as storms or, more generally “events.” Now, the storms will move
and have finite lifespans.
Specifically: fix positive numbers β ą 0 and δ ą 0 and a Borel probability measure πpdaq on a
Polish “attribute space” A. Define a σ-finite Borel measure
νpdωq :“ βu´2 du dξ dσ δe´δτ dτ πpdaq
(3)
`
˘
on the space Ω :“ R` ˆ Rd ˆ R ˆ R` ˆ A, and let N pdωq „ Po νpdωq be a Poisson random
measure on Ω with intensity νpdωq— i.e., a measure that assigns independent random variables
d
N pBj q „ PopνpBj qq to disjoint Borel sets
ş Bj Ă Ω. Fix a nonnegative function k : R ˆ R ˆ
d
R ˆ R ˆ R` ˆ A Ñ R` that satisfies Rd kpx, pqx, t; ξ, σ, τ, aq dξ “ 1 for each x, t, σ, τ, a. Using
the climatological heuristic, the support points tωj u “ tpuj , ξj , σj , τj , aj qu of N pdωq can be thought
of as representing storms of magnitudes uj ą 0, initiating at locations ξj P Rd at times σj P R,
with lifetimes τj P R` and attribute vectors aj that may include velocities vj , shapes Λj , or other
features. For any location x P Rd and time t P R, define
Y px, tq :“ suptuj kpx, t; ξj , σj , τj , aj qu
j

Y ˚ px, tq :“ sup Y px, sq,

tPR

(4a)

t ě 0.

(4b)

0ďsďt

We refer to (4a) as a max-stable velocity (MSV) process and (4b) as the corresponding maximal
MSV process.
In the sequel we will take A “ pRd ˆ P d q with elements aj “ pvj , Λj q P A that consist of a
d
velocity vector vj P Rd and a shape matrix Λj P P d , an element
ş of the cone P of positive-definite
dˆd matrices. Let ϕ : Rd Ñ R` be a continuous pdf satisfying Rd ϕpzqdz “ 1 that is non-increasing
in z 1 z and, for Λ P P d , set ϕΛ pzq :“ |Λ|1{2 ϕpΛ1{2 zq (here |Λ| denotes the determinant of Λ P P d ).
We take k to have the specific form
`
˘
kpx, t; ξ, σ, τ, aq “ ϕΛ x ´ ξ ´ vpt ´ σq 1tσďtăσ`τ u ,
(5)
the magnitude at time t and location x P Rd of a storm of unit severity that originated at location
ξ P Rd at a time σ ă t and moved at velocity v P Rd for time pt ´ σq. We write Y „ MSVpβ, δ, π, ϕq
to denote a process of the form (4a) with k defined by (5) and νpdωq as in (3). This includes and
generalizes the Gaussian and Student t kernels proposed as models for the generic max-stable process
by Smith [33].

2.2

Main results

Max-stable velocity processes are in fact space-time indexed forms of the max-stable process, as
shown in Theorem 2.2. All proofs are deferred to the Appendix or found in Supplementary Materials.
Theorem 2.2 (Max-stable velocity processes are max-stable). If Y „ MSVpβ, δ, π, ϕq, then Y is
max-stable jointly in px, tq on the index set Rd ˆ R` .
Theorem 2.2 is shown by verifying the definition in (1) directly— that is, by verifying that the
process is “stable” under finite maxima. It turns out that the maximal MSV process Y ˚ px, tq of
(4b) is also max-stable.
5

Tail waiting times

6

Theorem 2.3 (Maximal MSV process is max-stable). If Y „ MSVpβ, δ, π, ϕq, then Y ˚ px, tq is
max-stable in px, tq on the index set Rd ˆ R` .
Therefore, if data are generated from a MSV process and transformed by taking blockwise maxima, the transformed data are realizations of a max-stable process. This result provides a link
between traditional methods and the alternative we propose.
In addition to the marginal and joint distribution of the max-stable velocity process, we are
also interested in the distribution of waiting times between threshold exceedances at distinct points,
as these will eventually form the basis for inference. For two arbitrary locations x1 , x2 P Rd and
threshold levels y1 , y2 P R` , define the random variables
κi :“ inf tt : Y pxi , tq ą yi u,
tě0

i “ 1, 2

Then κ∆ :“ |κ1 ´ κ2 | is the waiting time between first exceedances of levels y1 at x1 and of y2 at
x2 , in either order. Theorem 2.4 gives the marginal distribution of the max-stable velocity process
Y px, tq and that of the waiting times κi until first exceedance.
Theorem 2.4 (Marginals and waiting times). The max-stable velocity process has unit Fréchet
marginals, with distribution function PrY px, tq ď ys “ e´β{δy for y ą 0. The marginal waiting time
distribution is given for t ě 0 by
"
*˙
ˆ
ż
β
β
δ`
ϕΛ pζq dζ |v| πpdv dΛq
,
(6)
´t
Prκi ą ts “ exp ´
δyi
δyi
Rd ˆP d ˆv K
where v K :“ tζ P Rd : v 1 Λζ “ 0u denotes the orthogonal complement of the span of v in the Λ
metric. This is a mixture of a point-mass at zero with weight 1 ´ exppβ{δyi q and an exponential
distribution with a shape-dependent rate constant that increases monotonically in the mean particle
speed.
Corollary 2.5 (Distribution of Y ˚ ). The distribution of Y ˚ pxi , tq “ PrY ˚ pxi , tq ă yi s is given by
(6).
Thus the waiting time distribution and the distribution of Y ˚ px, tq, unlike the marginal distribution of Y px, tq, depends on velocity and shape. For Gaussian kernels ϕ, (6) has a simple expression:
Corollary 2.6 (Waiting time distribution for Gaussian kernels). For Gaussian kernels ϕpzq “
p2πq´d{2 expp´z 1 z{2q,
ˆ
˙
”`
˘ 12 ı)
β
β !
1
Prκi ą ts “ exp ´
´t
δ ` Eπ v Λv{2π
.
(7)
δyi
δyi
“?
‰
Thus in the Gaussian case, the marginal waiting times converge to zero in distribution as Eπ v 1 Λv Ñ
8. Theorem 2.7 gives the joint distribution of the max-stable velocity process at finite collections
of locations txi u1ďiďn and times tti u1ďiďn .
Theorem 2.7 (Joint CDF). Let Y „ MSVpβ, δ, π, ϕq be a max-stable velocity process and let
txi u1ďiďn Ă Rd and tti u1ďiďn Ă R for some integer n P N. The joint CDF for tY pxi , ti qu1ďiďn is
given by
`
˘
`
˘
F py1 , . . . , yn q :“ P X1ďiďn rY pxi , ti q ď yi s “ exp ´ νpY1ďiďn Bi q
(8a)
6

J. E. Johndrow and R. L. Wolpert

7

for Bi :“ tω : Y pxi , ti q ą yi u, with
ÿ
ÿ
ÿ
`
˘
ν Y1ďiďn Bi “
νpBi q ´
νpBi X Bj q `
νpBi X Bj X Bk q ´ ¨ ¨ ¨ ,
1ďiďn

i‰j

(8b)

i‰j‰k

the alternating sum of terms
`
˘ β
J
ν XiPJ Bi “ e´δpt ´tJ q
δ

"

ż
min

Rd ˆA jPJ

*
ϕΛ pxj ´ vtj ´ zq
dz πpdv dΛq
yj

(8c)

for subsets J Ď t1, ¨ ¨ ¨ , nu. Here tJ :“ maxttj u denotes the maximum and tJ :“ minttj u the
minimum of ttj ujPJ .
While the univariate marginal distributions of Y px, tq do not depend on πpdv dΛq at all, higher
order marginal distributions do. For example,
Corollary 2.8 (Gaussian MSV process joint CDF at two points). For the Gaussian max-stable
´d{2
velocity
expp´z 1 z{2q, the bivariate distribution is given by F py1 , y2 q “
` process with˘ ϕpzq “ p2πq
exp ´ νpB1 Y B2 q with
¯´ 1
1 ¯ β ´δ|t2 ´t1 |
β´
1 ´ e´δ|t2 ´t1 |
`
` e
νpB1 Y B2 q “
δ
y1
y
δ
"
˙2
˙*
ˆ
ˆ
ż
1
1
SΛ pvq logpy1 {y2 q
SΛ pvq logpy2 {y1 q
ˆ
´
` Φ
´
πpdv dΛq
(9)
Φ
y1
2
SΛ pvq
y2
2
SΛ pvq
P d ˆRd
a
where Φp¨q is the standard Gaussian CDF and SΛ pvq :“ ∆1v Λ∆v for ∆v :“ px2 ´ x1 q ´ pt2 ´ t1 qv.
Equation (9) generalizes Smith [33, equation 3.1], and reduces to it if t1 “ t2 and if π is a unit
point mass. It leads to an explicit expression for the likelihood function for the Gaussian max-stable
velocity process with observations at two locations and times. This expression could be used to fit
the max-stable velocity process to data by maximum composite likelihood.
Theorem 2.9 gives an expression for the waiting times between first exceedances under independence, which provides a null distribution in the approach to inference proposed in Section 3.
Theorem 2.9 (Waiting times between exceedances under independence). Suppose that κ1 and κ2
are independent, each with the distribution given in (6) of Theorem 2.4. Then the distribution of
κ∆ is given for t ě 0 by
`
˘
`
˘
Prκ∆ ą ts “ 1 ´ e´a2 e´a1 ´b1 t
` 1 ´ e´a1 e´a2 ´b2 t
b2 e´a1 ´a2 ´b1 t
b1 e´a1 ´a2 ´b2 t
e
`
e
b `b
b1 ` b2
" 1 ´a21 ´a2
*
"
*
b1 e
p1 ´ e´a2 q
1 ´ e´a1
b2 e´a1 ´a2
´b1 t
“
`
e
`
`
e´b2 t
b1 ` b2
ea1
b1 ` b2
ea2

`

where
ai “

β
δyi

bi “

β
δyi

"
ż
δ`

(10)

*
ϕΛ pζq dζ |v| πpdv dΛq .

Rd ˆP d ˆv K

Thus, when κ1 K
K κ2 , κ∆ is distributed as a mixture with three components: a point mass at zero and
two exponential distributions with ”rates b1 andıb2 . For the Gaussian kernel, the integral expression
1
over P d ˆ Rd ˆ v K is given by Eπ pv 1 Λv{2πq 2 .
7

Tail waiting times

8

Theorem 2.10 gives the distribution of κ∆ in the special case of v “ 0. This result is useful in
understanding the effect of velocity on the waiting time distribution in the general case.
Theorem 2.10 (waiting times between exceedances when v ” 0). Suppose πpt0u ˆ P d q “ 1 in the
Gaussian max-stable velocity model, so v ” 0 almost surely. Then the distribution of κ∆ is given
for t ě 0 by
`
˘
`
λ0 ˘
Prκ∆ “ 0s “1 ´ e´λ0 e´λ1 ` e´λ2 ` e´λ` 1 `
λ`
“
‰
λ
`
λ
0
2 ´λ1
Prκ∆ ą ts “ e´pλ0 `λ2 qp1`δtq 1 ´
e
λ`
“
λ0 ` λ1 ´λ2 ‰
` e´pλ0 `λ1 qp1`δtq 1 ´
e
λ`

(11a)

(11b)

a mixture of an atom at zero and two exponential distributions with rates λ1 , λ2 . The quantities
tλj u Ă R` are given by λ` :“ λ0 ` λ1 ` λ2 and, for j “ 0, 1, 2,
„

ż
1
1
β
1 ϕΛ pzq ă y1 u ` 1t ϕΛ pzq ă y2 u ϕΛ pzq dz πpdaq
λ0 “
y2
δ Rd ˆA y1 t ϕΛ pz`∆
y2 ϕΛ pz´∆0 q y1
0q
„

ż
1
1
β
λ1 “
1 ϕ pz`∆0 q y ´ 1 ϕΛ pzq ă y2 u ϕΛ pzq dz πpdaq
y1
δ Rd ˆA y1 t ΛϕΛ pzq ă y12 u y2 t ϕΛ pz´∆
0q
„

ż
β
1
1
λ2 “
1 ϕ pz´∆0 q y ´ 1 ϕΛ pzq ă y1 u ϕΛ pzq dz πpdaq
y2
δ Rd ˆA y2 t ΛϕΛ pzq ă y21 u y1 t ϕΛ pz`∆
0q
„

ż
1
β
1
λ` “
1 ϕΛ pzq ă y1 u ` 1t ϕΛ pzq ą y2 u ϕΛ pzq dz πpdaq.
y2
δ Rd ˆA y1 t ϕΛ pz`∆
y2 ϕΛ pz´∆0 q y1
0q
with ∆0 :“ px2 ´ x1 q as in Theorem 2.7. The rates are available in closed form for the Gaussian
kernel:
*
ż "
1 ´ 1
1 ´ 1
β
y1
SΛ ¯
y2
SΛ ¯
` Φ
πpdaq
λ0 “
Φ
log
´
log
´
δ A y1
SΛ
y2
2
y2
SΛ
y1
2
*
ż "
1 ´ 1
1 ´ 1
β
y2
SΛ ¯
y2
SΛ ¯
´ Φ
πpdaq
λ1 “
Φ
log
`
log
´
δ A y1
SΛ
y1
2
y2
SΛ
y1
2
*
ż "
1 ´ 1
β
y1
SΛ ¯
1 ´ 1
y1
SΛ ¯
λ2 “
Φ
log
`
´ Φ
log
´
πpdaq
δ A y2
SΛ
y2
2
y1
SΛ
y2
2
*
ż "
β
1 ´ 1
y2
SΛ ¯
1 ´ 1
y1
SΛ ¯
λ` “
Φ
log
`
` Φ
log
`
πpdaq
δ A y1
SΛ
y1
2
y2
SΛ
y2
2
and simplify considerably if y1 “ y2 (see Supplementary Materials).
The next result gives a stochastic bound on the distribution of κ∆ in the general case. The result
is given here for the Gaussian kernel; expressions for the general case can be found in the Appendix.
Theorem 2.11 (Stochastic bound for survival function with nonzero velocity). Suppose Y „
MSVpβ, δ, π, ϕq is a Gaussian max-stable velocity process. Then for any two points x1 , x2 and thresholds y1 , y2 , κ∆ satisfies the stochastic ordering
(
Prκ∆ ą ts ď exp ´ ηpAq ,
8

J. E. Johndrow and R. L. Wolpert

9

with
β
ηpAq “ ?
2π

"

˙
ˆ
1
SΛK pvq logpy1 {y2 q
Φ
´
`
Rd ˆP d
y1
2
SΛK pvq
Xt´tď∆12 ďtu
˙*
ˆ
?
1
SΛK pvq logpy1 {y2 q
` Φ ´
´
v 1 ΛvπpdΛ dvq,
K
y2
2
SΛ pvq

ż

pt ´ ∆12 q e´δ∆12

(13)

where SΛK pvq :“ tpx2 ´ x1 q1 Λpx2 ´ x1 qK u1{2 , ∆12 :“ v 1 Λpx2 ´ x1 q{v 1 Λv, and px2 ´ x1 qK :“ px2 ´
x1 q ´ ∆12 v is the projection of px2 ´ x1 q into the orthogonal complement of v in the Λ metric.
Theorem 2.11 shows that for the Gaussian MSV process, κ∆ converges to zero in probability if
„
ˆ b
˙
?
Eπ
v 1 ΛvΦ ´ px2 ´ x1 q1 Λpx2 ´ x1 qK
Ñ 8.
?
This condition is only slightly stronger than the condition Eπ r v 1 Λvs Ñ 8, which is sufficient for
the marginal waiting times to converge to zero in probability based on (7). Informally, the difference
between these conditions is that πpdv dΛq cannot place too much mass on Λ with large eigenvalues,
which corresponds to extremely concentrated kernels. This is intuitive: it is difficult for the same
point to cause an exceedance at two different locations when the kernels are extremely concentrated,
since the volumes of space where exceedances can occur are very small.
It is instructive to compare the result in Theorems 2.10 and 2.11. When v ” 0, all exceedances
occur at birth times or at time zero, so points ω P Ω leading to exceedance when πpt0uˆP d q “ 1 will
also lead to exceedance for any other velocity distribution. Therefore, (11) also provides a stochastic
bound on the survival function with nonzero velocity. Taken together, these results imply that
expected waiting times decrease as mean speed increases, and the distribution of κ∆ is dominated
by a mixture of an atom at zero and exponential distributions. Thus, such mixtures provide a useful
model for κ∆ with arbitrary πpdv dΛq, and the model space of finite mixtures of exponentials with
an atom at zero contains the distribution of κ∆ in the max-stable velocity model.

3

Inference: Tail waiting times

Waiting times between exceedances in the max-stable velocity process are given by simple mixtures
of exponential distributions and a point mass at zero. Thus, we propose to use the waiting times
between exceedances as data to perform inference on tail dependence, bypassing pseudo-likelihood
procedures. In this section we propose a novel quantification of tail dependence based on waiting
times. We then propose a Bayesian approach to inference for this quantity.

3.1

Waiting times as a measure of extremal dependence

Classically, the tail dependence index of a stochastic process Y pxq at points x1 , x2 is defined as the
limit γ :“ limyÑ8 PrY px1 q ą y | Y px2 q ą ys of the conditional survival probability. A conceptually
similar approach is possible for tail waiting times κ∆ . Fix x1 , x2 P Rd and let tµ∆ u be the family of
Borel probability measures on R induced by the collection of random variables κ∆ , indexed by pairs
py1 , y2 q P R2` , so that for any Borel set A and any py1 , y2 q we have µ∆ pA; y1 , y2 q “ Prκ∆ P As. We
K
suppress dependence on py1 , y2 q to simplify notation. Similarly, let µK
∆ be the measure induced by
9

Tail waiting times

10

κ∆ under the assumption κ1 K
K κ2 , the product measure of the marginal distributions of κ1 and κ2 .
For any metric dp¨, ¨q on the space of probability measures, define a tail dependence index γd based
on tail waiting times by
˘
`
K
γd8 :“ lim γd .
(14)
γd :“ d µ∆ , µK
∆ ,
y1 ,y2 Ñ8

In this framework, a value of γd8 “ 0 corresponds to convergence (in the topology induced by d) as
y1 , y2 Ñ 8 of the distribution of κ∆ to its distribution under independence, while γd8 ą 0 indicates
a discrepancy between the measures that persists asymptotically.
Inference on γd requires specification of a metric dp¨, ¨q on probability measures, which we now
construct. Let pH, x¨, ¨yH q be a reproducing kernel Hilbert space on R with reproducing kernel
K : R ˆ R Ñ R and unit ball H1 :“ tf P H : xf, f yH ď 1u. For each θ P R denote by Kθ P H the
function Kpθ, ¨q, that satisfies xKθ , hyH “ hpθq for all h P H. ş a
The set PK :“ t Borel probability measures
µ on R s.t. R Kpθ, θq µpdθq ă 8u can be embedş
H
ded into H by the mapping µ ÞÑ µ :“ R Kθ µpdθq. This induces a pseudo-metric on PK by
ˇ
ˇż
ˇ
ˇ
H
H
ˇ
(15)
dK pµ1 , µ2 q :“ }µ1 ´ µ2 }H “ sup ˇ hpθqpµ1 ´ µ2 qpdθqˇˇ
hPH1

Θ

[17, §2.3]. For characteristic kernels, including strictly positive-definite ones like the Gaussian
kernel we use KG pθ, θ1 q :“ expp´|θ ´ θ1 |2 q, (15) is in fact a metric (i.e., satisfies dK pµ1 , µ2 q “ 0 if
and only if µ1 “ µ2 )— a simple consequence of Bochner’s and Plancherel’s Theorems for strictly
positive-definite kernels [25, §2.3]. It is this metric we use as “dp¨, ¨q” in (14).
This family of metrics has been studied extensively in the literature [16, 35, 36, 37, 17, 25], and
an estimator based on empirical measures is given by Minsker et al. [25, Equation (2.12)]. This
estimator is both computationally tractable and robust for small samples and high dimensions,
motivating our use of it.
Moreover, conditions similar to those required for convergence of the posterior in the Hellinger
metric [see 15] imply convergence in dK , at the same rate [25]. For simplicity, in the sequel we
refer to dK as the “RKHS distance,” and all of the estimates presented utilize the Gaussian kernel
Kpθ, θ1 q “ expp´|θ ´ θ1 |2 q.
One alternative
to the RKHS
ˇ
ˇ metric dK is the total variation distance (TV), given by }µ1 ´
µ2 }TV :“ supA ˇµ1 pAq ´ µ2 pAqˇ, where the supremum is taken over Borel sets A Ă R. For one
dimensional distributions, a point estimate of }µ1 ´ µ2 }TV is provided by the Kolmogorov-Smirnov
statistic. Another alternative is the Wasserstein W1 metric, to which dK bears some relation.

3.2

Estimating waiting times from observed data

Consider heavy-tailed data wpx, tq which feature (possibly after transformation) unit Fréchet marginals.
For each finite set txi u Ă Rd : i P 1, . . . , n of locations, select thresholds tyi u Ă R` corresponding
to high quantiles of the Fréchet distribution. Compute the marginal waiting times tκi u until first
exceedance using Algorithm 1.
The subtraction of ptl`1 ´tl q`ptf ´tf ´1 q from ptf ´tl q in Algorithm 1 is a convention that is easiest
to explain by example. Suppose data are observed at three consecutive time points tj´2 , tj´1 , tj ,
and that wpxi , tj´2 q ą yi , wpxi , tj q ą yi but wpxi , tj´1 q ď yi . We cannot know for  ă minptj´1 ´
tj´2 , tj ´ tj´1 q whether the value of wpxi , tj´1 ˘ q ą yi or not. Therefore, the waiting time
between first exceedances could be arbitrarily close to zero, or as long as tj ´ tj´2 . Subtraction
10

J. E. Johndrow and R. L. Wolpert

11

Data: times t1 , . . . , tm , data wpxi , tj q for j ď m, threshold yi
Result: Waiting times κi until first exceedances of yi at location xi
F “ tj : 1 ă j ď m, wpxi , tj q ą yi , wpxi , tj´1 q ď yi u; F “ sortpF, increasingq;
L “ tj : 1 ă j ď m, wpxi , tj q ď yi , wpxi , tj´1 q ą yi u; L “ sortpL, increasingq;
if min F ď min L then
F “ F2:|F | ;
end
if max L ě max F then
L “ L1:p|L|´1q ;
end
N “ |F |; κ “ vector of length N ;
for n “ 1, . . . , N do
f “ Fn ; l “ Ln ;
κn “ ptf ´ tl q ´ pptl`1 ´ tl q ` ptf ´ tf ´1 qq;
end
return κi “ κ;
Algorithm 1: Computation of marginal waiting times between first exceedances

of ptl`1 ´ tl q ` ptf ´ tf ´1 q results in an observed value of zero. In our approach to inference, we
appropriately adjust for the induced bias by imputing observed waiting times. Algorithm 2 returns
observed waiting times between first exceedances at pairs of sites pxi , xi1 q.

Data: times Ti , Ti1 , data wpxi , tq for t P Ti , wpxi1 , tq for t P Ti1 , thresholds yi , yi1
Result: Waiting times κ∆ between first exceedances of yi at location xi and yi1 at location xi1
Fi “ tt P Ti : wpxi , tq ą yi , wpxi , maxs ts P Ti , s ă tuq ď yi u; Fi “ sortpFi , increasingq;
Fi1 “ tt P Ti1 : wpxi1 , tq ą yi1 , wpxi1 , maxs ts P Ti1 , s ă tuq ď yi1 u;
N “ |Fi |; κ “ vector of length N ;
for n “ 1, . . . , N do
t “ Fi,n ; ∆ “ minti1 PFi1 |t ´ ti1 |;
κn “ ∆;
end
return κ∆ “ κ;
Algorithm 2: Computation of waiting times between first exceedances at two locations

The resulting data are n vectors κi of marginal waiting times until first exceedances and npn´1q{2
vectors of waiting times κ∆ between first exceedances at pairs of locations pxi , xi1 q. Clearly, the larger
yi , the fewer observations of first exceedances exist. The effect of the threshold on sample size and
the necessity of choosing yi large enough that values exceeding this threshold have approximately a
limiting max-stable velocity distribution are the two opposing considerations in choosing a threshold.
Our goal here is mainly to illustrate the potential of the method, so we take the simple approach of
choosing thresholds such that there are at least 100 observations of κi for each i.
11

Tail waiting times

12

3.3

Estimation of γd

K
We take a model-based approach to estimating µ∆ and µK
∆ , guided by Theorems 2.9 and 2.11. If
the data originate from a max-stable velocity process, then the marginal waiting time distribution
is the mixture of an atom at zero and an exponential distribution given by Theorem 2.4, and the
distribution of waiting times between exceedances is stochastically dominated by the expressions
in Theorems 2.10 and 2.11. As such, we model κi and κ∆ using a mixture of several exponential
components and an atom

κ „ η0 δ0 `

K´1
ÿ

ηj Exponentialpλj q.

(16)

j“1

We take a Bayesian approach to inference in this model by specifying the prior distributions
η „ Dirichletp1{K, . . . , 1{Kq, and λj „ Gammap1, 1q, and write ψ “ pη, λq to represent the model
parameters. Computation is via a straightforward Gibbs sampler, described in Supplementary Materials. All of the empirical results in the sequel were obtained using K “ 11 components, and point
estimates always correspond to the ergodic average from sample paths obtained from the Gibbs
sampler.
The motivation for a Bayesian approach is threefold. First, censoring of κ is common in applications because data are sampled at discrete intervals. This is easy to handle in a Bayesian model
by imputing uncensored waiting times. Second, by setting K to be relatively large, the posterior
is consistent for the true number of mixture components, and the extraneous components tend to
empty [29]. Thus, the Bayesian approach provides an easy way to estimate the number of mixture
components and assess model fit. Finally, a Bayesian approach provides estimates of uncertainty in
K
µ∆ and µK
∆ without the use of additional procedures.
Estimating models of the form in (16) by MCMC yields approximate posterior samples tψi u from
πpψi | κi q for i “ 1, . . . , n and tψ∆ u from πpψ∆ | κ∆ q for each pair of sites pi, i1 q P t1, . . . , nu2 . To
obtain approximate samples from the posterior distribution
πpγd | κq “ πpγd | κi , κi1 , κ∆ q,

(17)

for any pair of sites, follow the procedure:
(1) For each sampled value of tψi , ψi1 u, take M samples from the posterior predictive distribution
K
πpκK
∆ | ψi , ψi1 q, e.g., by taking M independent samples from the posterior predictive distributions
K
πpκi | ψi q and πpκi1 | ψi1 q and computing κK
∆ “ |κi ´ κi1 | for each sampled value.
(2) For each sampled value of ψ∆ , take M samples from the posterior predictive distribution πpκ∆ |
ψ∆ q.
K
(3) Compute an estimate of γd “ dpµ∆ , µK
∆ q based on the M samples taken in the previous two
steps.

(4) Repeat (1)–(3) for each retained sample of tψi , ψi1 , ψ∆ u.
The resulting samples are approximately samples from (17). All point estimates γ
bd in the empirical
results were obtained from ergodic averages of γd .
12

J. E. Johndrow and R. L. Wolpert

3.4

13

Posterior Inference

We now suggest an approach for characterizing the strength of tail dependence given the posterior
distribution in (17). For any ψ, let κ∆ pψq be a random variable with distribution given by (16),
and let µ∆ pψq be the measure induced by κ∆ pψq. Consider the random variable d˚ given by d˚ :“
ind
dpµ∆ pψ1 q, µ∆ pψ2 qq, where ψ1 , ψ2 „ πpψ | κ∆ q. The distribution of d˚ is the posterior distribution
of distances between measures induced by the likelihood in (16), and reflects uncertainty about the
parameters ψ after having observed the waiting times κ∆ . If there is strong tail dependence between
sites xi and xi1 , we would expect γ
bd to be large relative to d˚ . Thus, as a basic measure of the
relative magnitude of γd , we estimate
ı
”
pd “ Prγd ą d˚ | κs “ Eπpγd |κqˆπpd˚ |κ∆ q 1td˚ ăγd pµ∆ ,µK∆K qu ;
if pd is near 1, there is strong evidence for tail dependence.

4

Simulation

In this section, a simulation study is constructed to illustrate the method. The simulation is motivated by extreme weather events, where basic scientific knowledge allows informative choices. Data
are simulated from a Gaussian max-stable velocity process with attribute distribution
´1

πpdaq “ πpdr dφ dΛq 9 |Λ|pν´k´1q{2 e´trpΨ Λq{2
k´1
ź „1
qe´qφh 1tφh ą0u `
ˆ
2
h“1

r´3{2 e´papr´mq

2

q{p2m2 rq


1 ´qφh
qe
1tφh ă0u dr dφ dΛ,
2

where pr, φ1 , . . . , φk´1 q is the polar parametrization of the velocity v. This corresponds to a Wishart
distribution for Λ with degrees of freedom ν and shape Ψ, an inverse Gaussian distribution for
the magnitude of the velocity r with parameters m, a, and wrapped Laplace distributions with
parameter q for the angles pφ1 , . . . , φk´1 q defining the direction of the velocity in Rk . The storm
lifetimes τ „ Exponentialpδq and the support points ξ follow a homogeneous spatial Poisson process
with rate βdξ. The intensity measure in the specification of the process u 9 u´2 is improper.
For the simulation, we put u „ Paretopumin , 1q, the conditional distribution of u given u ą umin .
Hyperparameters for the simulation are shown in Table 1.
Table 1: Hyperparameter choices for simulations

value

β
1/600

umin
1

δ
1/120

Ψ
I2

ν
7

m
1/10

a
1/2

q
1/2

The data are simulated on a 10 ˆ 10 box B. In order to inform hyperparameter choices, this box
is taken to roughly represent a 5 000km2 area containing the continental United States and southern
Canada. As a result, the distributions of velocity and lifetimes of storms are chosen to approximate
the behavior of weather events in this geographic region. To set the time scale and allow easier
interpretation of results, one unit of time in the simulation is considered one hour. Support points
of the marginal process N pdξ dσq are sampled on B ˆ r0, T s, with T “ 50 ˆ 365 ˆ 24 “ 438 000, so
that the number of support points of N pdξ dσq are Poisson distributed with mean β ˆ T ˆ 100. The
13

Tail waiting times

14

choice of β “ 1{600 gives an average of four storms a day forming in the region. Storm lifetimes τj
are sampled for each support point from Exponentialp1{120q, which gives an average storm lifetime
of five days. Intensities, shapes, and velocities are sampled from the specified distributions with the
hyperparameters given in Table 1. The values m “ 0.1, φ “ 1{2 for the inverse Gaussian distribution
on r gives an average speed of about 30 miles per hour. The parameter q “ 0.5 places most of the
mass on easterly storm tracks.
The value of the process Y px, tq is recorded at one million homogeneously spaced time points
from r0, T s at the five locations tx1 “ p5, 5q, x2 “ p5, 5.5q, x3 “ p1, 1q, x4 “ p8, 8q, x5 “ p3, 5qu, as well
as twenty additional locations sampled uniformly on B. The five fixed points should result in the
process at some pairs of locations being highly tail dependent, some pairs weakly dependent, and
some nearly independent. After simulation, waiting times between exceedances of y “ Fbi´1 p0.99q
and y “ Fbi´1 p0.999q, where Fbi´1 p¨q is the empirical distribution function of yi , are calculated at the
every unique pair of points. Mixture models of the form in (16) are then fit to the waiting times
until first exceedance and the waiting times between exceedances at all pairs of locations.
Figure 2 shows results; some additional results are provided in Figure S2 in Supplementary
Materials. The posterior concentrates around a single exponential component in the models for κi ,
which is consistent with the marginal distribution when v “ 0 almost surely. Since most pairs of
points are more distant than the product of storm velocity times lifetime, this near-independence
result is sensible. Two components, and occasionally three, are required to model κ∆ . At both
thresholds, γ
bd has a prominent mode near 0.05, and regardless of the choice of d, γ
bd decreases with
the distance between points. When d “ dK , for distances greater than about 2 000 kilometers, γ
bd is
approximately 0.05. Notably, γ
bd decays more slowly and exhibits substantially more variation when
d “ TV, suggesting possibly lower power for tests based on that metric.

Figure 2: Left: histograms of γ
bd for all pairs of locations with d “ dK for thresholds yi “ Fbi´1 p0.999q and
´1
Fbi p0.99q. Center: plots of γ
bd for d “ dK and d “ TV for threshold yi “ Fbi´1 p0.999q versus Euclidean
distance between points. Right: Posterior estimates of ηi | κi and η∆ | κ∆ for y “ F ´1 p0.999q. Box plots
are over all recorded points or all pairs of points.

5

Applications

The method is applied to four real data sets: (1) Daily precipitation data for 25 weather stations in
the United States for the period 1940–2014; (2) Daily exchange rates for 12 currencies for the period
14

J. E. Johndrow and R. L. Wolpert

15

Figure 3: Examples of raw data wpx, tq. Dotted lines: more extreme threshold; dashed lines: less extreme
threshold for computing γd .

1986–1996; (3) Daily prices for the period 2000–2014 of the 30 stocks that made up the Dow Jones
Industrial Index as of January 2015; and (4) Electrical potential at 62 single neurons in the brain of
a mouse exploring a maze, sampled at 500 Hz.
The first dataset is similar to the simulation study and the physical heuristic that we introduced
in Section 2. The points in the max-stable velocity process can be thought of as storms or pressure
systems in this context. The fourth— the mouse electrophysiology data— retains an explicit spatial
component, as the physical distance between neurons is related to the speed at which signals can be
propagated between them. The second and third applications are financial. One way to interpret
the max stable velocity model in financial settings is to think of assets as embedded in latent space,
with the distances between them reflecting similarity in the factors that determine their price. The
support points of the process can be thought of as market sentiment, and their movement reflects the
evolution of sentiment through the asset space. Larger values of u reflect panics— sentiments that
affect many asset classes simultaneously— and large velocities reflects rapidly spreading sentiments.
Finally, in financial settings it is usually extremes in the left tail that are of particular interest, and
thus it is useful to model with max-stable velocity processes the negative of the observed data— so
exceedances of the negative of a low threshold are the relevant events.
In choosing thresholds for analysis, the heuristic used is that at least 100 first exceedances of
the highest threshold must be observed at all points. Thresholds correspond to empirical quantiles
of the observed data, and at least two thresholds are analyzed for every dataset. A consequence
of choosing thresholds in this way is that for some datasets, the empirical quantile of the chosen
threshold may be much more extreme than for others. For example, the electrophysiology data has
nearly two million observations, so we can choose thresholds of yi “ Fbi´1 p0.998q or yi “ Fbi´1 p0.99q
for analysis; the exchange rate data has only about two thousand observations, so the thresholds
chosen are yi “ Fbi´1 p0.05q and yi “ Fbi´1 p0.10q. For the Dow Jones data, the thresholds are
yi “ Fbi´1 p0.025q and Fbi´1 p0.05q, and for the precipitation data, the thresholds are yi “ Fbi´1 p0.99q
and yi “ Fbi´1 p0.95q. Examples of single components of the four datasets are shown in Figure 3. In
the case of the two financial datasets, the displayed series and data used for our analysis are the log
daily returns logrwpx, tq{wpx, t ´ 1qs, following standard practice in finance. The other two figures
show the raw data plotted in the time domain.
15

Tail waiting times

16

For comparison, we use the method described in [14]. When W is in the max domain-of-attraction
of a Brown-Resnick process, Engelke et al. [14] show that when the process is transformed to exponential margins, correlations in the increments of the process relative to its value any one location,
conditional on threshold exceedance at the reference location, have approximately a Gaussian distribution. Practically, this approach proceeds by first transforming the data to have exponential
margins. Label the transformed data w̃px, tq, then choose a reference location and label it location
1. Compute the increments w̃pxr´1s , tq ´ w̃px1 , tq relative to location 1, where w̃pxr´1s , tq is the
vector of process realizations at time t at locations other than location 1. Retain only data points
where w̃px1 , tq exceeds a high threshold. Now compute the correlations in w̃pxr´1s , tq ´ w̃px1 , tq. The
magnitude of the correlation provides inference on the strength of tail dependence; the closer the
absolute correlation is to 1 or -1, the stronger the tail dependence. We refer to this as the conditional
correlation method.

5.1

Precipitation

Results for analysis of the precipitation data are summarized in Figure 4. A map showing the
location of each station can be found in Figure S1 in Supplementary Materials. Clear geographic
structure is evident in the estimated values of γd . Similar geographic structure is seen in tail waiting
times and conditional correlations. Overall, tail dependence is evident at nearby sites but decays
with distance; for distances greater than about 500 km the estimated values of γd are all very similar.
The number of exponential components in the models for κ∆ and κi is almost always two, with κi
having higher median weight on the second exponential component. This departs somewhat from
the expected mixture of a single exponential component and an atom for a max-stable velocity
process, but the discrete mixture model of (16) is still a good fit to the data.

5.2

Dow Jones components

Figure 5 shows results for analysis of negative log daily returns ´ logrwpx, tq{wpx, t´1qs for the DJIA
data; the reference asset for the conditional correlations method is axp (American Express). Similar
dependence structure is evident using the conditional correlation and tail waiting times methods.
Most values of γ
bd cluster around 0.2 at both thresholds, with a long right tail in the posterior point
estimates. Here, unlike in the precipitation data, the posterior for the parameters of the mixture
models for both κ∆ and κi concentrated around three exponential components, albeit with most
of the weight on the first two. This suggests that the data generating process is somewhat more
complex than a simple max-stable velocity process. In particular, the max-stable velocity process
is stationary in time, while financial returns often exhibit temporal nonstationarity. About half
the assets pairs involving axp have values of pd ą 0.95. Notably, the largest values of γ
bd for pairs
involving axp are for jpm (JP Morgan Chase) and v (Visa), two other large credit card issuers.
There is considerable dispersion in imputed waiting times, with a median of around eight days. A
table showing the identity of the stock corresponding to each row/column in the colormaps in Figure
5 can be found in Table S2 in Supplementary Materials.

5.3

Exchange Rates

Exchange rate data for twelve currencies is described in [18] and [26]. Similar to the Dow Jones data,
these series are transformed to negative log daily returns (´ logrwpx, tq{wpx, t ´ 1qs) and analyzed
16

J. E. Johndrow and R. L. Wolpert

17

Figure 4: Results for daily precipitation data, d “ dK . Top Left: γ
bd for thresholds yi “ Fbi´1 p0.99q (below
´1
b
the main diagonal) and yi “ Fi p0.95q (above the main diagonal). The sites are arranged by geographic
distance. Top center: conditional correlations relative to site 1 yi “ Fbi´1 p0.99q (below the main diagonal)
bd for both thresholds. Bottom
and yi “ Fbi´1 p0.95q (above the main diagonal). Top right: histograms of γ
left: box plots of posterior means of η∆ | κ∆ (top) and ηi | κi for yi “ Fbi´1 p0.99q. Bottom center: plots
of γ
bd versus geographic distance for both thresholds. Bottom right: box plot of posterior samples of γd for
yi “ ´Fbi´1 p0.99q for pairs that include station 1, gray boxes indicate pd ą 0.95.

at two thresholds: yi “ ´Fbi´1 p0.05q and yi “ ´Fbi´1 p0.10q. A table giving the full name of each
currency and its corresponding row/column in the colormap is provided in Table S1 in Supplementary
Materials. Figure 6 shows some results. Here, there is very clear structure in the pattern of pairwise
dependence, with the European currencies (BEF, FRF, DEM, NLG, ESP, SEK, CHF, and GBP)
showing strong evidence of dependence while the other four currencies (AUD, CAD, JPY, and NZD)
show little evidence of tail dependence among themselves or with the European currencies. Results
for the conditional correlation method were similar (not shown). The posterior intervals for γd —
again, colored gray if pd ą 0.95— show the same pattern, though notably all of the currencies
have pd ą 0.95 for at least one pair. These results are broadly consistent with previous analysis of
dependence at the mean in this dataset [see 26].

5.4

Electrophysiology

Potential data recorded at single neurons in the brain of a mouse interacting with a maze are
described in [13] and analyzed here using two thresholds: yi “ Fbi´1 p0.998q and yi “ Fbi´1 p0.99q. In
electrophysiology, interest lies in modeling dependence between neuron “spikes” at different locations
in the brain. Neuronal voltage spikes indicate transmission of signals along axonal pathways, and
large potentials tend to cluster together in small time windows. These events are referred to as
17

Tail waiting times

18

Figure 5: Results for DJIA components. Top Left: γ
bd for thresholds yi “ Fbi´1 p0.025q (below the main
´1
b
diagonal) and yi “ Fi p0.05q (above the main diagonal). Top center: correlations conditional on exceedance
of yi “ Fbi´1 p0.025q (below the main diagonal) and yi “ Fbi´1 p0.05q (above the main diagonal) by axp. Top
right: histograms of γ
bd for both thresholds. Bottom left: posterior means of η∆ | κ∆ (top) and ηi | κi for
yi “ Fbi´1 p0.99q. Bottom center: posterior samples of γd for yi “ ´Fbi´1 p0.025q for pairs that include axp;
gray boxes indicate pd ą 0.95. Bottom right: samples of log imputed posterior waiting times (all pairs pooled)
for exceedance of yi “ ´Fbi´1 p0.025q; vertical line at ten days.

“spike trains.” Thus, in these data one expects to see extensive tail dependence, but the waiting
times between spikes at different neurons are relevant, since they inform about the pathway that
the signal takes through the brain.
The neurons are assigned to regions of the brain, which are shown in some of the subimages in
Figure 7. There is ample evidence of strong tail dependence for most pairs of neurons. Four neurons
show a markedly different pattern of dependence from the others: at times, electrodes are faulty,
resulting in data quality problems, which is the most likely explanation. This anomaly aside, strong
dependence is evident both from waiting times and conditional correlations approach. However, a
pattern revealed by the waiting times method that is not as clear with conditional correlations is that
γd increases markedly with y for almost every neuron pair. Moreover, clear differences in dependence
between and within regions are evident using the waiting times approach, but less clear with the
conditional correlation method. Mixture models for κi required three to four exponential components
and those for κ∆ between two and three, again suggesting that the data generating process may be
nonstationary in time or otherwise more complex than a generic max-stable velocity process (for
which two exponential components suffice). Posterior credible intervals for γd at y “ Fbi´1 p0.998q
have pd ą 0.95 for all pairs involving the first Accumbens region neuron; this was largely the case
for all 1891 pairs of neurons, indicating strong evidence of tail dependence across all brain regions.
The electrophysiology data are overall most similar to the Dow Jones data: values of γd are
18

J. E. Johndrow and R. L. Wolpert

19

Figure 6: Results for exchange rate data. Left: γ
bd for thresholds of (the negative of ) yi “ ´Fbi´1 p0.05q (below
´1
b
the main diagonal) and yi “ ´Fi p0.10q (above the main diagonal). Right: box plot of posterior samples of
γd for yi “ ´Fbi´1 p0.025q for all pairs of currencies. White boxes indicate pd ą 0.95

mostly large and increase as the threshold becomes more extreme. However, in these data the
shift toward one in γd as y increases is much more prominent, with γd appearing almost uniformly
distributed at yi “ Fbi´1 p0.99q but massed near 1 at yi “ Fbi´1 p0.999q. Again, this shift is only evident
using the waiting times method; it is not shown by conditional correlations.

6

Discussion

Characterizing tail dependence based on waiting times between peaks over thresholds has the advantage of greater flexibility and generality than existing alternatives in cases where temporal lags
in extreme events are possible. The method relies strictly on the waiting times and inference on
the parameter γd is relatively simple and computationally scalable, particularly when closed-form
K
expressions are available for dpµ∆ , µK
∆ q. The method has strong theoretical grounding in the class
of max-stable velocity processes, which is sufficiently general to apply to distinct applications where
extremal dependence is of interest. Broader classes of processes with non-trivial tail waiting times
and classical tail dependence certainly exist. Modifying the max-stable velocity process to include
time-varying magnitudes and velocities would create a much richer class of processes that retains
the basic characteristics of the simpler process described here. Moreover, mixtures of the process
could be devised to better model temporally nonstationary data that is common in financial and
neuroscience applications.
The inferential method based on waiting times is robust to misspecification of the underlying
process as long as the method selected for estimation of the waiting time distributions is sufficiently
flexible. Here we have employed a Bayesian approach to estimating these densities based on discrete
mixtures of exponentials and an atom at zero. This model appears to be sufficiently flexible in all of
the applications considered here. It is not hard to conceive of fully nonparametric alternatives for
situations in which the observed waiting time distributions are more complex, with the disadvantage
of a possible loss of power.
Like other peaks-over-thresholds methods, our approach requires the choice of appropriate thresholds. Substantial work has been done on threshold choice for standard peaks-over-thresholds methods. It is unclear whether this will translate directly to threshold choice in this novel context. Here,
we have taken the simpler approach of choosing multiple thresholds and attempting to infer a general pattern as the threshold changes. Further work on threshold choice is undoubtedly called for.
19

Tail waiting times

20

Figure 7: Results for electrophysiology data with d “ dK . Top left: γ
bd for thresholds yi “ Fbi´1 p0.998q (below
´1
b
the main diagonal) and yi “ Fi p0.99q (above the main diagonal). Top center: correlations in increments
relative to Accumbens 1, conditional on exceedance of yi “ Fbi´1 p0.998q (below the main diagonal) and
yi “ Fbi´1 p0.99q (above the main diagonal) at Accumbens 1. Top right: histograms of γd for both thresholds
across all pairs of neurons. Bottom left: η∆ | κ∆ and ηi | κi for y “ Fbi´1 p0.998q. Bottom right: posterior
samples of γd for yi “ Fbi´1 p0.998q for all pairs that include Accumbens 1.

Additionally, we have thus far modeled the pairwise waiting times entirely independently; clear gains
in estimation efficiency would result from sharing information across all pairs. These are promising
areas for future work.

Acknowledgments
The authors thank Dan Cooley and David Dunson for helpful comments on drafts of this manuscript,
and Kristian Lum for assistance in obtaining data for the climatological applications. James
Johndrow acknowledges funding from the United States National Institute of Environmental Health
Sciences, grant numbers ES017436, ES020619, and ES017240. Robert Wolpert acknowledges funding
from United States National Science Foundation, grant numbers SES-1521855 and ACI-1550225.
20

J. E. Johndrow and R. L. Wolpert

21

Appendix
A
A.1

Proof of Theorem 2.4 and Corollaries 2.6 and 2.5
Marginal Distribution of Y px, tq

Fix x P Rd , t P R, and y ą 0. Then the event
rY px, tq ď ys “ N pAq “ 0u
that Y px, tq does not exceed y is identical to the event that the Poisson random measure N pdωq
assigns zero points to the set
`
˘
(
A “ ω : σ ă t ă σ ` τ, uϕΛ x ´ ξ ´ pt ´ σqv ą y
of particles born before time t, surviving until time t, that move from their birth point ξ at velocity
v to a location ξt :“ ξ ` pt ´ σqv sufficiently close to x that their intensity u will lead to exceedance
of level y. We will use νp¨q to denote the measure of a set with respect to the intensity measure of
N . The probability of the event PpN pAq “ 0q for the Poisson measure N „ Popνpdωqq is
`
˘
PrY px, tq ď ys “ exp ´ νpAq
ż
ż
ż
´
¯
“ exp ´
u´2 du βδe´δτ dτ dσ dξπpdv dΛq
Rd ˆP d ˆRd σătăσ`τ uąy{ϕΛ px´ξt q
ż
¯
´ ż
1
ϕΛ px ´ ξt q βδe´δτ dτ dσ dξπpdv dΛq
“ exp ´
Rd ˆP d ˆRd σătăσ`τ y
´ ż
¯
`
˘
β
“ exp ´
ϕΛ x ´ ξ ´ pt ´ σqv dξπpdv dΛq
Rd ˆP d ˆRd δy
´
β¯
“ exp ´
,
δy
so Y px, tq „ Frp1, β{δq has the unit Fréchet distribution with scale β{δ for all locations x P Rd and
times t ą 0, for any probability distribution πpdv dΛq.

A.2

Marginal Distribution of κ and Y ˚ px, tq

Now suppose that ϕpzq is a monotonically decreasing function of squared Euclidean length, and
denote by ξs :“ ξ ` vps ´ σq the location at time s P R of the particle ω “ pu, ξ, σ, τ, Λ, vq, born at
location ξ at time σ. The event that the first exceedance time κpyq of level y ą 0 is later than any
specified t ą 0 is the event that the Poisson random measure N pdωq assigns zero points to the Borel
set
!
)
A :“ ω :
sup
uϕΛ px ´ ξs q ą y ,
(18)
sPr0,tsXrσ,σ`τ s

whose probability is Prκpyq ą ts “ expp´νpAqq. It is convenient for us to write A as the disjoint
union of several simpler pieces, and then sum their measures to find νpAq. For fixed Λ P P d and
21

Tail waiting times

22

v P Rd , denote by ∆ the time interval (positive or negative) between birth and arrival at the closest
point of approach to x (CPA), starting from ξ and traveling at velocity v, i.e.,
(
∆ :“ argmax ϕΛ px ´ ξ ´ svq
sPR
(
“ argmin px ´ ξ ´ svq1 Λpx ´ ξ ´ svq
sPR

“ v 1 Λpx ´ ξq{v 1 Λv.
Any vector z P Rd can be written uniquely as the sum z “ pz } ` z K q of its projections onto the
vector space spanned by the velocity vector v, and its orthogonal complement, in the inner product
Λ, given by:
z } :“ pv 1 Λz{v 1 Λvqv

z K :“ z ´ pv 1 Λz{v 1 Λvqv

(19)

The projection of px ´ ξq onto the span of v is px ´ ξq} “ ∆v, and ξ ` ∆v is the CPA to x.
The particle is initially approaching x if v 1 Λpx ´ ξq ą 0, i.e., ∆ ą 0; otherwise it is initially
receding from x. The supremum in (18) will be attained at time s˚ “ σ ` ∆, if that is in the interval
r0, ts X rσ, σ ` τ s. If not, the supremum will be attained at one of the endpoints of that interval.
We write A as the disjoint union of five sets, one for each possible time
(
s˚ “ argmax Y px, sq : s P r0, ts X rσ, σ ` τ s
at which Y px, sq attains its maximum in r0, ts: the beginning s “ 0 or end s “ t of the interval, the
time of the particle’s birth s “ σ or death s “ σ ` τ , or some intermediate point.
A1 : s˚ P p0, tq X pσ, σ ` τ q. The particle initially approaches x (i.e., ∆ ą“ 0) and reaches
CPA
‰
before its death or time t, whichever occurs first. The supremum sups ϕΛ px ´ ξs q occurs at
CPA time s˚ “ σ ` ∆, at which time ξs˚ “ ξ ` ∆v so
px ´ ξs˚ q “ x ´ ξ ´ ∆v “ px ´ ξqK ,
the projection of px ´ ξq onto the orthogonal complement of v. It follows that
`
˘(
A1 “ ω : uϕΛ px ´ ξs˚ q ą y,
p0 _ σq ď pσ ` ∆q ď t ^ pσ ` τ q

(20)

After integrating wrt u´2 du,
ż
1
νpA1 q “
1
βδe´δτ ϕΛ px ´ ξs˚ q dτ dσdξπpdv dΛq
y Rd ˆP d ˆRd ˆp´∆ďσďt´∆qˆp∆ďτ ă8q t∆ą0u
Integrating wrt τ , then σ,
ż
`
˘
β
“ t
1t∆ą0u e´δ∆ ϕΛ px ´ ξqK dξπpdv dΛq
y
Rd ˆP d ˆRd
1
1
}
Changing variables from ξ to p∆, ζq with ∆ “
? v Λpx ´ ξq{v Λv (so ∆v “ px ´ ξq ) and
K
K
d´1
1
ζ “ px ´ ξq P v ” R , with Jacobian |v| “ v v, and integrating wrt ∆,
ż
β
“
t
ϕΛ pζq dζ |v| πpdv dΛq.
(21)
δy
Rd ˆP d ˆtζPv K u

22

J. E. Johndrow and R. L. Wolpert

23

A2 : s˚ “ t. The particle “approaches‰x and survives until time t but but fails to reach CPA by t.
The supremum sups ϕΛ px ´ ξs q occurs at time s˚ “ t, at which time ξs˚ “ ξ ` pt ´ σqv:
`
˘
A2 “ ω : uϕΛ x ´ ξ ´ pt ´ σqv ą y,

(
p0 _ σq ď t ă σ ` pτ ^ ∆q

(22)

After integrating wrt u´2 du, and noting that px ´ ξs˚ q “ px ´ ξqK ` p∆ ´ pt ´ σqqv so the dξ
integral extends over the half-space on which v 1 Λpx ´ ξs˚ q ą 0,
ż
1
1
βδe´δτ ϕΛ px ´ ξs˚ q dξ dτ dσ πpdv dΛq
y Rd ˆP d ˆp´8,tsˆpt´σ,8qˆRd t∆ąt´σu
ż
1
“
βδe´δτ dτ dσ πpdv dΛq
2y Rd ˆP d ˆp´8,tsˆpt´σ,8q
ż
1
“
βδe´δτ dτ dσ
2y p´8,tsˆpt´σ,8q
β
.
(23)
“
2δy

νpA2 q “

A3 : s˚ “ σ ` τ . The particle lives to time“ 0 and approaches
x but dies before time t without
‰
reaching the CPA. The supremum sups ϕΛ px ´ ξs q occurs at death time s˚ “ σ ` τ , at which
time ξs˚ “ ξ ` τ v, so
`
˘
A3 “ ω : uϕΛ x ´ ξ ´ τ v ą y,

p0 _ σq ď σ ` τ ď t ^ pσ ` ∆q

(

(24)

After integrating wrt u´2 du, and noting that px ´ ξs˚ q “ px ´ ξqK ` p∆ ´ τ qv so the dξ integral
extends over the half-space on which v 1 Λpx ´ ξs˚ q ą 0,
νpA3 q “
“
“
“

ż
1
1
βδe´δτ ϕΛ px ´ ξ ´ τ vq dξ dσ dτ πpdv dΛq
y Rd ˆP d ˆp0,8qˆp´τ,t´τ qˆRd t∆ąτ u
ż
1
βδe´δτ dσ dτ πpdv dΛq
2y Rd ˆP d ˆp0,8qˆp´τ,t´τ q
ż
1
βδe´δτ dσ dτ
2y p0,8qˆp´τ,t´τ q
β
δt.
(25)
2δy

A4 : s˚ “ 0. The particle is born before time zero, and by time zero is still alive and is receding
from x (either because it receded initially, or because it passed the CPA before time zero).
Because this system is invariant under time-reversal, this set has the same ν measure as the
set A2 , giving νpA4 q “ βp2δyq´1 .
A5 : s˚ “ σ. The particle was born during the interval r0, ts and recedes from x. Again appealing
to time-reversibility, νpA5 q “ νpA3 q “ βp2δyq´1 .
23

Tail waiting times

24

A.3

Summary

The ν measures of the five pieces are:
ż
β
t
νpA1 q “
ϕΛ pζq dζ |v| πpdv dΛq
δy
Rd ˆP d ˆtζPv K u
β
β
β
νpA3 q “
δt
νpA4 q “
νpA2 q “
2δy
2δy
2δy

νpA5 q “

β
δt
2δy

whose sum is
β
νpAq “
δy

#

+

ż
1 ` δt ` t

ϕΛ pζq dζ |v| πpdv dΛq .

(26)

Rd ˆP d ˆtζPv K u

Thus, the probability distribution of the first time κpyq that Y px, tq exceeds any level y is of the
form
¯
´
`
˘
β
´ t ¨ const ,
(27)
Prκpyq ą ts “ exp ´ νpAq “ exp ´
δy
a mixture of a point-mass at zero of magnitude 1 ´ exppβ{δyq and an exponentially-distributed
random variable whose rate constant
"
*
ż
β
δ`
ϕΛ pζq dζ |v| πpdv dΛq
δy
Rd ˆP d ˆv K
depends on the mean particle speed and kernel shape. Corollary 2.5 follows immediately since
rY ˚ px, tq ă ys and rκpyq ą ts are identical events.

A.4

Special Case: Gaussian Kernel

Consider the Gaussian case of ϕpzq “ p2πq´d{2 expp´z 1 z{2q, and fix any Λ P P d and v P Rd . The
projection of ξ P Rd onto the span of v (see (19)) is ξ } “ pv 1 Λξ{v 1 Λvqv “ ∆v, where ∆ “ pv 1 Λξ{v 1 Λvq,
with squared Λ-length pξ } q1 Λpξ } q “ ∆2 v 1 Λv. Using (38)
ż
`
˘1
ϕΛ pζq dζ |v| “ v 1 Λv{2π 2 ,
vK

and νpAq from (26) is
β
νpAq “
δy

B

"
ż
1 ` δt ` t

`

1

v Λv{2π

˘ 21

*
πpdv dΛq .

(28)

Rd ˆP d

Proof of Theorem 2.7

Applying the inclusion-exclusion principle and change-of-variables from ξ to z :“ pξ ´ vσq to (8a)
gives the result.
24

J. E. Johndrow and R. L. Wolpert

C

25

Proof of Theorem 2.11

In this section we compute a lower bound for the probability of two exceedances Y pxj , sj q ą yj , j “
1, 2 of levels tyj u Ă R` at locations txj u Ă Rd within the interval tsj u Ă r0, ts for t ą 0. Some
calculations used in the sequel are found in Section D.

C.1

General case

We consider the case in which CPA is attained at both points during the interval r0, ts for the general
max-stable velocity process. For j “ 1, 2 set:
∆j :“

tϕΛ pxj ´ ξ ´ svqu “ v 1 Λpxj ´ ξq{v 1 Λv

argmax
sPR

s˚j :“

argmax

tϕΛ pxj ´ ξ ´ ps ´ σqvqu

sPr0,tsXrσ,σ`τ s

u˚j :“

yj
.
ϕΛ pxj ´ ξs˚ q
j

To achieve exceedance and CPA at both locations a particle must satisfy
u ą u˚1 _ u˚2

s˚j “ σ ` ∆j P r0, ts X rσ, σ ` τ s,

j “ 1, 2.

(29)

Write the set A Ă Ω of particles that satisfy these conditions as the union of four sets
Aij :“ tω P Ω : ∆i “ p∆1 ^ ∆2 q,

u˚j “ pu˚1 _ u˚2 qu

i, j “ 1, 2

(30)

characterized by which CPA occurs first and which exceedance requires the larger mass u. For
example,
A11 “ tω P Ω : u ą u˚1 ą u˚2 , p0 _ σq ă σ ` ∆1 ă σ ` ∆2 ă t ^ pσ ` τ qu
consists of the particles that initially approach both x1 and x2 and reach CPA for both before their
death. Both suprema sups rϕΛ px ´ ξs qs occur at the CPA times s˚j “ σ ` ∆j , at CPA locations
ξs˚ “ ξ ` ∆j v, so
j

pxj ´ ξs˚ q “ pxj ´ ξ ´ ∆j vq “ pxj ´ ξqK
j

is the projection of pxj ´ ξq onto the orthogonal complement of v for j “ 1, 2. After integrating wrt
u´2 du, we have
ż
`
˘
1
1t0ă∆1 ă∆2 u 1tu˚ ąu˚ u βδe´δτ ϕΛ px1 ´ ξqK dτ dσ dξ πpdv dΛq
νpA11 q “
Rd ˆP d ˆRd
1
2
y1 ˆp´∆1 ăσăt´∆2 q
ˆp∆2 ăτ ă8q

Integrating with respect to τ then σ,
“

β
y1

ż
pt ´ p∆2 ´ ∆1 qq e´δ∆2 ϕΛ ppx1 ´ ξqK q

ˆ

Rd ˆP d ˆRd

1t0ă∆1 ă∆2 u 1t∆2 ´∆1 ďtu 1tu˚ ąu˚ u dξ πpdv dΛq.
1

25

2

(31)

Tail waiting times

26

Change variables from ξ P Rd to pζ, γq P v K ˆ v } with
ζ :“

` x1 ` x2
2

´ξ

˘K

,

γ :“

` x1 ` x2
2

´ξ

˘}

and introduce the quantity
∆12 :“ p∆2 ´ ∆1 q “ v 1 Λpx2 ´ x1 q{v 1 Λv,
noting that it doesn’t depend on ξ. With this we can write
∆1 “ ´ 12 ∆12 ` v 1 Λγ{v 1 Λv

∆2 “ ` 21 ∆12 ` v 1 Λγ{v 1 Λv.

Introduce µ :“ 12 px2 ´ x1 qK and note that
`
˘
ϕΛ px1 ´ ξs˚ q “ ϕΛ px1 ´ ξqK “ ϕΛ pζ ´ µq
j
˘
`
ϕΛ px2 ´ ξs˚ q “ ϕΛ px2 ´ ξqK “ ϕΛ pζ ` µq

(32a)
(32b)

j

The limits imposed by the indicator functions in (31) are:
0 ď ∆2 ´ ∆ 1 ď t ô
0 ď v 1 Λpx2 ´ x1 q ď tv 1 Λv
0 ď ∆1
ô
∆12 v 1 Λv ď 2v 1 Λγ
˚
˚
u2 ď u1
ô ϕΛ pζ ´ µq{ϕΛ pζ ` µq ď y1 {y2
Rewriting (31) with this variable change, then integrating wrt γ, gives

νpA11 q “

β
y1

ż

β
“
δy1

1

1

pt ´ ∆12 q e´δ∆12 {2´δv Λγ{v Λv ϕΛ pζ ´ µq

ˆ

Rd ˆP d ˆv K ˆv }

1t0ďv1 Λpx2 ´x1 qďt v1 Λvu 1t∆12 v1 Λvď2v1 Λγu 1tϕΛ pζ´µq{ϕΛ pζ`µqďy1 {y2 u dγ dζ πpdv dΛq
ż
pt ´ ∆12 q e´δ∆12 1t0ďv1 Λpx2 ´x1 qďt v1 Λvu ˆ
Rd ˆP d
+
#ż
ζPv K :

ϕΛ pζ´µq
y1
ϕΛ pζ`µq ď y2

ϕΛ pζ ´ µq dζ

|v| πpdv dΛq.

(33)

Now, we find the measure of the other three sets. First
A21 “ tω P Ω : u ą u˚2 ą u˚1 , p0 _ σq ă σ ` ∆1 ă σ ` ∆2 ă t ^ pσ ` τ qu ,
giving
νpA21 q “

β
δy2

ż
Rd ˆP d

#ż

pt ´ ∆12 q e´δ∆12 1t0ďv1 Λpx2 ´x1 qďt v1 Λvu
+

ζPv K :

ϕΛ pζ´µq
y1
ϕΛ pζ`µq ě y2

ϕΛ pζ ` µq dζ

ˆ

|v| πpdv dΛq.

Finally, it is clear that with
A12 “ tω P Ω : u ą u˚1 ą u˚2 , p0 _ σq ă σ ` ∆2 ă σ ` ∆1 ă t ^ pσ ` τ qu
26

(34)

J. E. Johndrow and R. L. Wolpert

27

A22 “ tω P Ω : u ą u˚2 ą u˚1 , p0 _ σq ă σ ` ∆2 ă σ ` ∆1 ă t ^ pσ ` τ qu
we have
νpA12 q “

β
δy1

ż
Rd ˆP d

#ż

pt ` ∆12 q eδ∆12 1t0ď´v1 Λpx2 ´x1 qďt v1 Λvu
+

ζPv K :

β
νpA22 q “
δy2

ϕΛ pζ ´ µq dζ

|v| πpdv dΛq.

ż
Rd ˆP d

#ż

pt ` ∆12 q eδ∆12 1t0ď´v1 Λpx2 ´x1 qďt v1 Λvu
+

ζPv K :

C.2

ϕΛ pζ´µq
y1
ϕΛ pζ`µq ď y2

ˆ

ϕΛ pζ´µq
y1
ϕΛ pζ`µq ě y2

ϕΛ pζ ` µq dζ

ˆ

|v| πpdv dΛq.

Gaussian case

Now take ϕpzq “ p2πq´d{2 expp´z 1 z{2q and fix v, Λ P Rd ˆ P d . The set of ζ P v K over which the
bracketed integral in (33) is taken can be written as:
y1 `
ϕΛ ζ ` µq
y2
y1
´ 21 pζ ` µq1 Λpζ ` µq
´ 12 pζ ´ µq1 Λpζ ´ µq ď log
y2
1
y1
µ1 Λζ ď log
2
y2
ϕΛ pζ ´ µq ď

Writing ζ P v K as the sum ζ “ ζK ` ζ} of components orthogonal and parallel to px2 ´ x1 qK (in the
Λ metric), by (41),
ż
´ S K pvq logpy {y q ¯
a
β
1 2
`
πpdv dΛq.
νpA11 q “
pt ´ ∆12 q e´δ∆12 v 1 Λv{2πΦ ´ Λ
d
d
y1 R ˆP
2
SΛK pvq
0ď∆12 ďt

with SΛK pvq :“ tpx2 ´ x1 q1 Λpx2 ´ x1 qK u1{2 . The other three sets have measure
ż
´ S K pvq logpy {y q ¯
a
β
1 2
´δ∆12
1 Λv{2πΦ ´ Λ
νpA21 q “
v
´
pt
´
∆
q
e
πpdv dΛq.
12
y2 Rd ˆP d
2
SΛK pvq
0ď∆12 ďt
ż
´ S K pvq logpy {y q ¯
a
β
1 2
δ∆12
1 Λv{2πΦ ´ Λ
νpA12 q “
pt
`
∆
q
e
v
`
πpdv dΛq.
12
Rd ˆP d
y1
2
SΛK pvq
0ď´∆12 ďt
ż
´ S K pvq logpy {y q ¯
a
β
1 2
δ∆12
1 Λv{2πΦ ´ Λ
νpA22 q “
πpdv dΛq.
pt
`
∆
q
e
v
´
12
Rd ˆP d
y2
2
SΛK pvq
0ď´∆12 ďt

Recognizing a simple change of variables, we have
"
ˆ
˙
ż
β
1
SΛK pvq logpy1 {y2 q
´δ∆12
νpAq “ ?
pt
´
∆
q
e
Φ
´
`
12
Rd ˆP d
y1
2
SΛK pvq
2π
Xt´tď∆12 ďtu

27

Tail waiting times

28

`

˙*
ˆ
?
1
S K pvq logpy1 {y2 q
v 1 Λvπpdv dΛq
Φ ´ Λ
´
K
y2
2
SΛ pvq

where A “ A11 Y A12 Y A21 Y A22 . So
Pr|κ2 ´ κ1 | ą ts ď 1 ´ Prκ1 _ κ2 ď ts
ď 1 ´ Pr CPA Exceedances at x1 , x2 in r0, ts s
(
ď exp ´ νpAq ,
”?
´ a
¯ı
so the probability Pr|κ2 ´ κ1 | ą ts Ñ 0 as Eπ
v 1 ΛvΦ ´ px2 ´ x1 q1 Λpx2 ´ x1 qK Ñ 8.

D

Some Important Integrals

As before, fix v, x1 , x2 , ξ P Rd and Λ P P d , and set
∆j :“ v 1 Λpxi ´ ξq

∆12 :“ p∆2 ´ ∆1 q “ v 1 Λpx2 ´ x1 q

µ :“ 21 px2 ´ x1 qK

where, as before, for any z P Rd we denote the projections of z parallel and orthogonal to v in the
Λ metric by
z } :“ pv 1 Λz{v 1 Λvqv
zK “ z ´ z}
and the value of the kernel function at z by
`
˘
1
ϕΛ pzq :“ |Λ{2π| 2 exp ´ z 1 Λz{2
Then
ż
`
˘
1
exp ´ 21 ξ 1 Λξ dξ “ |Λ{2π| 2
Rd
ż
ż
` 1 1 ˘
`
˘
exp ´ 2 γ Λγ dγ “ exp ´ s2 v 1 Λv{2 |v| ds
v}

(35)

R

with the CoV γ “ sv with Jacobian dγ “

?

v 1 v ds
1

“ pv 1 Λv{2πq´ 2 |v|

(36)

ż
`

vK

˘

exp ´ 12 ζ 1 Λζ dζ “ the ratio of (35)/(36)
1

1

“ |Λ{2π| 2 pv 1 Λv{2πq 2 {|v|, so

(37)

ż
1

ϕΛ pζq dζ “ pv 1 Λv{2πq 2 {|v|.
ż
`
˘
`
˘
exp ´ 12 w1 Λw dw “ exp ´ s2 µ1 Λµ{2 |µ| ds

(38)

vK

ż
µ}

R
1

“ pµ1 Λµ{2πq´ 2 |µ|
28

(39)

J. E. Johndrow and R. L. Wolpert

29

ż
tv,µuK

`
˘
exp ´ 12 q 1 Λq dq “ the ratio of (37)/(39)
1

1

“ |Λ{2π| 2 pv 1 Λv µ1 Λµq 2 {2π|v| |µ|, so
ż
1

ϕΛ pζq dζ “ pv 1 Λv µ1 Λµq 2 {2π|v| |µ|.

(40)

tv,µuK

If the integral in (39) extends only over those ζ P µ} with µ1 Λζ ď
`
q´2µ1 Λµ ˘
?2 1
by a factor of Φ logpy12{y
, leading to
µ Λµ
ż
vK

ϕΛ pζ ´ µq 1tµ1 Λζě 12 logpy1 {y2 qu dζ “ Φ

1
2

logpy1 {y2 q, its value is reduced

´ 1 log y1 ´ µ1 Λµ ¯´ v 1 Λv ¯ 21
2
?y2 1
.
2πv 1 v
µ Λµ

(41)

References
[1] Balkema, G. and Embrechts, P. (2007). High risk scenarios and extremes: a geometric approach.
European Mathematical Society.
[2] Ballani, F. and Schlather, M. (2011). A construction principle for multivariate extreme value
distributions. Biometrika, 98(3):633–645.
[3] Beguerı́a, S., Angulo-Martı́nez, M., Vicente-Serrano, S. M., López-Moreno, J. I., and El-Kenawy,
A. (2011). Assessing trends in extreme precipitation events intensity and magnitude using nonstationary peaks-over-threshold analysis: a case study in northeast Spain from 1930 to 2006.
International Journal of Climatology, 31(14):2102–2114.
[4] Beguerı́a, S. and Vicente-Serrano, S. M. (2006). Mapping the hazard of extreme rainfall by
peaks over threshold extreme value analysis and spatial regression techniques. Journal of applied
meteorology and climatology, 45(1):108–124.
[5] Beirlant, J., Goegebeur, Y., Segers, J., and Teugels, J. (2006). Statistics of extremes: theory and
applications. John Wiley & Sons.
[6] Buishand, T. A., de Haan, L., and Zhou, C. (2008). On spatial extremes: with application to a
rainfall problem. The Annals of Applied Statistics, pages 624–642.
[7] Coles, S., Bawa, J., Trenner, L., and Dorazio, P. (2001). An introduction to statistical modeling
of extreme values, volume 208. Springer.
[8] Coles, S. G. and Tawn, J. A. (1991). Modelling extreme multivariate events. Journal of the
Royal Statistical Society. Series B (Methodological), 53(2):377–392.
[9] Das, B. and Resnick, S. I. (2011). Conditioning on an extreme component: Model consistency
with regular variation on cones. Bernoulli, 17(1):226–252.
[10] Davison, A. C. and Smith, R. L. (1990). Models for exceedances over high thresholds. Journal
of the Royal Statistical Society. Series B (Methodological), 52(3):393–442.
[11] de Haan, L. (1984). A spectral representation for max-stable processes. The Annals of Probability, 12(4):1194–1204.
29

Tail waiting times

30

[12] de Haan, L. and Ferreira, A. (2006). Extreme value theory: an introduction. Springer Science
& Business Media.
[13] Dzirasa, K., Phillips, H. W., Sotnikova, T. D., Salahpour, A., Kumar, S., Gainetdinov, R. R.,
Caron, M. G., and Nicolelis, M. A. L. (2010). Noradrenergic control of cortico-striato-thalamic
and mesolimbic cross-structural synchrony. The Journal of Neuroscience, 30(18):6387–6397.
[14] Engelke, S., Malinowski, A., Kabluchko, Z., and Schlather, M. (2014). Estimation of Hüsler–
Reiss distributions and Brown–Resnick processes. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 77(1):239–265.
[15] Ghosal, S., Ghosh, J. K., and Van Der Vaart, A. W. (2000). Convergence rates of posterior
distributions. Annals of Statistics, 28(2):500–531.
[16] Gretton, A., Borgwardt, K. M., Rasch, M., Schölkopf, B., and Smola, A. J. (2006). A kernel
method for the two-sample-problem. In Advances in neural information processing systems, pages
513–520.
[17] Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. (2012). A kernel
two-sample test. The Journal of Machine Learning Research, 13(1):723–773.
[18] Harrison, J. and West, M. (1999). Bayesian Forecasting & Dynamic Models. Springer.
[19] Heffernan, J. E. and Resnick, S. I. (2007). Limit laws for random vectors with an extreme
component. The Annals of Applied Probability, 17(2):537–571.
[20] Heffernan, J. E. and Tawn, J. A. (2004). A conditional approach for multivariate extreme values
(with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology),
66(3):497–546.
[21] Heffernan, J. E., Tawn, J. A., and Zhang, Z. (2007). Asymptotically (in) dependent multivariate
maxima of moving maxima processes. Extremes, 10(1-2):57–82.
[22] Katz, R. W., Parlange, M. B., and Naveau, P. (2002). Statistics of extremes in hydrology.
Advances in water resources, 25(8):1287–1304.
[23] Meinguet, T. (2012). Maxima of moving maxima of continuous functions. Extremes, 15(3):267–
297.
[24] Méndez, F. J., Menéndez, M., Luceño, A., and Losada, I. J. (2006). Estimation of the long-term
variability of extreme significant wave height using a time-dependent peak over threshold (POT)
model. Journal of Geophysical Research: Oceans, 111(C7).
[25] Minsker, S., Srivastava, S., Lin, L., and Dunson, D. B. (2014). Robust and scalable Bayes via
a median of subset posterior measures. arXiv preprint arXiv:1403.2660.
[26] Prado, R. and West, M. (2010). Time series: modeling, computation, and inference. CRC
Press.
[27] Resnick, S. I. and Roy, R. (1991). Random usc functions, max-stable processes and continuous
choice. The Annals of Applied Probability, pages 267–292.
30

J. E. Johndrow and R. L. Wolpert

31

[28] Rootzén, H. and Tajvidi, N. (2006). Multivariate generalized Pareto distributions. Bernoulli,
12(5):917–930.
[29] Rousseau, J. and Mengersen, K. (2011). Asymptotic behaviour of the posterior distribution in
overfitted mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710.
[30] Schlather, M. (2002). Models for stationary max-stable random fields. Extremes, 5(1):33–44.
[31] Smith, R. L. (1984). Threshold methods for sample extremes. In Statistical extremes and
applications, pages 621–638. Springer.
[32] Smith, R. L. (1989). Extreme value analysis of environmental time series: an application to
trend detection in ground-level ozone. Statistical Science, 4(4):367–377.
[33] Smith, R. L. (1990). Max-stable processes and spatial extremes. Unpublished manuscript,
University of Surrey.
[34] Smith, R. L. and Weissman, I. (1996). Characterization and estimation of the multivariate
extremal index. UNC working papers.
[35] Smola, A., Gretton, A., Song, L., and Schölkopf (2007). Hilbert space embeddings for distributions. In Proceedings of the 18th International Conference on Algorithmic Learning, pages 18–31.
Springer.
[36] Song, L., Huang, J., Smola, A., and Fukumizu, K. (2009). Hilbert space embeddings of conditional distributions with applications to dynamical systems. In Proceedings of the 26th Annual
International Conference on Machine Learning, pages 961–968. ACM.
[37] Sriperumbudur, B. K., Gretton, A., Fukumizu, K., Schölkopf, B., and Lanckriet, G. R. G.
(2010). Hilbert space embeddings and metrics on probability measures. The Journal of Machine
Learning Research, 11:1517–1561.
[38] Tawn, J. A. (1988). Bivariate extreme value theory: models and estimation. Biometrika,
75(3):397–415.
[39] Tawn, J. A. (1990). Modelling multivariate extreme value distributions. Biometrika, 77(2):245–
253.
[40] Zhang, Z. and Smith, R. L. (2010). On the estimation and application of max-stable processes.
Journal of Statistical Planning and Inference, 140(5):1135–1153.

31

Tail waiting times

32

Supplementary Materials
S1
S1.1

Proof of Theorems 2.2 and 2.3
Proof of Theorem 2.2

Theorem 2.2 is a corollary of the following result.
iid

Theorem S1.1 (Maxima of MSV Processes). If Yi „ MSVpβ, δ, π, ϕq are independent MSV processes for 1 ď i ď n then their maximum _1ďiďn Yi px, tq :“ max1ďiďn Yi px, tq is a max-stable velocity
process with the MSVpnβ, δ, π, ϕq distribution.
piq

Proof. Each of tYi u has a representation (4a) for the support points tωj ujPN of the ith of n inde`
˘
iid
pendent Poisson random measures N piq pdωq „ Po νpduq for the intensity measure νpdωq of (3).
`
˘
ř
Then their sum N ` pdωq :“ 1ďiďn N piq pdωq „ Po nνpduq is also a Poisson random measure, with
intensity measure n times νpdωq and with support equal to the union of the supports of tN piq u. The
MSV processes associated with N ` is precisely _1ďiďn Yi px, tq, since
P r_1ďiďn Yi px, tq ă ys “ P rXni“1 tYi ă yus
n
ź
“
e´νpBq “ e´nνpBq ,
i“1

for B “ tYi ă yu. But the measure n¨ν of (3) is exactly that of the MSVpnβ, δ, π, ϕq process.

We now prove Theorem 2.2.
iid

Proof. Fix Y „ MSVpβ, δ, π, ϕq and n iid max-stable copies tYi u1ďiďn „ MSVpβ, δ, π, ϕq. By Theorem S1.1 the maximum _1ďiďn Yi px, tq has the MSVpn β, δ, π, ϕq distribution. By (8) in Theorem 2.7, multiply β and each yj by n to see that all finite-dimensional marginal joint distributions
D 1
n

of t_1ďiďn Yi px, tqu are identical to those of tn Y px, tqu, so Y px, tq “
of Theorem 2.1 and Y is max-stable.

S1.2

_1ďiďn Yi px, tq satisfies (1)

Proof of Theorem 2.3

We first prove two lemmas used in obtaining the main result.
D

Lemma S1.1. If Y1 „ MSVpβ, δ, π, ϕq and Yc „ MSVpcβ, δ, π, ϕq for some c ą 0, then cY1˚ “ Yc˚ .
Proof.
Fix txi u1ďiďn Ă Rd and tti u1ďiďn Ă R. The joint CDF for tY1˚ pxi , ti qu at tyi u Ă R` is
“
‰
`
˘
P X1ďiďn rY1˚ pxi , ti q ď yi s “ exp ´ νpY1ďiďn Ai q
32

(S.1)

J. E. Johndrow and R. L. Wolpert

33

for the exceedance events Ai :“ rY1˚ pxi , ti q ą yi s, which may be written
,
$
/
’
.
&
yi
`
˘( .
Ai “ ω : u ą
’
sup
ϕΛ xi ´ ξ ´ vps ´ σq /
%
p0_σqďsďpti ^σ`τ q

`
˘
By the inclusion-exclusion principle, ν Y Ai can be evaluated as
ÿ
`
˘
`
˘
p´1q1`|J| ν XjPJ Aj ,
ν Y1ďiďn Ai “

(S.2)

H‰JĂt1,¨¨¨ ,nu

where |J| denotes the cardinality of J. For a finite set J Ă t1, ¨ ¨ ¨ , nu of indices, the intersection
$
»
fi,
’
/
&
.
č
yj
—
ffi
`
˘(
Aj “ ω : u ą max –
fl
jPJ
’
/
sup
ϕΛ xj ´ ξ ´ vps ´ σq
%
jPJ
p0_σqďsďptj ^σ`τ q

has ν-measure (after changing variables from ξ to z :“ ξ ` vσ)
»
˜
ν

¸
č

Aj

ż
“

sup

— p0_σqďsďptj ^σ`τ q
min –
jPJ

jPJ

`
˘( fi
ϕΛ xj ´ vs ´ z
ffi
fl dz βdσ δe´δτ dτ πpdv dΛq

yj

(S.3)

`
˘
The measure ν XjPJ Aj in (S.3)
if both β and each yj are multiplied by the same
` is unchanged
˘
constant c ą 0. By (S.2), the ν Y1ďiďn Ai is also unchanged.
Finally, by (S.1) and (S.3),
“
‰
“ č
‰
P X1ďiďn rc Y1˚ pxi , ti q ď yi s “ P
rY1˚ pxi , ti q ď yi {cs
1ďiďn

“ č
‰
“P
rYc˚ pxi , ti q ď yi s ,
1ďiďn

proving the lemma.
iid

Lemma S1.2. If tYi u1ďiďn „ MSVpβ, δ, π, ϕq are independent MSV processes with maximal processes
Yi˚ px, tq :“ sup Yi px, sq,
tą0
0ďsďt

_1ďiďn Yi˚ px, tq

then the maximum
Y „ MSVpnβ, δ, π, ϕq.

is the maximal process Y ˚ px, tq for a max-stable velocity process

Proof.
piq
Each of tYi u has a representation (4a) for the support points tωj ujPN of the ith of n independent
`
˘
iid
Poisson random measures N piq pdωq „ Po νpduq for the intensity measure νpdωq of (3). The
`
˘
ř
sum N ` pdωq :“ 1ďiďn N piq pdωq „ Po nνpduq is also a Poisson random measure, with intensity
33

Tail waiting times

34

piq

measure n νpdωq and with support equal to the union tωj u of the supports of tN piq u, and the maxstable velocity process Y px, tq associated with N ` pdωq by (4a) is _1ďiďn Yi px, tq. By Theorem S1.1,
Y „ MSVpnβ, δ, π, ϕq. But _1ďiďn Yi˚ px, tq “ Y ˚ px, tq and the Lemma follows.
We now prove Theorem 2.3.
Proof.
iid
Let tYi u1ďiďn „ MSVpβ, δ, π, ϕq and Y „ MSVpβ, δ, π, ϕq be independent MSV processes with
maximal processes
Yi˚ px, tq :“ sup Yi px, sq

Y ˚ px, tq :“ sup Y px, sq,

0ďsďt

tą0

0ďsďt

By Lemma S1.2, there exists a process Y0 „ MSVpn β, δ, π, ϕq whose maximal process Y0˚ is the
D
D
maximum _1ďiďn Yi˚ px, tq. By Lemma S1.1 with c “ n, Y0˚ “ n Y ˚ , i.e., Y ˚ px, tq “ n1 _1ďiďn
Yi˚ px, tq. By (1) in Theorem 2.1, Y ˚ px, tq is max-stable.

S2

Proof of Corollary 2.8

By Theorem 2.7, the joint CDF of Yt1 px1 q and Yt2 px2 q can be found as
`
˘
PrYt1 px1 q ď y1 , Yt2 px2 q ď y2 s “ exp ´ νpA1 Y A2 q
for the sets
`
˘
(
Ai :“ ω : σ ă ti ă σ ` τ, uϕΛ xi ´ ξ ´ pti ´ σqv ą yi
of particles leading to an exceedance
of˘yi at pxi , ti q for i “ 1, 2. In Theorem 2.4, we found νpAi q “
`
β{δyi . To find the required ν A1 Y A2 q “ νpA1 q ` νpA2 q ´ νpA1 X A2 q we need the measure of the
intersection. Putting ξti “ ξ ` pti ´ σqv,
ż
ż
`
˘
ν A1 X A2 q “
P d ˆRd ˆRd σďpt1 ^t2 q
ż
ż
´2
´δτ
dτ βdσ dξ πpdv dΛq
`
˘ u duδe
y
y
τ ěpt1 _t2 q´σ

ż

uě

1
2
ϕΛ px1 ´ξt q _ ϕΛ px2 ´ξt q
1
2

ż

“
P d ˆRd ˆRd

σďpt1 ^t2 q

ˆ

ż
τ ěpt1 _t2 q´σ

ż

ϕΛ px1 ´ ξt1 q ϕΛ px2 ´ ξt2 q
^
y1
y2

˙
δe´δτ dτ βdσ dξ πpdv dΛq

ż

“
P d ˆRd ˆRd

ˆ

σďpt1 ^t2 q

ϕΛ px1 ´ ξt1 q ϕΛ px2 ´ ξt2 q
^
y1
y2

˙
e´δppt1 _t2 q´σq βdσ dξ πpdv dΛq

Set
∆v :“ px2 ´ x1 q ´ pt2 ´ t1 qv

(S.4)
34

J. E. Johndrow and R. L. Wolpert

35

(which doesn’t depend on σ) and change variables from ξ to z :“ ξ ´
ż

ˆ

ż

“
P d ˆRd ˆRd

“

β ´δ|t2 ´t1 |
e
δ

σďpt1 ^t2 q

ϕΛ pz ´ ∆v {2q ϕΛ pz ` ∆v {2q
^
y1
y2

ˆ

ż
P d ˆRd ˆRd

x1 `x2
2

2
` vr t1 `t
´ σs:
2

˙

e´δppt1 _t2 q´σq βdσ dz πpdv dΛq
˙
ϕΛ pz ´ ∆v {2q ϕΛ pz ` ∆v {2q
^
dz πpdv dΛq,
(S.5)
y1
y2

so then
„

β 1
1
`
δ y1
y2
ˆ
˙
ż
β ´δ|t2 ´t1 |
ϕΛ pz ´ ∆v {2q ϕΛ pz ` ∆v {2q
´ e
^
dz πpdv dΛq
δ
y1
y2
P d ˆRd ˆRd

νpA1 Y A2 q “

In the Gaussian case with ϕpzq “ p2πq´d{2 expp´z 1 z{2q this simplifies. Set B1 :“ A1 X A2 X tω :
∆1 Λz ă log yy21 u, the set of ω P A1 X A2 where the minimum in (S.5) is attained at px1 , t1 q, and set
B2 “ pA1 X A2 qzB1 , the set on which it is attained at px2 , t2 q. Then
ż
β ´δ|t2 ´t1 |
νpB1 q “
e
ϕΛ pz ´ ∆v {2q1t∆1 Λzălog yy1 u dz πpdv dΛq
2
δy1
P d ˆRd ˆRd
For fixed v let ζ :“ z ´ p∆1v Λz{∆1v Λ∆v q∆v be the orthogonal projection onto the space ∆K
v perpendicular to ∆v in the Λ norm, and change variables from z to pζ, sq with z “ ζ ` s∆v and Jacobian
1
dz “ p∆1v ∆v q 2 dζ ds, to find
ż
`
˘ 2 1
β ´δ|t2 ´t1 |
1
νpB1 q “
e
|Λ{2π| 2 exp ´ ζ 1 Λζ{2 e´s ∆v Λ∆v {2
δy1
d
d
K
P ˆR ˆ∆v ˆR
a
1tsălog yy1 {∆1 Λ∆u ∆1v ∆v ds dζ πpdv dΛq
2
ż
¯
´ log y1
β ´δ|t2 ´t1 |
y2
´ 12 SΛ pvq πpdv dΛq
“
e
Φ
δy1
SΛ pvq
P d ˆRd
where SΛ pvq :“

a
∆1v Λ∆v . Similarly,

β ´δ|t2 ´t1 |
νpB2 q “
e
δy2

ż
Φ
P d ˆRd

νpA1 X A2 q “ νpB1 q ` νpB2 q “

´ log

y2
y1

SΛ pvq
ż

¯
´ 12 SΛ pvq πpdv dΛq,

"
ˆ
˙
log yy12
β ´δ|t2 ´t1 |
1
e
Φ
´ 12 SΛ pvq `
δ
y1
SΛ pvq
P d ˆRd
ˆ
˙*
y2
log y1
1
` Φ
´ 1 SΛ pvq
πpdv dΛq
y2
SΛ pvq 2

and
νpA1 Y A2 q “

so

β´ 1
1¯
`
´ νpA1 X A2 q
δ y1
y2
35

Tail waiting times

36

¯´ 1
β´
1¯
1 ´ e´δ|t2 ´t1 |
`
δ
y1
y2
ˆ
˙
"
ż
log yy12
1
β ´δ|t2 ´t1 |
1
Φ SΛ pvq ´
`
` e
δ
y1
2
SΛ pvq
P d ˆRd
ˆ
˙*
log yy21
1
1
` Φ SΛ pvq ´
πpdv dΛq.
y2
2
SΛ pvq
“

S3

(S.6)

Proof of Theorem 2.9

For ai , bi ě 0 let κ1 , κ2 be independent nonnegative random variables satisfying
`
˘
Prκi ą ts “ exp ´ ai ´ bi t ,
tą0
and set κ∆ :“ |κ1 ´ κ2 |. Then for t ą 0,
˘
`
Prκ∆ ą ts “ 1 ´ e´a1 e´a2 ´b2 t
˘
`
` 1 ´ e´a2 e´a1 ´b1 t
ż
`
˘
b1 exp ´ a1 ´ b1 x1 ´ a2 ´ b2 px1 ` tq dx1
`
R
ż `
`
˘
b2 exp ´ a2 ´ b2 x2 ´ a1 ´ b1 px2 ` tq dx2
`

κ1 “ 0 ă t ă κ2
κ2 “ 0 ă t ă κ1
0 ă κ1 ă κ1 ` t ă κ2
0 ă κ2 ă κ2 ` t ă κ1

R`

or, combining terms,
`
˘
Prκ∆ ą ts “ 1 ´ e´a2 e´a1 ´b1 t

`
˘
` 1 ´ e´a1 e´a2 ´b2 t

b2 e´a1 ´a2 ´b1 t
b1 e´a1 ´a2 ´b2 t
e
`
e
b `b
b `b
" 1 ´a21 ´a2
*
" 1 ´a12´a2
*
b1 e
b2 e
´a1
´a2
´b1 t
´a2
´a1
“
` e p1 ´ e q e
`
` e p1 ´ e q e´b2 t .
b1 ` b2
b1 ` b2

`

Thus κ∆ is distributed as a mixture with three components: a point mass at zero and two exponential
distributions with rates b1 and b2 . From (6) we have
"
*
ż
β
β
ai “
bi “
δ`
ϕΛ pζq dζ |v| ΠpdΛ dvq
δyi
δyi
Rd ˆP d ˆv K

S4

Proof of Theorem 2.10

We first calculate the measures of three sets:
B0 : the set where one or more support points is alive between time 0 and time t that causes an
exceedance of y1 at x1 and y2 at x2 .
36

J. E. Johndrow and R. L. Wolpert

37

B1 : the set where one or more support points is alive between time 0 and time t that causes an
exceedance of y1 at x1 but no support points are alive during this time that cause an exceedance
of y2 at x2 .
B2 : the set where one or more support points is alive between time 0 and time t that causes an
exceedance of y2 at x2 but no support points are alive during this time that cause an exceedance
of y1 at x1 .
These measures will then be used to derive the waiting time distribution with zero velocity. As
usual ω “ pu, ξ, σ, τ, aq P Ω will be used to denote a generic support point, with attribute a “ pv, Λq.
Recall that, since v ” 0 and πpt0u ˆ P d q “ 1, we are abusing notation and writing the marginal
distribution of Λ as “πpdΛq”.

S4.1

Evaluating νpB0 q

The set B0 consists of those ω P Ω for which Y px1 , sq ą y1 and Y px2 , sq ą y2 for some 0 ď s ď t,
leading to a simultaneous exceedance and κ∆ “ 0. The time of the simultaneous occurrence will
necessarily be s “ 0, if Y px1 , 0q ą y1 and Y px2 , 0q ą y2 , and otherwise s “ σ P p0, ts. The probability
of a contemporaneous exceedance in the interval r0, ts is expp´νpB0 qq, with
*
"
y2
y1
_
, ´τ ăσ ďt
νpB0 q “ ν ω : u ą
ϕΛ px1 ´ ξq ϕΛ px2 ´ ξq
ż
“
1tuąr ϕ pxy1 ´ξq _ ϕ pxy2 ´ξq su 1t´τ ăσďtu u´2 du dξ βdσ δe´δτ dτ πpdΛq
Λ 1
Λ 2
Ω
ż
“
1t´τ ăσďtu βdσ δe´δτ dτ ˆ
RˆR`

„


ϕΛ px1 ´ ξq ϕΛ px2 ´ ξq
ˆ
^
dξ πpdΛq
y1
y2
Rd ˆP d
„

ż
ϕΛ px1 ´ ξq ϕΛ px2 ´ ξq
β
^
dξ πpdΛq
“ p1 ` δtq
δ
y1
y2
d
d
R ˆP
„

ż
β
1
1
“ p1 ` δtq
1t ϕΛ pzq ă y1 u ` 1t ϕΛ pzq ă y2 u ϕΛ pzq dz πpdΛq
y2
ϕΛ pz`∆0 q
δ
y2 ϕΛ pz´∆0 q y1
Rd ˆP d y1
ż

for ∆0 “ px2 ´ x1 q as in (S.4), with v “ 0.

S4.2

Evaluating νpB1 q and νpB2 q

The set B1 consists of those ω P Ω for which Y px1 , sq ą y1 for some 0 ď s ď t but Y px2 , sq ď y2 for
all 0 ď s ď t. The probability that this exceedance pattern does not occur is expp´νpB1 qq, with
*
"
y1
y2
νpB1 q “ ν ω :
ăuď
, ´τ ăσ ďt
ϕΛ px1 ´ ξq
ϕΛ px2 ´ ξq
ż
“
1t ϕ pxy1 ´ξq ăuď ϕ pxy2 ´ξq su 1t´τ ăσďtu u´2 du dξ βdσ δe´δτ dτ πpdΛq
Λ 1
Λ 2
Ω
ż
“
1t´τ ăσďtu βdσ δe´δτ dτ ˆ
RˆR`

37

Tail waiting times

38


ϕΛ px1 ´ ξq ϕΛ px2 ´ ξq
ˆ
´
1t ϕΛ px2 ´ξq ă y2 u dξ πpdΛq
y1
ϕΛ px1 ´ξq
y1
y2
Rd ˆP d

„
ż
β
ϕΛ px1 ´ ξq ϕΛ px2 ´ ξq
“ p1 ` δtq
´
1t ϕΛ px2 ´ξq ă y2 u dξ πpdΛq
y1
ϕΛ px1 ´ξq
δ
y1
y2
Rd ˆP d
„

ż
β
1
1
“ p1 ` δtq
1t ϕΛ pz`∆0 q ă y2 u ´ 1t ϕΛ pzq ă y2 u ϕΛ pzq dz πpdΛq.
y1
ϕΛ pzq
δ
y2 ϕΛ pz´∆0 q y1
Rd ˆP d y1
ż

„

Similarly, or by symmetry,
*
y2
y1
ěuą
, ´τ ăσ ďt
ϕΛ px1 ´ ξq
ϕΛ px2 ´ ξq
„

ż
1
β
1
“ p1 ` δtq
1t ϕΛ pz´∆0 q ă y1 u ´ 1t ϕΛ pzq ă y1 u ϕΛ pzq dz πpdΛq.
y2
ϕΛ pzq
δ
y1 ϕΛ pz`∆0 q y2
Rd ˆP d y2

"
νpB2 q “ ν ω :

For Gaussian kernels with ϕpzq9 expp´z 2 {2q more explicit expressions are available:
*
ż "
1 ´ 1
1 ´ 1
β
y1
SΛ ¯
y2
SΛ ¯
` Φ
πpdΛq
νpB0 q “ p1 ` δtq
Φ
log
´
log
´
δ
y1
SΛ
y2
2
y2
SΛ
y1
2
Pd
*
ż "
1 ´ 1
β
1 ´ 1
y2
SΛ ¯
y2
SΛ ¯
νpB1 q “ p1 ` δtq
´ Φ
πpdΛq
Φ
log
`
log
´
δ
y1
SΛ
y1
2
y2
SΛ
y1
2
Pd
*
ż "
1 ´ 1
β
y1
SΛ ¯
1 ´ 1
y1
SΛ ¯
νpB2 q “ p1 ` δtq
Φ
log
`
´ Φ
log
´
πpdΛq
δ
y2
SΛ
y2
2
y1
SΛ
y2
2
Pd
*
ż "
1 ´ 1
1 ´ 1
β
y2
SΛ ¯
y1
SΛ ¯
` Φ
πpdΛq
νpB` q “ p1 ` δtq
Φ
log
`
log
`
δ
y1
SΛ
y1
2
y2
SΛ
y2
2
Pd
a
where SΛ :“ px2 ´ x1 q1 Λpx2 ´ x1 q (this would be SΛ p0q, for v “ 0 in the notation introduced in
the proof of Corollary 2.8).
For equal thresholds y1 “ y2 “ y, these have the simple form
ż
`
˘(
1β
1β
p1 ` δtq
2Φ ´ SΛ {2 πpdΛq “
p1 ` δtq pZq
(S.7)
νpB0 q “
yδ
y
δ
d
żP
`
˘
(
1β
1β
νpB1 q “ νpB2 q “
2Φ SΛ {2 ´ 1 πpdΛq “
p1 ` δtq
p1 ` δtqp1 ´ Zq
yδ
y
δ
d
P
ż
`
˘(
1β
1β
νpB` q “
p1 ` δtq
2Φ SΛ {2 πpdΛq
“
p1 ` δtqp2 ´ Zq
yδ
yδ
Pd
for
ż
`
˘
2Φ ´ SΛ {2 πpdΛq.

Z :“
Pd

S4.3

Computing the distribution of |κ2 ´ κ1 |

Set
κ1 “ time until first exceedance of y1 at x1
38

(S.8)

J. E. Johndrow and R. L. Wolpert

39

κ2 “ time until first exceedance of y2 at x2
κ˚1 “ time until first exceedance of y1 at x1 without exceeding y2 at x2
κ˚2 “ time until first exceedance of y2 at x2 without exceeding y1 at x1
κ˚0 “ time until first common exceedance,
Note that κ1 “ pκ˚0 ^ κ˚1 q and κ2 “ pκ˚0 ^ κ˚2 q. Set
m˚ :“ min κ˚j
0ďjď2

for the minimum of tκ˚j u, i.e., the time until the first exceedance.
Set λj :“ νpBj q{p1 ` δtq for j “ 0, 1, 2, so
νpBj q “ λj p1 ` δtq.
Then λ0 ` λ1 “ δyβ1 and λ0 ` λ2 “ δyβ2 . In this notation, the distributions of κ1 and κ2 (given in
(7)) and of m˚ are given for t ě 0 by
Prκi ą ts “ e´pλ0 `λj qp1`δtq

Prm˚ ą ts “ e´λ` p1`δtq ,

each a mixture distribution with a point mass at zero and an exponential distribution. Next we
consider the distribution of |κ2 ´ κ1 |.
For t ě 0, the event |κ2 ´κ1 | ą t can occur in four different ways, depending on which exceedance
occurs first and whether or not the first occurs at time 0. First, the exceedances at time 0:
• κ1 “ 0 ă t ă κ2 :
Prκ1 “ 0 ă t ă κ2 s “ Prκ˚1 “ 0 ă t ă pκ˚0 ^ κ˚2 qs
`
˘
“ 1 ´ e´λ1 e´pλ0 `λ2 qp1`δtq

(S.9a)

Prκ2 “ 0 ă t ă κ1 s “ Prκ˚2 “ 0 ă t ă κ˚0 ^ κ˚1 s
`
˘
“ 1 ´ e´λ2 e´pλ0 `λ1 qp1`δtq

(S.9b)

• κ2 “ 0 ă t ă κ1 :

Next, those where the first exceedance m˚ ą 0 occurs after time zero, and the second after an
additional time t.
• 0 ă κ1 ă κ1 ` t ă κ2 :
Pr0 ă κ1 ă κ1 ` t ă κ2 s “ Prp0 ă m˚ q X pκ˚1 “ m˚ q X pκ2 ą m˚ ` tqs
“ Prm˚ ą 0sPrκ˚1 “ m˚ | m˚ ą 0s
ˆ Prpκ˚0 ^ κ˚2 q ą κ˚1 ` t | pκ˚0 ^ κ˚2 q ą κ˚1 “ m˚ ą 0s
λ1 ´λ` ´pλ0 `λ2 qδt
“
e
λ`
39

Tail waiting times

40

• 0 ă κ2 ă κ2 ` t ă κ1 :
Pr0 ă κ2 ă κ2 ` t ă κ1 s “ Prp0 ă m˚ q X pκ˚2 “ m˚ q X pκ1 ą m˚ ` tqs
λ2 ´λ` ´pλ0 `λ1 qδt
“
e
λ`

(S.9c)

Summing the four cases, for t ě 0 we have
“
λ0 ` λ2 ´λ1 ‰
Prκ∆ ą ts “ e´pλ0 `λ2 qp1`δtq 1 ´
e
λ`
“
λ0 ` λ1 ´λ2 ‰
` e´pλ0 `λ1 qp1`δtq 1 ´
e
,
λ`

(S.10)

a mixture of two exponential distributions and a point mass at zero of mass
`
˘
`
λ0 ˘
Prκ∆ “ 0s “1 ´ e´λ0 e´λ1 ` e´λ2 ` e´λ` 1 `
.
λ`

S5

Gibbs sampler for mixture models

Consider the Bayesian model
κ „ q0 δ0 `

K´1
ÿ

qj Exponentialpλj q

(S.11)

j“1
iid

q „ Dirichletpαq,

λj „ Gammapa, bq.

(S.12)

For observed waiting times κi for i P t1, . . . , N u, a Gibbs sampler for this model, including
imputation of continuous waiting times from censored waiting times, cycles through the following
steps
1. For each observation i P t1, . . . , N u, sample a latent variable zi „ Categoricalpq̃q, where
$
qj λj e´λκ̃i
&
ř
ją0
´λj κ̃i
q0 1tκ̃i “0u ` K
j“1 λj e
q̃j “
q0 1tκ̃i “0u
%
řK
j“0
´λ κ̃
q0 1tκ̃i “0u `

j“1

λj e

j i

2. Sample λj from
˜

¸

λj „ Gamma a `

ÿ

1tzi “ju , b `

i

ÿ

κ̃i 1tzi “ju

.

i

3. Sample q from
˜
q „ Dirichlet α `

¸
ÿ

1tzi “1u , . . . , α `

i

ÿ
i

40

1tzi “Ku

J. E. Johndrow and R. L. Wolpert

41

4. Sample wi „ Categoricalpq̃q where
$
qj pFλj pκi `1q´Fλj pκi qq
&
ř
q0 1tκi “0u ` j qj pFλj pκi `1q´Fλj pκi qq
q̃j “
q 1
%
ř 0 tκi “0u
q0 1tκi “0u `

j

qj pFλj pκi `1q´Fλj pκi qq

ją0
j“0

then sample
#
κ̃i „

δ0
wi “ 1
Exponentialrκi ,κi `1s pλj q wi “ j ` 1

where Exponentialrl,us is an exponential distribution truncated to the interval rl, us.

S6

Supplemental Figures

Table S1: Key indicating identity of currencies in Figure 6.
col./row #
1
2
3
4
5
6

symbol
AUD
BEF
CAD
FRF
DEM
JPY

name
Australian Dollar
Belgian Franc
Canadian Dollar
French Franc
German Deutschmark
Japanese Yen

41

col./row #
7
8
9
10
11
12

symbol
NLG
NZD
ESP
SEK
CHF
GBP

name
Dutch Guilder
New Zealand Dollar
Spanish Peseta
Swedish Kroner
Swiss Franc
British Pound

Tail waiting times

42

Table S2: Key indicating identity of stocks in the images in Figure 5
col./row #
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

symbol
axp
ba
cat
csco
cvx
dd
dis
ge
gs
hd
ibm
intc
jnj
jpm
ko

name
American Express
Boeing
Caterpillar
Cisco Systems
Chevron
DuPont
Disney
General Electric
Goldman Sachs
Home Depot
IBM
Intel
Johnson & Johnson
J.P. Morgan Chase
Coca-Cola

col./row #
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

symbol
mcd
mmm
mrk
msft
nke
pfe
pg
t
trv
unh
utx
v
vz
wmt
xom

name
McDonald’s
3M
Merck
Microsoft
Nike
Pfizer
Proctor & Gamble
AT&T
Travelers
United Healthcare
United Technologies
Visa
Verizon
Wal-Mart
Exxon-Mobil

Figure S1: Map with numbers labeling locations of weather stations; the numbers correspond to the order in
which the stations appear in the colormap images in Figure 4.

42

J. E. Johndrow and R. L. Wolpert

43

Figure S2: Left: plots of γ
bd for d “ W1,ϕ and d “ T V for threshold yi “ Fbi´1 p0.99q versus Euclidean distance
between points. Right: Posterior estimates of ηi | κi and η∆ | κ∆ for y “ F ´1 p0.99q.

43

