An improved spectral turning-bands algorithm for simulating stationary vector Gaussian random fields

Abstract

We propose a spectral turning-bands approach for the simulation of second-order stationary vector Gaussian random fields. The approach improves existing spectral methods through coupling with importance sampling techniques. A notable insight is that one can simulate any vector random field whose direct and cross-covariance functions are continuous and absolutely integrable, provided that one knows the analytical expression of their spectral densities, without the need for these spectral densities to have a bounded support. The simulation algorithm is computationally faster than circulant-embedding techniques, lends itself to parallel computing and has a low memory storage requirement. Numerical examples with varied spatial correlation structures are presented to demonstrate the accuracy and versatility of the proposal.

Appendix: proof of equations (14), (16) and (17)

Let $Y_0$ be a scalar random field. The random field regularized by a sampling function $\omega _r(\cdot )$ being the indicator function of the ball of ${\mathbb{R}}^d$ with arbitrary radius r is defined as (Chilès and Delfiner 2012)

$$\begin{aligned} Y_{r}({\mathbf{x}}) = \int _{{\mathbb{R}}^d} Y_0({\mathbf{x}}+{\mathbf{t}})\omega _r({\mathbf{t}})\, {\mathrm{d}}{\mathbf{t}}. \end{aligned}$$

If $Y_0$ is a second-order stationary random field with covariance $C_0({\mathbf{h}})$, then $Y_{r}$ is also a second-order stationary random field, therefore one can define the following covariances that depend only on ${\mathbf{h}}$.

The covariance between $Y_0({\mathbf{x}}+{\mathbf{h}})$ and $Y_{r}({\mathbf{x}})$ is

$$\begin{aligned} C_{0r}({\mathbf{h}})= & {} \int _{{\mathbb{R}}^d} {\mathbb{E}}\{Y_0({\mathbf{x}}+{\mathbf{h}})\cdot Y_0({\mathbf{x}}+{\mathbf{t}})\}{\omega }_r({\mathbf{t}})\, {\mathrm{d}}{\mathbf{t}}\\= & {} \int _{{\mathbb{R}}^d}C_0({\mathbf{h}}-{\mathbf{t}}){\omega }_r({\mathbf{t}})\, {\mathrm{d}} {\mathbf{t}}, \end{aligned}$$

that is

$$\begin{aligned} C_{0r} = C_0*\omega _r. \end{aligned}$$

(18)

Similarly, the covariance between $Y_r(\mathbf {x}+\mathbf {h})$ and $Y_{r'}(\mathbf {x})$ is (Chilès and Delfiner 2012)

$$\begin{aligned} C_{rr'} = C_0*(\breve{\omega }_r*\omega _{r'}) \end{aligned}$$

(19)

with $\breve{\omega }({\mathbf{t}})=\omega (-{\mathbf{t}})$.

Since the Fourier transformation exchanges convolution and multiplication, Eqs. (18) and (19) in terms of Fourier transforms follow respectively, as:

$$\begin{aligned} f_{0r}= & {} f_{00}\cdot \xi _{\omega _r}\\ f_{rr^\prime}= & {} f_{00}\cdot \xi _{\breve{\omega }_r}\cdot \xi _{\omega _{r^\prime}}, \end{aligned}$$

where $f_{00}$ is the spectral density of $Y_0$, and $\xi _{\omega _r}$, $\xi _{\breve{\omega }_r}$ and $\xi _{\omega _{r'}}$ are the Fourier transforms of $\omega _r$, $\breve{\omega }_r$ and $\omega _{r'}$, respectively. As $\omega _r$ is the indicator of a ball with radius r, one has $\xi _{\omega _r} = \xi _{\breve{\omega }_r} = \xi (\cdot ,r)$ (Eq. (15)) (Gradshtein and Ryzhik 1965); likewise, $\xi _{\omega _{r'}} = \xi (\cdot ,r')$.