2.1. Poisson Point Processes

Let Ω be an Euclidean space (or subspace or domain). A random countable collection of points Π ⊂ Ω is called a point process [17]. We denote by Π(B) the number of points of Π residing in the subset B:

Hence, the point process Π ⊂ Ω induces a counting measure on Ω given by (4).

A point process Π ⊂ Ω is said to be a Poisson point process [17] with rate r(ω) if the following pair of conditions hold:

• If B is a subset of Ω, then the random variable Π(B) is Poisson distributed with mean ∫_B r(ω) dω.

• If {B_k}_k is a finite collection of disjoint subsets of Ω, then {Π(B_k)}_k is a finite collection of independent random variables.

The Poisson point process Π can also be described by its finite-dimensional distributions: The probability that a given set of points {ω_j}_j=1ⁿ belongs to Π is

Informally, Ω is divided into infinitesimal cells. Each cell contains either a single point or no points at all. The cells are independent, and the probability that the cell dω contains a point is r(ω) dω.

The best known example of a Poisson point process is the standard Poisson process, where Ω = [0,∞) and the rate function r(ω) is constant.

2.2. The Event Process

Equipped with the notion of Poisson point processes, we can now rigorously define the point process

presented in Section 1. Recall that we considered a generic continuous-time physical system whose “history” is a random collection

of points in the plane—the point (t,x) represents an event of magnitude x taking place at time t. In Section 1, we said, informally, that “events of magnitude belonging to the infinitesimal range (x,x + dx) arrive at rate λ(x) dx” and, “the occurrences of events of different magnitudes are independent.” Put rigorously, the process

(henceforth referred to as our event process) is taken to be a time-homogeneous Poisson point process on

with rate

Thus,

We set Λ(x) to be the Poissonian rate at which samples of size greater than x arrive, namely

The function Λ(x) is smooth, monotone nonincreasing, and fully characterizes the event process

. This function will turn out to play a key role in the sequel. We henceforth refer to the function Λ(x) as the characteristic of the process

and assume that

We will elaborate on this assumption in the following subsection. Furthermore, in Subsection 3.1 we shall show that this assumption is in fact an essential requirement.

Finally, note that the random variable

—counting the number of events of size x ∈ [a,b] occurring during the time interval [t,t + Δ] —is Poisson distributed with mean

2.3. From Discrete to Continuous and Back

In the discrete-time setting described in Section 1, the underlying time series {ξ_n}_n=1^∞ is an i.i.d. sequence of random variables. If we take the random samples {ξ_n}_n=1^∞ to arrive according to some continuous-time counting process (N(t))_t≥0, then the sample set, at time t, would be

(the sample set being empty in the case N(t) = 0). The setting in which (N(t))_t≥0 is a renewal process was introduced in [10] and coined “Random Record Models” (see also [4,22]).

If (N(t))_t≥0 is a standard Poisson process with rate ρ, then the “sample process” (11) is in fact a time-homogeneous Poisson point process on

with rate

where f (x) is the probability density function of the ξ's.

On the other hand, if the rate function λ(x) of the event process

is integrable (i.e., if ∫λ(x) dx < ∞), then

can be described in the form of the sample set (11), where (i) the Poissonian arrival rate is ρ := ∫λ(x) dx and (ii) the probability density function of the ξ's is f (x) := λ(x)/ρ.

Hence, in the case of integrable rate functions, there is a one-to-one correspondence, via embedding, between discrete-time i.i.d. sequences {ξ_n}_n=1^∞ and continuous-time event processes

. However, the truly interesting case is where the rate function is actually not integrable: ∫λ(x) dx = ∞. In this case,

has infinitely many events occurring on all time scales (i.e., on any time interval), and no correspondence between

and any discrete-time i.i.d. sequence {ξ_n}_n=1^∞ can be established. Thus, in the case of a nonintegrable rate function, the event process

is a truly continuous-time “creature.”

In this article, we focus on the “truly continuous-time case,” where the rate functions are nonintegrable. This is the reason for assumption (9). The limit lim_x→−∞ Λ(x) = +∞ ensures that the rate function λ(x) is indeed nonintegrable. The limit lim_x→+∞ Λ(x) = 0, on the other hand, ensures that events of infinitely large magnitude do not occur. We will return to discuss assumption (9) in Subsection 3.1.

3.1. The Maximum

We compute the distribution of the maximum X₁(t). The maximum at time t is less or equal to x if and only if no events of magnitude greater than x occurred during the time period [0,t]; that is,

However, the random variable

is Poisson distributed with mean tΛ(x) (recall (10)) and, hence,

The probability density function of X₁(t) is given, in turn, by

The assumption of (9): From (14), we see that the distribution of the maximum is proper if and only if the assumption of (9) holds. Indeed, in order for the distribution of X₁(t) to be proper, the limit of (14) must equal zero when x → −∞ and must equal one when x → +∞. This, however, takes place if and only if lim_x→−∞ Λ(x) = +∞ and lim_x→+∞ Λ(x) = 0. Hence, the assumption of (9) is not an impeding restriction, but an essential requirement.

3.2. Beating the Maximum

Assume that the current maximum level is x. How long will we have to wait until this maximum level is beat? Since events of magnitude greater than x arrive at rate Λ(x), the answer is straightforward: The waiting time is exponentially distributed with rate Λ(x) (mean 1/Λ(x)). Let us consider the two following variations of this question.

To prove (16), use the probability density function of (15) and the change of variables u = tΛ(x):

3.3. Fréchet, Weibull, and Gumbel

Equation (14) implies that the function Λ(x) fully characterizes the distribution of the maximum. On the other hand, the function Λ(x) also fully characterizes the underlying event process

. Hence, there is a one-to-one correspondence—conveyed by the characteristic Λ(x)—between event processes and the distribution of their maxima. In particular, the event processes corresponding to the “Central Limit Theorem” distributions of EVT (viz. the Fréchet, Weibull, and Gumbel distributions) are given, respectively, by

(the exponent α appearing in the Fréchet and Weibull distributions is a positive parameter).

Could these extreme value laws be obtained (as in the “classic” EVT) as the only possible maxima scaling limits? The answer, as explained below, is affirmative.

The continuous-time scaling of the maximum X₁(t), analogous to the discrete-time scaling of (2), is

where a(t) and b(t) are the scaling functions (a(t) being positive valued). The distribution of the scaled maxima

, using (14), is hence given by

On the other hand, the distribution of the scaled maxima

of (2) is given by

where L(x) := −ln(P(ξ₁ ≤ x)). Moreover, the function L(x) satisfies the very same properties the function Λ(x) does; namely it is monotone decreasing, with lim_x→−∞ L(x) = +∞ and lim_x→+∞ L(x) = 0.

Thus, (17) (as t → ∞) and (18) (as n → ∞) must yield the same distributional limits. The “classic” EVT asserts that the three possible limits of (18) are the extreme value distributions: Fréchet, Weibull, and Gumbel. Hence, these probability laws are also the only possible limits of (17).

The scaling functions in the Fréchet, Weibull, and Gumbel cases are

It is interesting to note that in these three special cases, the scaling yields

.

For a detailed analysis regarding the “basins of attraction” of the Fréchet, Weibull, and Gumbel extreme value distributions, we refer the readers to [3].

3.4. The nth “Runner-up” and the “Top n”

The distribution of X_n+1(t)—the nth “runner-up” is computed analogously to the distribution of the maximum X₁(t). The (n + 1)st order statistic at time t is less or equal to x if and only if no more than n events of magnitude greater than x occurred during the time period [0,t]; that is,

Hence, since

is Poisson distributed with mean tΛ(x), we have

The probability density function of X_n+1(t) is given, in turn, by

Note that (19) and (20) indeed coincide with (14) and (15) in the case n = 0.

Equation (19) gives the one-dimensional distribution of the order statistics. What about the joint, multidimensional, probability distribution of the order statistics; namely the joint distribution of the “top n”—the vector (X₁(t),…,X_n(t))?

Well, the joint probability density function of the vector (X₁(t),···,X_n(t)) is given by

where x₁ > x₂ > ··· > x_n. The explanation follows.

In order for the points x₁ > x₂ > ··· > x_n to be the “top n” points of the sample process

at time t, we need (for j = 1,…,n) the following: (i) an event of magnitude x_j takes place during the time period [0,t]; this occurs with probability tλ(x_j) dx_j; and (ii) no events of magnitude ∈ (x_j,x_j−1) (where x₀ is set to equal +∞) take place during the time period [0,t]; this occurs with probability exp{−t(Λ(x_j) − Λ(x_j−1))}. Multiplying these probabilities together yields the multidimensional density function (21).

4.1. Markovian Structure

The conditional distribution of a (n + 1)st order statistic, given the value of the nth order statistic, is

(y < x). The explanation follows.

Given that the nth order statistic equals x, the (n + 1)st order statistic will be less or equal to y (where y < x) if and only if no events of magnitude ∈ (y,x) occur during the time interval [0,t]; that is, if and only if

. Since

is a Poisson point process, the random variable

is Poisson distributed with mean t(Λ(y) − Λ(x)) and is independent of the points of

residing in [0,t] × [x,∞). Hence, the left-hand side of (22) equals the probability that

, which, in turn, is given by the right-hand side of (22).

Since (22) is a Markovian recursion, it implies that the sequence of order statistics {X_n(t)}_n=1^∞ is a Markov chain (in the variable n, for t fixed). The initial condition of this Markov chain is given by the distribution of the maximum X₁(t); see (14).

4.2. The Hidden Poissonian Structure

When viewed in the proper perspective, the sequence of order statistics {X_n(t)}_n=1^∞ conceals a hidden Poissonian structure. The “proper perspective” is to consider the sequence {Λ(X_n(t))}_n=1^∞, rather than the original sequence {X_n(t)}_n=1^∞. Let us begin with the distribution of Λ(X₁(t)).

Using (14), we have

(u ≥ 0). Hence, Λ(X₁(t)) is exponentially distributed with parameter t (mean 1/t). In a similar way, (19) implies that Λ(X_n(t)) is Gamma distributed with parameters (t,n).

More informative, however, is to analyze the conditional distribution of the increment Λ(X_n+1(t)) − Λ(X_n(t)) given Λ(X_n(t)). Indeed, using (22), we have

that is, the increment Λ(X_n+1(t)) − Λ(X_n(t)) is independent of Λ(X_n(t)) and is exponentially distributed with parameter t (mean 1/t). Hence, we have obtained the following proposition:

Proposition 1: Let

be an event process with characteristic Λ(x) and order statistics {X_n(t)}_n=1^∞. Then the increasing sequence of points {Λ(X_n(t))}_n=1^∞ forms a standard Poisson process with rate t.

4.3. Simulation and Gamma Embedding

Proposition 1 provides us with a remarkably simple simulation algorithm for the entire sequence of order statistics {X_n(t)}_n=1^∞:

where {Z_n}_n=1^∞ is an i.i.d. sequence of exponentially distributed random variables with unit mean. In particular, the simulation algorithms for the Fréchet, Weibull, and Gumbel order statistics are as follows:

Equation (23) is in fact a manifestation of a Gamma embedding as we now explain.

The Gamma process (G(s))_s≥0 is a stochastic process starting at the origin (G(0) = 0) whose increments are independent, stationary, and Gamma distributed:

(b > a ≥ 0). The Gamma process is a special example of one-sided Lévy processes (Lévy subordinators) [1].

At the point s = 1, the value of the Gamma process G(1) is exponentially distributed with unit mean. Hence, since the Gamma process has independent and stationary increments, the increments {G(n) − G(n − 1)}_n=1^∞ form an i.i.d. sequence of exponentially distributed random variables with unit mean. Combining this observation with (23), we obtain the following Gamma embedding of the sequence of order statistics {X_n(t)}_n=1^∞:

that is, the discrete-parameter sequence {X_n(t)}_n=1^∞ is embedded in the continuous-time process (Λ⁻¹(G(s)/t))_s≥0, which, in turn, is a transformation of the Gamma process (G(s))_s≥0.

Furthermore, the random variable Λ⁻¹(G(s)/t), for noninteger parameter values s, may be considered a “virtual” s-order statistic. Continuing yet another step in this direction, we can define the Fréchet, Weibull, and Gumbel virtual order processes as follows (s ≥ 0):

5.1. The Fréchet and Weibull Cases

In the Fréchet case, Λ(x) = x^−α, and in the Weibull case, Λ(x) = |x|^α (the exponent α being a positive parameter). Hence, Proposition 2 implies the following:

Most interesting is the first ratio—X₂(t)/X₁(t) in the Fréchet case and |X₁(t)|/|X₂(t)| in the Weibull case—which measures the relative magnitudes of the “winner” and the “first runner-up.” The ratio distribution, in both the Fréchet and Weibull cases, is governed by the probability density function f (u) = αu^α−1 (0 < u < 1). This density function undergoes a phase transition when crossing the parameter value α = 1: (i) When 0 < α < 1, the density f (u) is monotone decreasing, starting at f (0) = ∞ and decreasing to f (1) = α; (ii) when α = 1 (the “critical” parameter value), the density f (u) is uniform; and (iii) when α > 1, the density f (u) is monotone increasing, starting at f (0) = 0 and increasing to f (1) = α. Hence, the “first runner-up” trails considerably behind the “winner” when α is small, and it follows the “winner” closely when α is large.

At the other end of the order spectrum, the distribution function u^αn becomes more and more concentrated around the value u = 1 as n → ∞. This implies that as n → ∞, the order statistics get closer and closer to each other in ratio. In the Fréchet case, as n → ∞, the order statistics converge to zero and hence get closer and closer to each other, also absolutely. In the Weibull case, on the other hand, the order statistics diverge to −∞ and hence remain apart.

In the Fréchet and Weibull cases, Proposition 2 also provides us with an algorithm for the simulation of the relative contributions of the “top n” order statistics X₁(t),X₂(t),…,X_n(t) to their aggregate magnitude. The explanation follows.

Let {U_k}_k=1ⁿ⁻¹ be an i.i.d. sequence of uniformly distributed random variables, and set {M_k}_k=0ⁿ⁻¹ to be given by

Proposition 2 implies (all equalities below are equalities in law):

5.2. The Gumbel Case

In the Gumbel case Λ(x) = exp{−x} and hence proposition 2 implies:

Alternatively, the random variables {n(X_n(t) − X_n+1(t))}_n=1^∞ are i.i.d. and exponentially distributed with unit mean.

Let us study the following “asymptotically centered” sequence of differences {D_n}_n=1^∞ defined by

Using (30), the sequence {D_n}_n=1^∞ admits the probabilistic representation

where {Z_n}_n=1^∞ is an i.i.d. sequence of exponentially distributed random variables with unit mean. The representation (31), in turn, yields that the mean and variance of D_n are

where γ ≃ 0.577 is Euler's constant.

Moreover, the representation (31) further yields that the Laplace transform of D_n is

(θ ≥ 0). Hence, using Stirling's formula, we have

The function Γ(1 + θ), however, is the Laplace transform of the standard Gumbel distribution (see, for example, [5]), and hence we conclude the following: