Simulating space-time random fields with nonseparable Gneiting-type covariance functions

Abstract

Two algorithms are proposed to simulate space-time Gaussian random fields with a covariance function belonging to an extended Gneiting class, the definition of which depends on a completely monotone function associated with the spatial structure and a conditionally negative definite function associated with the temporal structure. In both cases, the simulated random field is constructed as a weighted sum of cosine waves, with a Gaussian spatial frequency vector and a uniform phase. The difference lies in the way to handle the temporal component. The first algorithm relies on a spectral decomposition in order to simulate a temporal frequency conditional upon the spatial one, while in the second algorithm the temporal frequency is replaced by an intrinsic random field whose variogram is proportional to the conditionally negative definite function associated with the temporal structure. Both algorithms are scalable as their computational cost is proportional to the number of space-time locations that may be irregular in space and time. They are illustrated and validated through synthetic examples.

Appendices

A Proofs for the spectral approach

1.1 A.1 Proof of Theorem 1

We use the notations defined in Sect. 3. Substituting (7) in (15), one obtains

$$\begin{aligned}&{\mathbb {E}}\left[ e^{ \textstyle i \, u \, \Omega ^T(\varvec{\omega },r) } \right] \nonumber \\&\quad = \exp \bigl ( - \lambda (\varvec{\omega },r) \gamma (u) \bigr ) \nonumber \\&\quad = \exp \left( - \lambda (\varvec{\omega },r) \int _{{\mathbb {R}}}\left( 1 - \cos ( u x) \right) {\mathcal {X}}(dx)\right) \nonumber \\&\quad = \exp \left( \lambda (\varvec{\omega },r) \int _{{\mathbb {R}}}\left( e^{i u x}-1-\frac{i u x}{1+x^2}\right) {\mathcal {X}}(dx)\right) ,\nonumber \\&\quad u \in {\mathbb {R}}, \end{aligned}$$

(22)

where $\lambda (\varvec{\omega },r)=\frac{|\varvec{\omega }|^2}{4r}$ and where the last equality stems from the symmetry of ${\mathcal {X}}$ and the integrability condition (8). Therefore, the distribution $F_T ( d \tau \mid \varvec{\omega }, r)$ is an infinitely divisible distribution with Lévy measure ${\nu }_{\varvec{\omega },r}=\lambda (\varvec{\omega },r){\mathcal {X}}$.

Assume first that ${\mathcal {X}}({\mathbb {R}})=+\infty$. In this case, ${\nu }_{\varvec{\omega },r}({\mathbb {R}})=+\infty$ for $\varvec{\omega }\ne 0$ and the Lévy measure ${\nu }_{\varvec{\omega },r}$ is absolutely continuous, since ${\mathcal {X}}$ is absolutely continuous. Then, by applying Lemma 1 of Sato (1982) the infinitely divisible distribution $F_T ( d \tau \mid \varvec{\omega }, r)$ is absolutely continuous for $\varvec{\omega }\ne 0$ and for a.e. $\varvec{\omega }$. As a consequence, since $F_S ( d\varvec{\omega }\mid r)$ is also absolutely continuous, it follows that

$$\begin{aligned} F(d\varvec{\omega },d\tau )=\int _0^{+\infty } F_S ( d\varvec{\omega }\mid r)F_T (d \tau \mid \varvec{\omega }, r)\mu (dr) \end{aligned}$$

is absolutely continuous.

Let us now assume $\vartheta ={\mathcal {X}}({\mathbb {R}})<+\infty$. In this case, (22) can be rewritten as

$$\begin{aligned}&{\mathbb {E}}\left[ e^{ \textstyle i \, u \, \Omega ^T(\varvec{\omega },r) } \right] = \exp \left( \lambda (\varvec{\omega },r) \int _{{\mathbb {R}}}\left( e^{i u x}-1\right) {\mathcal {X}}(dx)\right) , \end{aligned}$$

by symmetry of the finite measure ${\mathcal {X}}$. Since $0< \vartheta < +\infty$, $\Omega ^T(\varvec{\omega },r)$ has a compound Poisson distribution with ${\mathbb {P}}\bigl (\Omega ^T(\varvec{\omega },r)=0)\bigr ) =\exp \bigl (-\vartheta \lambda (\varvec{\omega },r)\bigr )$, which implies that

$$\begin{aligned} {\mathbb {P}}\left( \Omega ^T=0\right) ={\mathbb {E}}\left[ \exp \left( -\frac{\vartheta \vert \varvec{\Omega }^S\vert ^2 }{4R}\right) \right] >0. \end{aligned}$$

In conclusion, when $\vartheta < +\infty$, the distribution of $\Omega ^T$ has an atom at 0. Hence, F is not absolutely continuous in this case.

1.2 A.2 Proof of Theorem 2

Recall that

$$\begin{aligned}&C (\varvec{h}, u \mid r ) = \frac{1}{\bigl (\gamma (u)+1 \bigr )^{k/2}} \, \exp \left( - \frac{r |\varvec{h}|^2}{\gamma (u)+1} \right) \nonumber \\&:= \int _{{\mathbb {R}}^k} \int _{\mathbb {R}}e^{ \textstyle i \langle \varvec{h}, \varvec{\omega }\rangle + i u \tau } \, F_T ( d \tau \mid \varvec{\omega }, r) \, F_S ( d \varvec{\omega }\mid r) . \end{aligned}$$

(23)

If $u=0$, then

$$\begin{aligned} C (\varvec{h}, 0 \mid r ) = \exp \left( - r |\varvec{h}|^2 \right) := \int _{{\mathbb {R}}^k} e^{ \textstyle i \langle \varvec{h}, \varvec{\omega }\rangle } \, F_S ( d \varvec{\omega }\mid r), \end{aligned}$$

which is nothing but (14). This shows that $F_S$ possesses the density

$$\begin{aligned} f_S(\varvec{\omega }\mid r ) = \frac{1}{(4\pi r)^{k/2}} \exp \left( - \frac{ |\varvec{\omega }|^2}{4 r} \right) . \end{aligned}$$

(24)

Plugging (24) into (23), one obtains

$$\begin{aligned} C (\varvec{h}, u \mid r )&= \frac{1}{(4\pi r)^{k/2}} \int _{{\mathbb {R}}^k} e^{ \textstyle i \langle \varvec{h}, \varvec{\omega }\rangle } \nonumber \\&\quad \times \exp \left( - \frac{ |\varvec{\omega }|^2}{4 r} \right) \, \int _{\mathbb {R}}e^{ \textstyle i u \tau } \, F_T ( d \tau \mid \varvec{\omega }, r) \, d \varvec{\omega }. \end{aligned}$$

(25)

On the other hand, up to a multiplicative factor, the function $\varvec{\omega }\, \mapsto \exp \left( -\, r |\varvec{h}|^2 / (\gamma (u)+1) \right)$ is the Fourier transform of a Gaussian random vector:

$$\begin{aligned} C (\varvec{h}, u \mid r )&= \frac{1}{(4\pi r)^{k/2}} \int _{{\mathbb {R}}^k} e^{ \textstyle i \langle \varvec{h}, \varvec{\omega }\rangle } \nonumber \\&\quad \times \, \exp \left( - \frac{ |\varvec{\omega }|^2 \, (\gamma (u)+1)}{4 r} \right) \, d \varvec{\omega }. \end{aligned}$$

(26)

Comparing (25) and (26), the injectivity of the Fourier transform implies

$$\begin{aligned}&\exp \left( - \frac{ |\varvec{\omega }|^2}{4 r} \right) \, \int _{\mathbb {R}}e^{ \textstyle i u \tau } \, F_T ( d \tau \mid \varvec{\omega }, r) \\&\quad = \exp \left( - \frac{ |\varvec{\omega }|^2 \, (\gamma (u)+1) }{4 r} \right) \qquad \varvec{\omega }\text {-a.e.}, \end{aligned}$$

or equivalently

$$\begin{aligned} \int _{\mathbb {R}}e^{ \textstyle i u \tau } \, F_T ( d \tau \mid \varvec{\omega }, r) = \exp \left( - \frac{ |\varvec{\omega }|^2 \gamma (u)}{4 r} \right) \qquad \varvec{\omega }\text {-a.e.}, \end{aligned}$$

which is precisely (15). $\square$

1.3 A.3 On the generic approach

The aim of this section is to show that $\sum _{n \ge 1} X_{\mathrm{T}_n}$ and $\Omega ^T (\varvec{\omega },r)$ have the same distribution. This is done by comparing their Fourier transforms. Remind that the spectral measure of the variogram is positive, symmetric, without an atom at the origin, and satisfies the integrability property

$$\begin{aligned} \int _{{\mathbb {R}}} \frac{x^2 \, {\mathcal {X}}(dx)}{1 + x^2} = A < +\infty . \end{aligned}$$

(27)

Let us start with

$$\begin{aligned} {\mathcal {X}}(dx) = \int _{{\mathbb {R}}_+} \exp \left( - t \frac{x^2}{1 + x^2} \right) \, \frac{x^2 \, {\mathcal {X}}(dx)}{1 + x^2} \, dt. \end{aligned}$$

Because of (27), the positive function $\theta$ defined on ${\mathbb {R}}_+$ by

$$\begin{aligned} \theta (t) = \int _{{\mathbb {R}}} \exp \left( - t \frac{x^2}{1 + x^2} \right) \, \frac{x^2 \, {\mathcal {X}}(dx)}{1 + x^2} \end{aligned}$$

is upper bounded by A. It follows that, for each $t>0$, the measure

$$\begin{aligned} {\mathcal {X}}_t (dx) = \frac{1}{\theta (t)} \exp \left( - t \frac{x^2}{1 + x^2} \right) \, \frac{x^2 \, {\mathcal {X}}(dx)}{1 + x^2} \end{aligned}$$

is a probability measure on ${\mathbb {R}}$. This measure is symmetric, and satisfies

$$\begin{aligned} {\mathcal {X}}(dx) = \int _{{\mathbb {R}}_+} {\mathcal {X}}_t (dx) \, \theta (t) \, dt. \end{aligned}$$

(28)

Consider now a Poisson point process $\bigl (\mathrm{T}_n: n \ge 1 \bigr )$ with intensity $\lambda (t) = \lambda \, \theta (t)$ on ${\mathbb {R}}_+$ ($\lambda$ is put here as a short notation for $\frac{|\varvec{\omega }|^2}{4r}$). Since $\theta$$\lambda (t)$ is upper bounded by $\lambda A$, this process has no accumulation point. Consider also a family $\bigl ( X_t: t \in {\mathbb {R}}_+ \bigr )$ of independent random variables, with $X_t$ being distributed as ${\mathcal {X}}_t$. Because ${\mathcal {X}}_t$ is symmetric, the Fourier transform of $X_t$ can be written as

$$\begin{aligned} {\mathbb {E}}\bigl [ \exp ( i u X_t ) \bigr ] = \int _{{\mathbb {R}}} \cos (u x) \, {\mathcal {X}}_t (dx). \end{aligned}$$

(29)

In what follows, we calculate the Fourier transform of $\mathrm{T}= \sum _{n \ge 1} X_{\mathrm{T}_n}$. Denoting by $\Lambda (t_0)$ the integral of $\lambda (t)$ on $]0,t_0[$, we have

$$\begin{aligned}&{\mathbb {E}}\bigl [ \exp ( i u \mathrm{T}) \bigr ]\\&\quad = \lim _{t_0 \longrightarrow +\infty } \sum _{n=0}^{+\infty } \exp \bigl ( - \Lambda (t_0) \bigr ) \, \frac{\Lambda ^n (t_0)}{n!} \\&\qquad \times \left[ \int _0^{t_0} \frac{\lambda (t)}{\Lambda (t_0)} \, {\mathbb {E}}\bigl [ \exp ( i u X_t ) \bigr ] \, dt\right] ^n \\&\quad = \lim _{t_0 \longrightarrow +\infty } \exp \left( \int _0^{t_0} {\mathbb {E}}\bigl [ \exp ( i u X_t ) - 1 \bigr ] \, \lambda (t)\, dt \right) \\&\quad = \exp \left( \int _0^{+\infty } {\mathbb {E}}\bigl [ \exp ( i u X_t ) - 1 \bigr ] \, \lambda (t) \, dt \right) . \end{aligned}$$

This implies, owing to (29)

$$\begin{aligned} {\mathbb {E}}\bigl [ \exp ( i u \mathrm{T}) \bigr ] = \exp \left( \int _0^{+\infty } \int _{{\mathbb {R}}} \bigl [ \cos (u x) - 1 \bigr ] \, {\mathcal {X}}_t (dx) \, \lambda (t) \, dt \right) . \end{aligned}$$

Permuting the integrals and replacing $\lambda (t)$ by its expression, we obtain

$$\begin{aligned} {\mathbb {E}}\bigl [ \exp ( i u \mathrm{T}) \bigr ] =\exp \left( \frac{ \vert \varvec{\omega }|^2}{4 r} \, \int _{{\mathbb {R}}} \bigl [ \cos (u x) - 1 \bigr ] \, {\mathcal {X}}(dx) \right) . \end{aligned}$$

Finally, the spectral representation (7) of $\gamma$ gives

$$\begin{aligned} {\mathbb {E}}\bigl [ \exp ( i u \mathrm{T}) \bigr ] =\exp \left( - \frac{ \vert \varvec{\omega }|^2}{4 r} \gamma (u) \right) , \end{aligned}$$

which is precisely the Fourier transform (15) of $\Omega ^T (\varvec{\omega },r)$. $\Box$

1.4 A.4 Implementing the generic approach for logarithmic variograms

The construction of $\theta$ and ${\mathcal {X}}_t$ proposed in Appendix A.3 is not necessarily unique. Starting from

$$\begin{aligned} {\mathcal {X}}(d x)= & {} \frac{\exp (- a |x |)}{ |x |\, \ln a^2} \\= & {} \frac{1}{\ln a^2} \int _a^{+\infty } \exp \bigl (- t \, |x |\bigr ) \, d t , \end{aligned}$$

it appears that a possible decomposition such as (16) can be obtained by taking

$$\begin{aligned} \theta (t) = \frac{2}{t \, \ln a^2} \, 1_{t \ge a} =\frac{1}{t \, \ln a} \, 1_{t \ge a} \end{aligned}$$

and

$$\begin{aligned} {\mathcal {X}}_t (dx) = \frac{1}{2} t \, e^{ - t |x \vert } \, dx, \qquad t > a . \end{aligned}$$

Consider now a Poisson point process $\bigl ( \mathrm{T}_n : \, n \ge 1 \bigr )$ on ${\mathbb {R}}_+$ with intensity function

$$\begin{aligned} \lambda (t) = \frac{ |\varvec{\omega }|^2}{4 \, r } \, \theta (t) := \frac{\lambda }{t} \, 1_{t \ge a}. \end{aligned}$$

A simple approach to simulate a nonhomogeneous Poisson point process is recommended by Devroye (1986, Chapter 6). It consists in calculating the inverse of the primitive of the intensity function $\lambda (t)$ that vanishes at a, i.e. $\Lambda (t) = \lambda \, \ln ( t / a)$, at the points of a homogeneous Poisson process with a unit intensity on the positive half-line, as shown in Fig. 3.

In this figure, the $U_i$’s are independent standard uniform variables and are related to the Poisson time $\mathrm{T}_n$ by the formula $\Lambda ( \mathrm{T}_n ) = \lambda \, \ln ( \mathrm{T}_n / a ) = - \ln \bigl ( U_1 \cdots U_n \bigr )$, which gives

$$\begin{aligned} \mathrm{T}_n = \frac{a}{(U_1 \cdots U_n )^{1 / \lambda }} . \end{aligned}$$

(30)

Fig. 3

Simulation of a nonhomogeneous Poisson point process by inversion of a homogeneous Poisson process with a unit intensity

Now, recall that $\Omega ^T(\varvec{\omega },r)$ has the same distribution as $\sum _{n=1}^{+\infty } X_{\mathrm{T}_n}$, where each $X_t$ is distributed as ${\mathcal {X}}_t$. Because the $X_t$’s are independent, we have

$$\begin{aligned} \mathrm {Var} \Bigl [ \sum _{n \ge 1} X_{\mathrm{T}_n} \Bigr ] =\sum _{n \ge 1} \mathrm {Var} \bigl [ X_{\mathrm{T}_n} \bigr ]. \end{aligned}$$

Moreover, (30) implies

$$\begin{aligned} \mathrm {Var} \bigl [ X_{\mathrm{T}_n} \bigr ]= & {} {\mathbb {E}}\Bigl [ \mathrm {Var} \bigl [ X_{\mathrm{T}_n} \vert \mathrm{T}_n \bigr ] \Bigr ] = {\mathbb {E}}\Bigl [ \frac{2}{\mathrm{T}^2_n} \Bigr ]\\= & {} \frac{2}{a^2} \, {\mathbb {E}}\Bigl [ (U_1 \cdots U_n)^{2 \lambda } \Bigr ] =\frac{2}{a^2} \, \left( \frac{\lambda }{\lambda + 2}\right) ^n. \end{aligned}$$

Consequently

$$\begin{aligned} \mathrm {Var} \Bigl [ \sum _{n \ge 1} X_{\mathrm{T}_n} \Bigr ] =\frac{2}{a^2} \, \sum _{n \ge 1} \left( \frac{\lambda }{\lambda + 2}\right) ^n = \frac{\lambda }{a^2}. \end{aligned}$$

Similarly, if the series is truncated at order $n_0$, then the same calculation leads to the residual variance

$$\begin{aligned} \mathrm {Var} \Bigl [ \sum _{n \ge n_0+1} X_{\mathrm{T}_n} \Bigr ]= & {} \frac{2}{a^2} \sum _{n \ge n_0+1} \left( \frac{\lambda }{\lambda + 2}\right) ^n \\= & {} \left( \frac{\lambda }{\lambda + 2} \right) ^{n_0} \frac{\lambda }{a^2}. \end{aligned}$$

Let $\varepsilon > 0$ be arbitrarily small. From the previous calculations, it follows that

$$\begin{aligned} \frac{ \mathrm {Var} \Bigl [ \sum _{n \ge n_0+1} X_{\mathrm{T}_n} \Bigr ]}{\mathrm {Var} \Bigl [ \sum _{n \ge 1} X_{\mathrm{T}_n} \Bigr ] }< \varepsilon \\Longleftrightarrow & {} \ \left( \frac{\lambda }{\lambda + 2} \right) ^{n_0} < \varepsilon \\\Longleftrightarrow & {} \ n_0 > \frac{- \ln \epsilon }{\ln ( 1 + 2 / \lambda ) }. \end{aligned}$$

$\square$

B Proofs for the substitution approach

1.1 Proof of Theorem 3

For $(\varvec{x},t) \in {\mathbb {R}}^k \times {\mathbb {R}}$, $Z(\varvec{x},t)$ conditional on $(R,\tilde{\varvec{\Omega }},U,W)$ (i.e. only letting $\Phi$ vary randomly) has a zero expectation, insofar as it is proportional to the cosine of a random variable uniformly distributed on an interval of length $2 \pi$. The prior expectation of $Z(\varvec{x},t)$ is therefore zero:

$$\begin{aligned} {\mathbb {E}}[Z(\varvec{x},\varvec{t}) ] = {\mathbb {E}}[ {\mathbb {E}}[ Z(\varvec{x},t) \mid R, \tilde{\varvec{\Omega }}, U, W ] ] = 0. \end{aligned}$$

Let us now calculate the covariance between the random variables $Z(\varvec{x},t)$ and $Z(\varvec{x}^{\prime },t^{\prime })$, with $(\varvec{x},t) \in {\mathbb {R}}^k \times {\mathbb {R}}$ and $(\varvec{x}^{\prime },t^{\prime }) \in {\mathbb {R}}^k \times {\mathbb {R}}$:

$$\begin{aligned}&{\mathbb {E}}[ Z(\varvec{x},t) Z(\varvec{x}^{\prime },t^{\prime })] \\&\quad = 2 {\mathbb {E}}\bigg [-\ln (U) \cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}\rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} W(t) + \Phi \right) \\&\qquad \times \cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}^{\prime } \rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} W(t^{\prime }) + \Phi \right) \bigg ] \\&\quad = 2 {\mathbb {E}}\bigg [\cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }},\varvec{x}\rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} W(t) + \Phi \right) \\&\qquad \times \cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}^{\prime } \rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} W(t^{\prime }) + \Phi \right) \bigg ]. \end{aligned}$$

The last equality stems from the fact that $-\ln (U)$ is an exponential random variable with mean 1 and is independent of $(R,\tilde{\varvec{\Omega }},W)$. Using the product-to-sum trigonometric identities, one can write the product of cosines as half the sum of two cosines, namely:

• the cosine of the difference: $\cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}-\varvec{x}^{\prime } \rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} (W(t)-W(t^{\prime })) \right)$

• the cosine of the sum: $\cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}+\varvec{x}^{\prime } \rangle + \frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} (W(t)+W(t^{\prime })) + 2 \Phi \right)$.

However, because $\Phi$ is uniformly distributed on $(0, 2\pi )$ and independent of $(R, \tilde{\varvec{\Omega }},W)$, the expectation of the cosine of the sum is 0. It remains

$$\begin{aligned}&{\mathbb {E}}[Z(\varvec{x},t) Z(\varvec{x}^{\prime },t^{\prime })]\\&\quad = {\mathbb {E}}\bigg [\cos \left( \sqrt{2R} \, \langle \tilde{\varvec{\Omega }}, \varvec{x}-\varvec{x}^{\prime } \rangle +\frac{|\tilde{\varvec{\Omega }}|}{\sqrt{2}} (W(t)-W(t^{\prime })) \right) \bigg ]. \end{aligned}$$

The increment $W(t)-W(t^{\prime })$ is a Gaussian random variable with zero mean and variance $2 \gamma (t-t^{\prime })$ and is independent of $(R,\tilde{\varvec{\Omega }})$, i.e.:

$$\begin{aligned} W(t)-W(t^{\prime }) = \sqrt{2 \gamma (t-t^{\prime })} Y, \end{aligned}$$

with $Y \sim {\mathcal {N}} (0,1)$ independent of $(R,\tilde{\varvec{\Omega }})$. Defining $\varvec{h}= \varvec{x}- \varvec{x}^{\prime }$ and $u=t - t^{\prime }$ and denoting by g the standard Gaussian probability density, one therefore obtains:

$$\begin{aligned}&{\mathbb {E}}[ Z(\varvec{x},t) Z(\varvec{x}^{\prime },t^{\prime })] \nonumber \\&\quad = \int _{{\mathbb {R}}_+} \int _{{\mathbb {R}}^k} \int _{{\mathbb {R}}} \cos \left( \sqrt{2r} \, \langle \tilde{\varvec{\omega }}, \varvec{h}\rangle + |\tilde{\varvec{\omega }} |\sqrt{\gamma (u)} y \right) \nonumber \\&\qquad \times g(y) dy \frac{1}{(2\pi )^{k/2}} \, \exp \left( -\frac{|\tilde{\varvec{\omega }} |^2}{2} \right) d \tilde{\varvec{\omega }} \mu (dr) \nonumber \\&\quad = \frac{1}{(2\pi )^{k/2}} \int _{{\mathbb {R}}_+} \int _{{\mathbb {R}}^k} \cos \left( \sqrt{2r} \, \langle \tilde{\varvec{\omega }}, \varvec{h}\rangle \right) \nonumber \\&\qquad \times \int _{{\mathbb {R}}} \cos \left( |\tilde{\varvec{\omega }}|\sqrt{\gamma (u)} y \right) g(y) dy \, \exp \left( -\frac{|\tilde{\varvec{\omega }} |^2}{2} \right) d\tilde{\varvec{\omega }} \mu (dr) . \end{aligned}$$

(31)

The last equality in (31) stems from the angle-sum trigonometric identity and the fact that $\langle \tilde{\varvec{\omega }}, \varvec{h}\rangle$ is an odd function of $\tilde{\varvec{\omega }}$ and $|\tilde{\varvec{\omega }}|\sqrt{\gamma (u)} y$ is an even function of $\tilde{\varvec{\omega }}$.

The simulated random field Z is therefore second-order stationary, since its expectation is identically zero and the covariance between any two variables $Z(\varvec{x},t)$ and $Z(\varvec{x}^{\prime },t^{\prime })$ only depends on $\varvec{h}= \varvec{x}- \varvec{x}^{\prime }$ and $u = t - t^{\prime }$. Up to a multiplicative factor, the last integral in (31) appears as the Fourier transform of the standard Gaussian probability density g(y) on ${\mathbb {R}}$. Specifically:

$$\begin{aligned} \int _{{\mathbb {R}}} \cos \left( |\tilde{\varvec{\omega }}|\sqrt{\gamma (u)} y \right) g(y) dy = \exp \left( -|\tilde{\varvec{\omega }}|^2 \frac{\gamma (u)}{2} \right) . \end{aligned}$$

Hence:

$$\begin{aligned}&{\mathbb {E}}[ Z(\varvec{x},t) Z(\varvec{x}^{\prime },t^{\prime })]\nonumber \\&\quad = \frac{1}{(2\pi )^{k/2}} \int _{{\mathbb {R}}_+} \int _{{\mathbb {R}}^k} \cos (\sqrt{2r} \, \langle \tilde{\varvec{\omega }},\varvec{h}\rangle ) \nonumber \\&\qquad \times \exp \left( - |\tilde{\varvec{\omega }}|^2 \frac{\gamma (u)+1}{2} \right) d\tilde{\varvec{\omega }} \mu (dr) \nonumber \\&\quad = C(\varvec{h},u). \end{aligned}$$

(32)

The last equality in (32) stems from (11) and completes the proof. $\square$