Sampling Sup-Normalized Spectral Functions for Brown-Resnick Processes

TL;DR

This paper introduces improved simulation methods for Brown-Resnick processes using spectral functions, optimizing proposal densities for efficiency in MCMC and rejection sampling, with demonstrated performance benefits.

Contribution

It generalizes existing simulation approaches for Brown-Resnick processes by developing new classes of proposal densities and optimizing their efficiency.

Findings

Enhanced simulation efficiency demonstrated in examples

Optimal proposal densities identified for MCMC and rejection sampling

Generalized approaches applicable to spectral functions of max-stable processes

Abstract

Sup-normalized spectral functions form building blocks of max-stable and Pareto processes and therefore play an important role in modeling spatial extremes. For one of the most popular examples, the Brown-Resnick process, simulation is not straightforward. In this paper, we generalize two approaches for simulation via Markov Chain Monte Carlo methods and rejection sampling by introducing new classes of proposal densities. In both cases, we provide an optimal choice of the proposal density with respect to sampling efficiency. The performance of the procedures is demonstrated in an example.

Figures

Equations

$$Z(t){=}_d\underset{i\in \mathbb{N}}{\max }\left\{U_i\cdot V_i(t)\right\},t\in K,$$

$$\displaystyle Z(t)=_{d}\max_{i\in\mathbb{N}}\Big{\{}U_{i}\cdot c_{\infty}\cdot\frac{V^{\max}_{i}(t)}{\|V^{\max}_{i}\|_{\infty}}\Big{\}},\quad t\in K,$$

$$\mathbb{P}(V^{\max }\in B)=c_{\infty }^{-1}\cdot {\int }_B\Vert v{\Vert }_{\infty }\mathbb{P}(V\in \mathrm{d}v),B\in \mathcal{C}(K),$$

$$\mathcal{L}\left(x^{-1}X(\cdot )\mid \Vert X{\Vert }_{\infty }>x\right){⟶}_w\mathcal{L}\left(P\frac{V^{\max }(\cdot )}{\Vert V^{\max }{\Vert }_{\infty }}\right)$$

$$\gamma (h)=\frac{1}{2}\mathbb{E}\left\{(G(t+h)-G(t){)}^2\right\},t,h\in {\mathbb{R}}^d.$$

$$\frac{V^{\max }(\cdot )}{\Vert V^{\max }{\Vert }_{\infty }}=\exp \left(W^{\max }(\cdot )-\Vert W^{\max }{\Vert }_{\infty }\right)$$

$$\displaystyle\mathbf{Q}=\mathbf{C}^{-1}-\frac{\mathbf{C}^{-1}\mathbf{1}_{N}\mathbf{1}_{N}^{\top}\mathbf{C}^{-1}}{\mathbf{1}_{N}^{\top}\mathbf{C}^{-1}\mathbf{1}_{N}}\qquad\text{and}\qquad\mathbf{m}=-\Big{(}\frac{1}{2}\bm{\sigma}+\frac{1-\frac{1}{2}\bm{\sigma}^{\top}\mathbf{C}^{-1}\mathbf{1}_{N}}{\mathbf{1}_{N}^{\top}\mathbf{C}^{-1}\mathbf{1}_{N}}\mathbf{1}_{N}^{\top}\Big{)}\mathbf{C}^{-1}$$

$$\mathbb{P}\left({\mathbf{V}}^{\max }\in S_i\right)=\frac{\det ({\mathbf{Q}}_{-i}{)}^{-1/2}\exp \{\frac{1}{2}{\mathbf{m}}_{-i}^{\top }{\mathbf{Q}}_{-i}^{-1}{\mathbf{m}}_{-i}\}{\Phi }_{N-1}(0_{N-1};{\mathbf{Q}}_{-i}^{-1}{\mathbf{m}}_{-i},{\mathbf{Q}}_{-i}^{-1})}{{\sum }_{j=1}^N\det ({\mathbf{Q}}_{-j}{)}^{-1/2}\exp \{\frac{1}{2}{\mathbf{m}}_{-j}^{\top }{\mathbf{Q}}_{-j}^{-1}{\mathbf{m}}_{-j}\}{\Phi }_{N-1}(0_{N-1};{\mathbf{Q}}_{-j}^{-1}{\mathbf{m}}_{-j},{\mathbf{Q}}_{-j}^{-1})},$$

$$f_{\max }(\mathbf{w})=$$

$$\alpha (\tilde{\mathbf{w}},\mathbf{w})=\min \left\{\frac{f_{\max }(\mathbf{w})/f_{\mathrm{prop}}(\mathbf{w})}{f_{\max }(\tilde{\mathbf{w}})/f_{\mathrm{prop}}(\tilde{\mathbf{w}})},1\right\},$$

$$\Vert P^n({\mathbf{w}}^{(0)},\cdot )-\mathbb{P}({\mathbf{W}}^{\max }\in \cdot ){\Vert }_{\mathrm{TV}}\stackrel{n\to \infty }{⟶}0,$$

$$\alpha (\tilde{\mathbf{w}},\mathbf{w})=\min \left\{\frac{{\max }_{i=1}^N{e}^{w_i}}{{\max }_{i=1}^N{e}^{{\tilde{w}}_i}},1\right\},\mathbf{w},\tilde{\mathbf{w}}\in {\mathbb{R}}^N.$$

$${\mathrm{essinf}}_{\mathbf{w}\in {\mathbb{R}}^N}\frac{f(\mathbf{w})}{f_{\max }(\mathbf{w})}={\mathrm{essinf}}_{\mathbf{w}\in {\mathbb{R}}^N}\left(c_{\infty }\cdot {\min }_{i=1}^N{e}^{-w_i}\right)=0$$

$$f_{\mathrm{prop}}={\sum }_{i=1}^N{p}_i{f}_i$$

$$\tilde{\alpha }(\tilde{\mathbf{w}},\mathbf{w})=\min \left\{\frac{{\max }_{i=1}^N{e}^{w_i}/{\sum }_{i=1}^N{p}_i{e}^{w_i}}{{\max }_{i=1}^N{e}^{{\tilde{w}}_i}/{\sum }_{i=1}^N{p}_i{e}^{{\tilde{w}}_i}},1\right\}.$$

$$\underset{\mathbf{w}\in {\mathbb{R}}^N}{\inf }\frac{f_{\mathrm{prop}}(\mathbf{w})}{f_{\max }(\mathbf{w})}=\underset{\mathbf{w}\in {\mathbb{R}}^N}{\inf }\frac{{\sum }_{i=1}^N{p}_i{e}^{w_i}}{c_{\infty }^{-1}{\max }_{i=1}^N{e}^{w_i}}=c_{\infty }\cdot \underset{i=1}{\overset{N}{\min }}{p}_i>0,$$

$$\displaystyle D(p_{1},\ldots,p_{N})=\int\nolimits_{\mathbb{R}^{N}}\bigg{(}\frac{f_{\rm prop}(\mathbf{w})}{f_{\max}(\mathbf{w})}-1\bigg{)}^{2}\,f_{\max}(\mathbf{w})\,\mathrm{d}\mathbf{w},$$

$$D(p_1,\dots ,p_N)=$$

$$\displaystyle\left(\begin{array}[]{cc}\bm{\Sigma}&\mathbf{1}_{N}\\ \mathbf{1}_{N}^{\top}&0\end{array}\right)\,\left(\begin{array}[]{c}\mathbf{p}\\ \lambda\end{array}\right)=\left(\begin{array}[]{c}\mathbf{0}_{N}\\ 1\end{array}\right),$$

$$\displaystyle\textstyle\sigma_{ik}=\mathbb{E}\big{\{}e^{W(t_{i})+W(t_{k})-\max_{j=1}^{N}W(t_{j})}\big{\}}.$$

$$\mathbf{p}=\frac{{\mathbf{\Sigma }}^{-1}{1}_N}{1_N^{\top }{\mathbf{\Sigma }}^{-1}{1}_N}$$

$$\min$$

$$\text{s.t.}{1}_N^{\top }\mathbf{p}$$

$$p_i$$

$$f_{\max }(\mathbf{w})=c_{\infty }^{-1}\underset{i=1}{\overset{N}{\max }}\exp (w_i)f(\mathbf{w})=\frac{c_{\infty }^{-1}{\max }_{i=1}^N\exp (w_i)}{(2\pi {)}^{N/2}\det (\mathbf{C}{)}^{1/2}}\exp \left\{-\frac{1}{2}{\left(\mathbf{w}+\frac{\mathbf{\sigma }}{2}\right)}^{\top }{\mathbf{C}}^{-1}\left(\mathbf{w}+\frac{\mathbf{\sigma }}{2}\right)\right\},\mathbf{w}\in {\mathbb{R}}^N,$$

$$f_i(\mathbf{w})=\frac{\exp (w_i)}{(2\pi {)}^{\frac{N}{2}}\det (\mathbf{C}{)}^{\frac{1}{2}}}\exp \left\{-\frac{1}{2}{\left(\mathbf{w}+\frac{\mathbf{\sigma }}{2}\right)}^{\top }{\mathbf{C}}^{-1}\left(\mathbf{w}+\frac{\mathbf{\sigma }}{2}\right)\right\}=\frac{1}{(2\pi {)}^{\frac{N}{2}}\det (\mathbf{C}{)}^{\frac{1}{2}}}\exp \left\{-\frac{1}{2}{\left(\mathbf{w}-{\mathbf{C}}_{\cdot i}+\frac{\mathbf{\sigma }}{2}\right)}^{\top }{\mathbf{C}}^{-1}\left(\mathbf{w}-{\mathbf{C}}_{\cdot i}+\frac{\mathbf{\sigma }}{2}\right)\right\},$$

$$\mathcal{L}({\mathbf{W}}^{(i)})=\mathcal{L}(\mathbf{W}+C_{\cdot i})$$

$${\sigma }_{ik}=\mathbb{E}\left\{e^{W(t_i)+W(t_k)-{\max }_{j=1}^{N}W(t_j)}\right\}=$$

$$=$$

$$f_{\max }(\mathbf{w})\le (c_{\infty }\cdot C{)}^{-1}{\tilde{f}}_{\mathrm{prop}}(\mathbf{w}),\text{for all }\mathbf{w}\in {\mathbb{R}}^N,$$

$$C=\underset{\mathbf{w}\in {\mathbb{R}}^N}{\inf }\frac{{\tilde{f}}_{\mathrm{prop}}(\mathbf{w})}{c_{\infty }{f}_{\max }(\mathbf{w})}.$$

Full text

1 Introduction

Spatial and spatio-temporal extreme value analysis aims at investigating extremes of quantities described by stochastic processes. In the classical setting, the real-valued process of interest $X=\{X(t),\,t\in K\}$ is sample-continuous on a compact domain $K\subset\mathbb{R}^{d}$. Analysis of its extremes is often based on results of the limiting behavior of maxima of independent copies $X_{i}$, $i\in\mathbb{N}$. Provided that there exist continuous normalizing functions $a_{n}:K\to(0,\infty)$ and $b_{n}:K\to\mathbb{R}$ such that the process of normalized maxima $\{\max_{i=1}^{n}a_{n}^{-1}(t)\cdot(X_{i}(t)-b_{n}(t)),\,t\in K\}$ converges in distribution to some sample-continuous process $Z$ with nondegenerate margins as $n\to\infty$, the limit process $Z$ is necessarily max-stable and we say that $X$ is in the max-domain of attraction of $Z$.

From univariate extreme value theory, it follows that the marginal distributions of $Z$ are necessarily generalized extreme value (GEV) distributions cf. 4 for instance. As max-stability is preserved under marginal transformations between different GEV distributions, without loss of generality, it can be assumed that $Z$ has standard Fréchet margins, i.e. $\mathbb{P}(Z(t)\leq z)=\exp(-1/z)$, $z>0$, for all $t\in K$. By de Haan 3, any sample-continuous max-stable process with standard Fréchet margins can be represented as

[TABLE]

where the so-called spectral processes $V_{i}$, $i\in\mathbb{N}$, are independent copies of a nonnegative sample continuous stochastic process $V$ on $K$ satisfying $\mathbb{E}\{V_{i}(t)\}=1$ for all $t\in K$, and $\sum_{i\in\mathbb{N}}\delta_{U_{i}}$ is a Poisson point process on $(0,\infty)$ which is independent of the $V_{i}$ and has intensity measure $\Lambda$ given by $\Lambda\{(u,\infty)\}=u^{-1}$ for all $u>0$.

Due to its complex structure, many characteristics of the max-stable process $Z$ in (1) cannot be calculated analytically, but need to be assessed via simulations. In order to simulate $Z$ efficiently, Oesting \BOthers. 15 suggest to make use of the sup-normalized spectral representation

[TABLE]

where the $U_{i}$ are the same as above, the processes $V^{\max}_{i}$ are independently and identically distributed, independently of the $U_{i}$, with distribution $\mathbb{P}(V^{\max}\in\cdot)$ given by

[TABLE]

and $\|f\|_{\infty}=\sup_{t\in K}f(t)$ for every $f\in C(K)$, where $C(K)$ denotes the set of all real-valued continuous functions on $K$ equipped with the supremum norm $\|\cdot\|_{\infty}$ and corresponding $\sigma$-algebra $\mathcal{C}(K)$. Here, the normalizing constant $c_{\infty}=\mathbb{E}\{\|V\|_{\infty}\}$ is the so-called extremal coefficient of the max-stable process $Z$ over the domain $K$. In a simulation study, Oesting \BOthers. 15 demonstrate that simulation based on the sup-normalized spectral representation is competitive to other state-of-the-art algorithms such as simulation based on extremal functions 7 provided that the normalized spectral process $V^{\max}$ can be simulated efficiently.

The law of the processes $V^{\max}$ also occurs when analyzing the extremes of a stochastic process $X$ in an alternative way focusing on exceedances over a high threshold: If $X$ is in the max-domain of attraction of the max-stable process $Z$ in (1), we have

[TABLE]

as $x\to\infty$, where $P$ is a standard Pareto random variable and $V^{\max}$ is an independent process with law given in (3). The limit process $PV^{\max}(\cdot)/\|V^{\max}\|_{\infty}$ is called Pareto process cf. 9, 8.

Arising as sup-normalized spectral process for both max-stable and Pareto processes, the process $V^{\max}$ plays an important role in modelling and analyzing spatial extremes. As a crucial building block of spatial and spatio-temporal models, this process needs to be simulated in an efficient way. Due to the measure transformation in (3), however, sampling of $V^{\max}$ is not straightforward even in cases where the underlying spectral process $V$ can be simulated easily.

In the present paper, we focus on the simulation of $V^{\max}$ for the very popular class of log Gaussian spectral processes, i.e. $V(t)=\exp(W(t))$ for some Gaussian process $W$ such that $\mathbb{E}\{\exp(W(t))\}=1$ for all $t\in K$. The resulting subclass of max-stable processes $Z$ in (1) comprises the only possible nontrivial limits of normalized maxima of rescaled Gaussian processes, the class of Brown–Resnick processes 13, 12. In order to obtain Brown–Resnick processes that can be extended to stationary processes on $\mathbb{R}^{d}$, 13 consider $W(t)=G(t)-\mathrm{Var}\{G(t)\}/2$, $t\in K$, with $G$ being a centered Gaussian process on $\mathbb{R}^{d}$ with stationary increments and variogram

[TABLE]

It is important to note that the law of the resulting max-stable process $Z$ does not depend on the variance of $W$, but only on $\gamma$. Therefore, $Z$ is called Brown–Resnick process associated to the variogram $\gamma$.

Recently, Ho \BBA Dombry 11 introduced a two-step procedure to simulate the corresponding sup-normalized process

[TABLE]

efficiently if the finite domain $K=\{t_{1},\ldots,t_{N}\}$ is of small or moderate size:

Sample the index $i$ of the component where the vector $\mathbf{V}^{\max}=(V^{\max}(t_{k}))_{k=1,\ldots,N}$ assumes its maximum, i.e. select one of the events $\mathbf{V}^{\max}\in S_{i}=\{\mathbf{s}\in(0,\infty)^{N}:\ \|\mathbf{s}\|_{\infty}=s_{i}\}$, $i=1,\ldots,N$. Provided that the covariance matrix $\mathbf{C}$ of the Gaussian vector $\mathbf{W}=(W(t_{k}))_{i=k}^{N}$ is nonsingular, we have that this index is a.s. unique and that the probabilities of the corresponding events can be calculated in terms of the matrix $\mathbf{Q}\in\mathbb{R}^{N\times N}$ and the vector $\mathbf{m}\in\mathbb{R}^{N}$ given by

[TABLE]

where $\bm{\sigma}=(\mathrm{Var}(W(t_{k})))_{k=1,\ldots,N}$ is the variance vector of $\mathbf{W}$ and $\mathbf{1}_{N}=(1,\ldots,1)^{\top}\in\mathbb{R}^{N}$. More precisely, by Ho \BBA Dombry 11,

[TABLE]

where $\mathbf{m}_{-j}$ denotes the vector $\mathbf{m}$ after removing the $j$th component, $\mathbf{Q}_{-j}$ denotes the matrix $\mathbf{Q}$ after removing the $j$th row and $j$th column and $\Phi_{N-1}(\mathbf{0}_{N-1};\bm{\mu},\bm{\Sigma})$ is the distribution function of an $(N-1)$-dimensional Gaussian distribution with mean vector $\bm{\mu}\in\mathbb{R}^{N-1}$ and covariance matrix $\bm{\Sigma}\in\mathbb{R}^{(N-1)\times(N-1)}$ evaluated at $\mathbf{0}_{N-1}=(0,\ldots,0)\in\mathbb{R}^{N-1}$. 2. 2.

Conditional on $\mathbf{V}^{\max}\in S_{i}$, we have $V^{\max}(t_{i})/\|\mathbf{V}^{\max}\|_{\infty}=1$ and the distribution of the vector $\mathbf{M}=(\log(V^{\max}(t_{j})))_{j\neq i}-\log(\|\mathbf{V}^{\max}\|)$ is an $(N-1)$-dimensional Gaussian distribution with mean vector $\mathbf{Q}_{-i}^{-1}\mathbf{m}_{-i}$ and covariance matrix $\mathbf{Q}_{-i}^{-1}$ conditional on $\bm{M}$ being nonpositive.

However, the first step includes computationally expensive operations such as the evaluation of $(N-1)$-dimensional Gaussian distribution functions and the inversion of matrices of sizes $N\times N$ and $(N-1)\times(N-1)$. Furthermore, an efficient implementation of the second step is not straightforward. Thus, the procedure is feasible for small or moderate $N$ only.

In this paper, we will introduce alternative procedures for the simulation of $V^{\max}$, or, equivalently, $W^{\max}=\log V^{\max}$, that are supposed to work for larger $N$, as well. To this end, we will modify a Markov Chain Monte Carlo (MCMC) algorithm proposed by Oesting \BOthers. 15 and a rejection sampling approach based on ideas of de Fondeville \BBA Davison 2. Both procedures have originally been designed to sample sup-normalized spectral functions in general. Here, we will adapt them to the specific case of Brown–Resnick processes.

2 Simulating $W^{\max}$ via MCMC algorithms

Based on the Brown–Resnick process as our main example, we consider a max-stable process $Z$ with spectral process $V=e^{W}$ for some sample-continuous process $W$. Henceforth, we will always assume that the simulation domain $K=\{t_{1},\ldots,t_{N}\}\subset\mathbb{R}^{d}$ is finite and that the corresponding spectral vector $\mathbf{W}$ possesses a density $f$ w.r.t. some measure $\mu$ on $\mathbb{R}^{N}$. Then, by (3), the transformed spectral vector $\mathbf{W}^{\max}=\log(\mathbf{V}^{\max})$, where the logarithm is applied componentwise, has the multivariate density

[TABLE]

which obviously has the same support as $f$, i.e. $\mathrm{supp}(f_{\max})=\mathrm{supp}(f)$.

As direct sampling from the density $f_{\max}$ is rather sophisticated and the normalizing constant $c_{\infty}$ is not readily available, it is quite appealing to choose an MCMC approach for simulation. In the present paper, we focus on Metropolis-Hastings algorithms with independence sampler cf. 17 for example. Denoting the strictly positive proposal density on $\mathrm{supp}(f)$ by $f_{\rm prop}$, the algorithm is of the following form:

Algorithm 1**.**

*MCMC APPROACH (METROPOLIS–HASTINGS)

Input: proposal density $f_{\rm prop}$

Simulate $\mathbf{w}^{(0)}$ according to the density $f_{\rm prop}$.

for $k=1,\ldots,n_{MCMC}$ {

  for Sample $\mathbf{w}$ from $f_{\rm prop}$ and set

  for $\mathbf{w}^{(k)}=\begin{cases}\mathbf{w}&\text{with probability }\alpha(\mathbf{w}^{(k-1)},\mathbf{w}),\\ \mathbf{w}^{(k-1)}&\text{with probability }1-\alpha(\mathbf{w}^{(k-1)},\mathbf{w}),\end{cases}$

  for where the acceptance probability $\alpha(\cdot,\cdot)$ is given by (4).

}

Output: Markov chain $(\mathbf{w}^{(1)},\ldots,\mathbf{w}^{(n_{MCMC})})$.

Here, the acceptance ratio $\alpha(\widetilde{\mathbf{w}},\mathbf{w})$ for a new proposal $\mathbf{w}\in\mathrm{supp}(f)$ given a current state $\widetilde{\mathbf{w}}\in\mathrm{supp}(f)$ is

[TABLE]

using the convention that a ratio is interpreted as [math] if both the enumerator and the denominator are equal to [math]. This choice of $\alpha(\widetilde{\mathbf{w}},\mathbf{w})$ ensures reversibility of the resulting Markov chain $\{\mathbf{w}^{(k)}\}_{k\in\mathbb{N}}$ with respect to the distribution of $\mathbf{W}^{\max}$. Further, it allows for a direct transition from any state $\widetilde{\mathbf{w}}\in\mathrm{supp}(f)$ to any other state $\mathbf{w}\in\mathrm{supp}(f)$. Consequently, the chain is irreducible and aperiodic and, thus, its distribution converges to the desired stationary distribution, that is, for a.e. initial state $\mathbf{w}^{(0)}\in\mathrm{supp}(f)$, we have that

[TABLE]

where $P^{n}(\mathbf{w}^{(0)},\cdot)$ denotes the distribution of the $n$-th state of a Markov chain with initial state $\mathbf{w}^{(0)}$ and $\|\cdot\|_{\rm TV}$ is the total variation norm.

As a general approach for the simulation of sup-normalized spectral processes of arbitrary max-stable processes, Oesting \BOthers. 15 propose to use Algorithm 1 with the density $f$ of the original spectral vector $\mathbf{W}$ as proposal density (Algorithm 1A) and the Metropolis-Hastings acceptance ratio in (4) simplifies to

[TABLE]

As the proposal density $f_{\rm prop}=f$ is strictly positive on $\mathrm{supp}(f)$, convergence of the distribution of the Markov chain to the distribution of $\mathbf{W}^{\max}$ as in (5) is ensured. If the support of $f$ is unbounded, however, there is no uniform geometric rate of convergence of the chain in (5), as we have

[TABLE]

14. In particular, this holds true for the case of a Brown–Resnick process where $\mathbf{W}$ is a Gaussian vector.

Furthermore, due to the structure of the acceptance ratio in (6), the Markov chain may get stuck, once a state $\widetilde{\mathbf{w}}$ with a large maximum $\max_{i=1}^{N}\exp(\widetilde{w}_{i})$ is reached. This might lead to rather poor mixing properties of the chain. Even though independent realizations could still be obtained by starting new independent Markov chains cf. 15, such a behavior is undesirable having chains in high dimension $N$ with potentially long burn-in periods in mind.

While the algorithm in Oesting \BOthers. 15 is designed to be applicable in a general framework, we will use a specific transformation to construct a Markov chain with stronger mixing and faster convergence to the target distribution. For many popular models such as Brown–Resnick processes, this transformation is easily applicable. More precisely, we consider the related densities $f_{i}$, $i=1,\ldots,N$, with $f_{i}(\mathbf{w})=\exp(w_{i})f(\mathbf{w})$. These densities are closely related to the distributions $P_{i}$ that have been studied in Dombry \BOthers. 7.

Hence, we propose to approach the target distribution with density $f_{\max}=c_{\infty}^{-1}\max_{i=1}^{N}f_{i}$ by Algorithm 1 using a mixture

[TABLE]

as proposal density, where the weights $p_{i}\geq 0$, $i=1,\ldots,N$, are such that $\sum_{i=1}^{N}p_{i}=1$. The corresponding acceptance probability in (4) is then

[TABLE]

With the proposal density being strictly positive on $\mathrm{supp}(f)$, we can see that the distribution of the Markov chain again converges to its stationary distribution with density $f_{\max}$. As we further have

[TABLE]

provided that $p_{i}>0$ for $i=1,\ldots,N$, the results found by Mengersen \BBA Tweedie 14 even ensure a uniform geometric rate of convergence for any starting value $\mathbf{w}^{(0)}\in\mathrm{supp}(f)$ in contrast to the case where $f_{\rm prop}=f$.

In order to obtain a chain with good mixing properties, we choose $p_{i}$ such that the acceptance rate in Algorithm 1 is high provided that the current state $\mathbf{w}^{(k)}$ is (approximately) distributed according to the stationary distribution. To this end, we minimize the relative deviation between $f_{\rm prop}$ and $f_{\max}$ under $f_{\max}$, i.e. we minimize

[TABLE]

under the constraint $\sum_{i=1}^{N}p_{i}=1$. Introducing a Lagrange multiplier $\lambda\in\mathbb{R}$, minimizing

[TABLE]

results in solving the linear system

[TABLE]

where $\mathbf{p}=(p_{1},\ldots,p_{N})^{\top}$ and $\bm{\Sigma}=(\sigma_{ik})_{1\leq i,k\leq N}$ with

[TABLE]

Provided that the matrix $\bm{\Sigma}$ is nonsingular, the solution of (16) is given by

[TABLE]

cf. 1 for instance. This solution does not necessarily satisfy the additional restriction $p_{i}\geq 0$ for all $i=1,\ldots,N$. In case that $\Sigma$ is singular or the vector $\mathbf{p}$ has at least one negative entry, the full optimization problem

[TABLE]

has to be solved. Using the Karush–Kuhn–Tucker optimality conditions, the quadratic program (2) can be transformed into a linear program with additional (nonlinear) complementary slackness conditions. It can be solved by modified simplex methods. Alternatively, the problem (2) can be solved by the dual method by Goldfarb \BBA Idnani 10.

Remark 2.1**.**

In order to ensure a geometric rate of convergence of the distribution of the Markov chain, we might replace the condition $p_{i}\geq 0$ for all $i=1,\ldots,N,$ in (2) by $p_{i}\geq\varepsilon$ for some given $\varepsilon>0$. Then, a geometric rate of convergence follows from (9) as described above.

In the case of Brown–Resnick processes, for simplicity, we consider the case that $\mathbf{W}$ possesses a full Lebesgue density. In particular, the covariance matrix $\mathbf{C}=\left(\mathrm{Cov}(W(t_{i}),W(t_{j}))\right)_{1\leq t_{i},t_{j}\leq N}$ of $\mathbf{W}$ is assumed to be nonsingular. Then, the target density is

[TABLE]

where $\bm{\sigma}=(\mathrm{Var}(W(t_{k})))_{k=1,\ldots,N}$ is again the variance vector of $\mathbf{W}$. Now, the densities $f_{i}(\mathbf{w}))=e^{w_{i}}f(\mathbf{w})$ which form the proposal density are just shifted Gaussian distributions:

[TABLE]

cf. Lemma 1 in the Supplementary Material of Dombry \BOthers. 7, i.e. we have

[TABLE]

where the Gaussian vectors $\mathbf{W}^{(i)}$ and $\mathbf{W}$ possess densities $f_{i}$ and $f$, respectively. The calculation of the optimal weights $p_{i}$ is based on the expectation in (17) which typically cannot be calculated analytically, but needs to be assessed numerically via simulations. Such a numerical evaluation, however, is challenging as the random variable $\exp(W(t_{i})+W(t_{k})-\max_{j=1}^{N}W(t_{j}))$ is unbounded. To circumvent these computational difficulties, we make use of the identity

[TABLE]

This expression can be conveniently assessed numerically as the random variable $\exp(W^{(i)}(t_{k})-\max_{j=1}^{N}W^{(i)}(t_{j}))$ is bounded by $1$.

Remark 2.2**.**

Note that both (19) and the final result in (20) still hold true if $\mathbf{W}$ does not possess a full Lebesgue density, but exactly one component $W_{i^{*}}$ is degenerate and the reduced covariance matrix $(C_{ij})_{i,j\neq i^{*}}$ is nonsingular. This situation appears in several examples such as $W$ being a fractional Brownian motion where $W(0)=0$ a.s.

In summary, we propose the procedure below to simulate the normalized spectral vector $\mathbf{W}^{\max}$ for the Brown–Resnick process (Algorithm 1B):

Calculate $\mathbf{p}$ by solving the quadratic program (2) where the entries of the matrix $\bm{\Sigma}$ are given by (20). Provided that all its components are nonnegative, the solution $\mathbf{p}$ has the form (18). 2. 2.

Run Algorithm 1 with proposal density $f_{\rm prop}=\sum_{i=1}^{N}p_{i}f_{i}$ and acceptance probability given by (8). The output of the algorithm is a Markov chain whose stationary distribution is the distribution of $\mathbf{W}^{\max}$.

3 Exact Simulation via Rejection Sampling

In this section, we present an alternative procedure to generate samples from $\mathbf{W}^{\max}$ with probability density $f_{\max}$. In contrast to Section 2 where we generated a Markov chain consisting of dependent samples with the desired distribution as stationary distribution, here, we aim to produce independent realizations from the exact target distribution. To this end, we make use of a rejection sampling approach cf. 5 for instance based on a proposal density $\tilde{f}_{\rm prop}$ satisfying

[TABLE]

for some $C>0$.

Algorithm 2**.**

*REJECTION SAMPLING APPROACH

Input: proposal density $\tilde{f}_{\rm prop}$ and constant $C>0$ satisfying (21)

repeat {

  repeat Simulate $\mathbf{w}^{*}$ according to the density $\tilde{f}_{\rm prop}$.

  repeat Generate a uniform random number $u$ in $[0,1]$.

} until $u\cdot\tilde{f}_{\rm prop}(\mathbf{w}^{*})\leq C\cdot c_{\infty}\cdot f_{\max}(\mathbf{w}^{*})$

Output: exact sample $\mathbf{w}^{*}$ from distribution with density $f_{\max}$

Thus, on average, $(c_{\infty}\cdot C)^{-1}$ simulations from the proposal distribution are needed to obtain an exact sample from the target distribution. Of course, to minimize the computational burden, for a given proposal density $\tilde{f}_{\rm prop}$, the constant $C$ should be chosen maximal subject to (21), i.e.

[TABLE]

Recently, de Fondeville \BBA Davison 2 followed a similar idea and suggested to base the simulation of a general sup-normalized spectral process $V(\cdot)^{\max}/\|V^{\max}\|_{\infty}$ on the relation

[TABLE]

where $\tilde{V}$ is a spectral process normalized with respect to another homogeneous functional $r$ instead of the supremum norm, i.e. $r(\tilde{V})=1$ a.s. If $\|\tilde{V}\|_{\infty}$ is a.s. bounded from above by some constant, from the relation (22), we obtain an inequality of the same type as (21) for the densities of $V(\cdot)^{\max}/\|V^{\max}\|_{\infty}$ and $\tilde{V}(\cdot)/\|\tilde{V}\|_{\infty}$ instead of $f_{\max}$ and $f_{\rm prop}$, respectively. Thus, samples of $\tilde{V}(\cdot)/\|\tilde{V}\|_{\infty}$ can be used as proposals for an exact rejection sampling procedure. For instance, the sum-normalized spectral vector $\tilde{\mathbf{V}}$, i.e. the vector which is normalized w.r.t. the functional $r(f)=\|f\|_{1}=\sum_{k=1}^{N}|f(t_{k})|$, can be chosen as it is easy to simulate in many cases cf. 7 and satisfies $\|\tilde{\mathbf{V}}\|_{\infty}\leq 1$ almost surely.

For a Brown–Resnick process, it is well-known that the sum-normalized process $\tilde{\mathbf{V}}$ has the same distribution as $\exp(\mathbf{W}^{\rm prop})/\|\exp(\mathbf{W}^{\rm prop})\|_{1}$ where $\mathbf{W}^{\rm prop}$ has the density $f_{\rm prop}$ from (7) in Section 2 with equal weights $p_{1}=\ldots=p_{N}=1/N$ see also 6. Thus, in this case, the procedure proposed by de Fondeville \BBA Davison 2 with $r(f)=\|f\|_{1}$ is equivalent to performing rejection sampling for $\mathbf{W}^{\max}$ with $\tilde{f}_{\rm prop}=\sum_{i=1}^{N}\frac{1}{N}f_{i}$ as proposal distribution (Algorithm 2A). From Equation (9), it follows that rejection sampling can also be performed with $\tilde{f}_{\rm prop}=\sum_{i=1}^{N}p_{i}f_{i}$ and arbitrary positive weights $p_{1},\ldots,p_{N}$ summing up to $1$, since we have (21) with $C=\min_{i=1}^{N}p_{i}$. Thus, accepting a proposal $\mathbf{w}^{*}$ in the rejection sampling procedure with probability

[TABLE]

we will obtain a sample of independent realizations from the exact target distribution $f_{\max}$. The rejection rate, however, is pretty high. In order to obtain one realization of $\mathbf{W}^{\max}$, on average $\left(c_{\infty}\cdot\min_{i=1}^{N}p_{i}\right)^{-1}$ simulations from $f^{*}_{\rm prop}$ are required. It can be easily seen that the computational costs are indeed minimal for the choice $p_{1}=\ldots=p_{N}=1/N$, i.e. the choice in the approach based on the sum-normalized representation. In this case, one realization of $\mathbf{W}^{\max}$ on average requires to sample $c_{\infty}^{-1}N$ times from $f_{\rm prop}$. Therefore, this approach becomes rather inefficient if we have a large number $N$ of points on a dense grid.

In order to reduce the large computational costs of rejection sampling which are mainly due to the fact that $\min_{i=1}^{N}f_{i}(\mathbf{w})/f_{\max}(\mathbf{w})$ gets small as $\|\mathbf{w}\|\to\infty$ , we replace each density $f_{i}$ by the modified multivariate Gaussian density $g_{i,\varepsilon}$ whose variance is increased by the factor $(1-\varepsilon)^{-1}\geq 1$ for some $\varepsilon\in[0,1)$:

[TABLE]

Analogously to $f_{\rm prop}$ for the MCMC approach in Section 2, we propose a mixture

[TABLE]

with $p_{i}\geq 0$ and $\sum_{i=1}^{N}p_{i}=1$ as proposal density for the rejection sampling algorithm. A proposal $\mathbf{w}^{*}$ is then accepted with probability

[TABLE]

where

[TABLE]

Thus, to summarize, for appropriately chosen $\varepsilon>0$ and $\mathbf{p}\geq 0$ such that $\|\mathbf{p}\|_{1}=1$, we propose to run Algorithm 2 with proposal density $\tilde{f}_{\rm prop}=\sum\nolimits_{i=1}^{N}p_{i}g_{i,\varepsilon}$ and $C=C(\mathbf{p},\varepsilon)$ according to (3).

Remark 3.1**.**

To further reduce the computational costs in the simulation, we might even choose a more flexible approach. For instance, instead of using a mixture of a finite number of functions $g_{1,\varepsilon},\ldots,g_{N,\varepsilon}$, one could consider arbitrary mixtures

[TABLE]

where $g_{t,\varepsilon}(\mathbf{w})=\int_{\mathbb{R}}g_{N+1,\varepsilon}(\mathbf{w},w_{N+1})\mathrm{d}w_{N+1}$ on the enlarged domain $\{t_{1},\ldots,t_{N},t\}$ and $\nu$ is a probability measure on $\mathbb{R}^{d}$. Furthermore, depending on $t\in\mathbb{R}^{d}$, different values for $\varepsilon=\varepsilon(t)\in[0,1)$ might be chosen. However, due to the complexity of the optimization problems involved, we restrict ourselves to the situation above where $\nu$ is a probability measure on $K=\{t_{1},\ldots,t_{N}\}$ and $\varepsilon$ is constant in space.

Using the procedure described above, on average, $(c_{\infty}\cdot C(\mathbf{p},\varepsilon))^{-1}$ simulations from the proposal distribution are needed to obtain one exact sample from the target distribution, i.e. the computational complexity of the algorithm depends on the choices of $\mathbf{p}$ and $\varepsilon$. The remainder of this section will be devoted to this question.

Choice of $\mathbf{p}$ and $\varepsilon$

For a given $\varepsilon\geq 0$, the computational costs of the algorithm can be minimized by choosing $\mathbf{p}=\mathbf{p}^{*}(\varepsilon)$ such that the constant $C(\mathbf{p},\varepsilon)$ given in (3) is maximal, i.e. by choosing $\mathbf{p}$ as the solution of the nonlinear optimization problem

[TABLE]

Optimizing further w.r.t. $\varepsilon\in[0,1)$, we obtain the optimal choice $(\mathbf{p},\varepsilon)=(\mathbf{p}^{*}(\varepsilon^{*}),\varepsilon^{*})$ where $\varepsilon^{*}=\operatorname*{argmax}_{\varepsilon\in[0,1)}C(\mathbf{p}^{*}(\varepsilon),\varepsilon)$.

As the above optimization problem includes optimization steps w.r.t. $\mathbf{w}\in\mathbb{R}^{N}$, $\mathbf{p}\in\{\mathbf{x}\in[0,1]^{N}:\ \|\mathbf{x}\|_{1}=1\}$ and $\varepsilon\in[0,1]$, none of which can be solved analytically, the solution is quite involved. In order to reduce the computational burden, we simplify the problem by maximizing an analytically simpler lower bound. To this end, we decompose the convex combination $\sum_{i=1}^{N}p_{i}g_{i,\varepsilon}/f_{j}$ into sums over disjoint subsets of the form $I=\{i_{1},\ldots,i_{m}\}\subset\{1,\ldots,N\}$. For a convex combination of $(g_{i_{k},\varepsilon})_{k=1,\ldots,m}$ with weight vector $\bm{\lambda}=(\lambda_{k})_{k=1,\ldots,m}\in[0,1]^{m}$, we obtain the lower bound

[TABLE]

where we made use of the convexity of the exponential function. Setting $\kappa_{I}^{(j)}(\varepsilon,\bm{\lambda})=(1-\varepsilon)\sum_{k=1}^{m}\lambda_{k}\mathbf{C}_{\cdot i_{k}}-\mathbf{C}_{\cdot j}$, this bound can be calculated explicitly:

[TABLE]

Hence,

[TABLE]

where $\mathbf{p}_{I}=(p_{i})_{i\in I}$ for every subset $I\subset\{1,\ldots,N\}$. Now, for each $j\in\{1,\ldots,N\}$, let $J^{(j)}$ be a partition of $\{1,\ldots,N\}$, so that

[TABLE]

Thus, the RHS of (3) provides an explicit lower bound for the average acceptance probability for any choice of the $J^{(j)}$.

Remark 3.2**.**

Assume that, for some $j\in\{1,\ldots,N\}$, there is some index set $I=\{i_{1},\ldots,i_{m}\}\subset\{1,\ldots,N\}$ such that $\gamma(t_{i}-t_{j})=\Gamma$ for all $i\in I$. Then, Equation (3) provides the bound

[TABLE]

for all $\mathbf{w}\in\mathbb{R}^{N}$. Alternatively, for the same index set $I$, we could bound each summand separately, i.e.

[TABLE]

Note that, for all $\mathbf{w}\in\mathbb{R}^{N}$, the RHS of (28) is less than the RHS of (27), i.e. the lower bound is less sharp. Therefore, we prefer pooling locations with the same distance to $t_{j}$ rather than considering them separately in order to have the bound in (3) as sharp as possible.

In view of (3), instead of considering the exact value $C(\mathbf{p},\varepsilon)$ which is needed to calculate the minimal rejection rate, but cannot be given explicitly, we might maximize the function

[TABLE]

for fixed partitions $J^{(1)},\ldots,J^{(N)}$. Due to the complex dependence of $c_{I}^{(j)}$ on $\mathbf{p}$, the resulting optimization problem is nonlinear in $\mathbf{p}$ even for fixed $\varepsilon$. To circumvent this difficulty, for each $I\in J^{(j)}$, we fix $|I|$-dimensional weight vectors $\bm{\lambda}(I)$ and consider the function

[TABLE]

with $\mathbf{c}(\varepsilon;\bm{\lambda})=\{c_{ij}(\varepsilon;\bm{\lambda})\}_{1\leq i,j\leq N}$ and $c_{ij}(\varepsilon;\bm{\lambda})=c_{I}^{(j)}(\varepsilon;\bm{\lambda}(I))$ for the unique set $I\in J^{(j)}$ such that $i\in I$.

Analogously to the solution above, we first maximize $C_{\mathrm{groups}}^{\mathrm{fix}}(\cdot,\varepsilon;\bm{\lambda})$ for fixed $\varepsilon\in[0,1)$ and $\bm{\lambda}$, i.e. we consider the optimization problem

[TABLE]

To convert the linear program to standard form, we introduce an additional variable $z\in\mathbb{R}$, unconstrained in sign, leading to the equivalent program

[TABLE]

The standard form of (3) is then given by

[TABLE]

Such a linear program in standard form can be solved by standard techniques such as the simplex algorithm. Compared to the optimization of $C_{\mathrm{groups}}(\cdot,\varepsilon)$, the complementary one-dimensional problem of maximizing $C_{\mathrm{groups}}(\mathbf{p},\cdot)$ for fixed $\mathbf{p}$ can be solved rather easily.

To summarize, starting from some $\varepsilon^{*}>0$ and $\mathbf{p}^{*}=N^{-1}\mathbf{1}_{N}$, we propose to apply the following two steps repeatedly (Algorithm 2B):

Define

[TABLE]

and set

[TABLE]

i.e. the solution of the optimization problem (3) (or (3) or (3), equivalently). 2. 2.

Set $\varepsilon^{*}=\operatorname*{argmax}_{\varepsilon}C_{\mathrm{groups}}(\mathbf{p}^{*},\varepsilon)$.

Even though $C_{\mathrm{groups}}(\mathbf{p},\varepsilon)$ might be significantly smaller than $C(\mathbf{p},\varepsilon)$, in some cases, this bound is already sufficient to improve the results for $\varepsilon=0$ that have been discussed in the beginning of this section, where we have already seen that the corresponding optimal weight vector equals $\mathbf{p}^{*}=N^{-1}\mathbf{1}_{N}$ and that $C(\mathbf{p}^{*},0)=1/N$. We show an example to illustrate that this choice is not necessarily optimal, i.e. there is some $\varepsilon>0$ and a vector $\mathbf{p}$ of weights such that $C(\mathbf{p},\varepsilon)\geq C_{\mathrm{groups}}(\mathbf{p},\varepsilon)>1/N$.

Example 3.3** (Fractional Brownian Motion).**

Let $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ be $N$ equidistant locations in $[0,1]$ and $\gamma(h)=|h|^{\alpha}$ for some $\alpha>1$. Choose $\mathbf{p}=N^{-1}\mathbf{1}_{N}$, $\varepsilon=\sqrt{2}N^{-1}$ and set $J^{(1)}=\ldots=J^{(N)}=\{\{1\},\ldots,\{N\}\}$. Then, for every $x_{i}$ there are at least $\lfloor N^{-1/\alpha}\cdot N\rfloor$ locations $\mathbf{x}_{j}$ such that $\gamma(\mathbf{x}_{i}-\mathbf{x}_{j})\leq 1/N$. Thus, we obtain for large $N$ that

[TABLE]

which is eventually larger than $1/N$ as $\alpha>1$.

4 Illustration

Finally, we illustrate the performance of Algorithm 1 and the rejection sampling algorithm in an example. Taking up Example 3.3 in higher dimension, we consider the case that $Z$ is a Brown–Resnick process associated to the variogram

[TABLE]

on the grid $K=\{0,0.2,\ldots,5\}\times\{0,0.2,\ldots,5\}$ ($N=676$ points). We run four different algorithms:

1. 1A.

Algorithm 1 with proposal density $f_{\rm prop}=f$ as proposed by Oesting \BOthers. 15 2. 1B.

Algorithm 1 with proposal density $f_{\rm prop}=\sum_{i=1}^{N}p_{i}f_{i}$ where $\mathbf{p}$ is given as the solution of (2) (cf. Section 2). 3. 2A.

Algorithm 2 with proposal density $f_{\rm prop}=\frac{1}{N}\sum_{i=1}^{N}f_{i}$ and $C=1/N$ (equivalent to the procedure proposed in de Fondeville \BBA Davison 2 based on sum-normalized spectral functions) 4. 2B.

Algorithm 2 with proposal density $\tilde{f}_{\rm prop}=\sum_{i=1}^{N}p_{i}^{*}(\varepsilon^{*})g_{i,\varepsilon^{*}}$ and $C=C(\mathbf{p}^{*},\varepsilon^{*})$ in (3) where $\mathbf{p}^{*}(\varepsilon^{*})\in[0,1]^{N}$ and $\varepsilon^{*}\in(0,1)$ are obtained as described in Section 3

Even though the laws of the Brown–Resnick process $Z$ and the normalized spectral process $V^{\max}$ do not depend on the variance, but only on the variogram of the underlying Gaussian process $G$, the choice of the Gaussian process may affect the performance of the algorithms. Here, we choose the Gaussian process $G$ whose law is uniquely defined via the construction

[TABLE]

where $G_{0}$ is an arbitrary centered Gaussian process with variogram $\gamma$. Oesting \BBA Strokorb 16 show that this process has a smaller maximal variance and is thus preferable in the context of simulation.

We first calculate the optimal weights $\mathbf{p}=(p_{1},\ldots,p_{676})^{\top}$ as a solution of (2) (used in Algorithm 1B) as well as the optimal weights $\mathbf{p}^{*}(\varepsilon^{*})$ as a solution of (3) and the optimal variance modification $\varepsilon^{*}$ (used in Algorithm 2B). The results for $\mathbf{p}$ and $\mathbf{p}^{*}(\varepsilon^{*})$ are displayed in Figure 1. It can be seen that, in both cases, the weights are not spatially constant, but are larger on the boundary of the convex hull $\mathrm{conv}(K)=[0,5]\times[0,5]$ with the maximum in the corners of the square. This observation is well in line with the fact that these points have the largest contribution to $\max_{t\in K}\exp(G(t)-\mathrm{Var}(G(t))/2)$ since the variance of $G$ attains its maximum there see also 16.

We run Algorithm 1 with both proposal densities as specified above (Algorithms 1A and 1B, respectively) to obtain two Markov chains $\{W^{(k)}_{1}\}_{k=1,\ldots,1\,000\,000}$ and $\{W^{(k)}_{2}\}_{k=1,\ldots,1\,000\,000}$ of length $n_{MCMC}=1\,000\,000$. It can be seen that the empirical acceptance rate of the second chain ($0.855$) is remarkably higher than the one of the first chain ($0.656$) which already indicates stronger mixing. This impression is confirmed by analyzing the empirical autocorrelation functions of the time series $\{\|\exp(W^{(k)}_{1})\|_{\infty}\}_{k=1,\ldots,1\,000\,000}$ and $\{\|\exp(W^{(k)}_{2})\|_{\infty}\}_{k=1,\ldots,1\,000\,000}$ which are shown in Figure 2. Here, the empirical autocorrelation in the second chain is drastically reduced in comparison with the first chain, indicating that two states of the chain can be regarded as nearly uncorrelated after roughly five steps.

The rejection sampling algorithm (Algorithm 2A and Algorithm 2B) automatically generates independent realizations from the multivariate target density $f_{\max}$. Therefore, we will compare them with respect to their computational complexity. To this end, we run them to generate a sample of size $100\,000$ and count the average number of simulations of Gaussian vectors from the proposal density to generate one realization of $\mathbf{W}^{\max}$. In case of Algorithm 2A, this number is $203.1$ which is close to the theoretical expression $c_{\infty}^{-1}\cdot N$. For Algorithm 2B, the number is improved by a factor of approximately $4.4$, leading to an average number of $45.9$ Gaussian vectors to be simulated to obtain one realization from the target distribution. This improvement is well in line with the corresponding value $C(\mathbf{p}^{*},\varepsilon^{*})\approx 0.0065\approx 4.4\cdot N^{-1}$.

As the example illustrates, the two modifications we suggested may lead to significant improvements of MCMC and rejection algorithms that have been proposed so far. Here, only the modified rejection sampling algorithm ensures independence of exact samples from the target distribution. However, as the example indicates, the MCMC algorithm might be particularly attractive in practice as a thinned chain results in nearly independent samples even if the thinning rate is rather small. Note that we also tried other examples such as a Brownian sheet ($\alpha=1$). However, we found that significant improvements in the rejection sampling procedure become apparent only for $\alpha>1$, see also Example 3.3.

Acknowledgements

The authors are grateful to Kirstin Strokorb, Dimitri Schwab and Jonas Brehmer for numerous valuable comments. M. Schlather has been financially supported by Volkswagen Stiftung within the project ‘Mesoscale Weather Extremes – Theory, Spatial Modeling and Prediction (WEX-MOP)’.