Vol.:(0123456789) https://doi.org/10.1007/s11009-022-09958-x 1 3 On Some Distributional Properties of Subordinated Gaussian Random Fields Robin Merkle1  · Andrea Barth1 Received: 28 May 2021 / Revised: 3 March 2022 / Accepted: 15 March 2022 © The Author(s) 2022 Abstract Motivated by the subordinated Brownian motion, we define a new class of (in general dis- continuous) random fields on higher-dimensional parameter domains: the subordinated Gaussian random field. We investigate the pointwise marginal distribution of the con- structed random fields, derive a Lévy-Khinchin-type formula and semi-explicit formulas for the covariance function. Further, we study the pointwise stochastic regularity and pre- sent various numerical examples. Keywords Subordinated Gaussian random fields · Lévy fields · Pointwise distribution · Pointwise stochastic regularity Mathematics Subject Classification (2010) 60G60 · 60G51 · 60G55 · 60G15 · 60E05 1 Introduction In many applications of stochastic modeling, it is meaningful to consider random fields which are discontinous in space (e.g. in fractured porous media modeling). In the situation of a one-dimensional parameter space, like financial modeling, Lévy processes turned out to be a very powerful class of (in general) discontinuous stochastic processes, combined with useful properties, see for example Schoutens (2003), Applebaum (2009), Sato (2013). Whereas the extension of ℝ-valued Lévy processes with one-dimensional parameter space to Hilbert space H -valued Lévy processes is straight forward (see for example Barth and Stein (2018b)), the extension of Lévy processes to higher-dimensional param- eter spaces is more challenging. The reason can be found at the very starting point of the definition of Lévy processes where time increments are considered: In fact, the definition of Lévy processes makes explicitly use of the total ordered structure underlying the con- sidered time interval. The absence of such a structure on a higher-dimensional parameter * Robin Merkle robin.merkle@mathematik.uni-stuttgart.de Andrea Barth andrea.barth@mathematik.uni-stuttgart.de 1 IANS SimTech, University of Stuttgart, Stuttgart, Germany Published online: 6 April 2022 Methodology and Computing in Applied Probability (2022) 24:2661–2688 / http://orcid.org/0000-0002-8120-689X http://crossmark.crossref.org/dialog/?doi=10.1007/s11009-022-09958-x&domain=pdf 1 3 space makes it difficult to extend the definition of a standard Lévy process to higher-dimen- sional parameter spaces. Subordinated fields did receive only little attention in the recent literature. In some classical papers on generalized random fields, of which Dobrushin (1979) is an important representative (see also the references therein), subordinated fields are defined in terms of iterated Itô-integrals. In the recent article, Makogin and Spodarev (2021), the authors investigate deterministic transformations of Gaussian random fields, so called Gaussian subordinated fields, and study excursion sets. The Rosenblatt distributions and long-range dependence of (subordinated) fields are looked into in Leonenko et  al. (2017). The arti- cle Barndorff-Nielsen et al. (2001) presents an extension of the concept of subordination to multivariate Lévy processes and investigates self-decomposibility of the resulting pro- cesses defined on a one-dimensional parameter domain. In the recent paper Buchmann et al. (2019), the authors define the so called weak subordination of multivariate Lévy pro- cesses as a generalization of the classical subordination. The resulting multivariate pro- cesses depend on a one-dimensional time parameter and the authors prove that weak subor- dination is an extension of the classical (strong) subordination. In Buchmann et al. (2017) and Buchmann et al. (2020), the authors consider multivaritate Brownian motions (weakly) subordinated by multivariate Thorin subordinators, investigate self-decomposibility as well as the existence of moments of the resulting distributions and present some applications in mathematical finance. The considered subordinated multivariate processes are defined on a one-dimensional parameter domain. In contrast, the main contribution of our work is to prove properties of the (discontinu- ous) subordinated random fields on higher-dimensional parameter domains and of their pointwise distributions, which are important in applications (see for example Zhang and Kang (2004), Bastian (2014) and Barth and Stein (2018a)). We present an approach for an extension of a subclass of Lévy processes to more gen- eral parameter spaces: Motivated by the subordinated Brownian motion, we employ a higher-dimensional subordination approach using a Gaussian random field together with Lévy subordinators. Figure  1 illustrates the approach with samples of a Gaussian random field (GRF) on [0, 1]2 with Matérn-1.5 covariance function and the corresponding subordinated field, where we used Poisson and Gamma processes on [0, 1] to subordinate the GRF. These examples illustrate how the jumps of the Lévy subordinators produce jumps in the two-dimensional subordinated GRF. The flexibility of the resulting random fields make them attractive for a variety of applications. In a recent article by Barth and Merkle (2020), the authors consider a randomized elliptic partial differential equation, where the subordi- nated GRF occur in the diffusion coefficient, to name just one possible application. Fig. 1 Sample of Matérn-1.5-GRF (left), Poisson-subordinated GRF (middle) and Gamma-subordinated GRF (right) 2662 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 The question arises whether it is possible to transfer some theoretical results of one- dimensional Lévy processes to these random fields on higher-dimensional parameter spaces. In particular, a Lévy-Khinchin-type formula to access the pointwise distribution of the constructed random field is of great interest (see Sect. 4). In Sect. 5 we investigate the covariance structure of the subordinated fields and show how it is influenced by the choice of subordinators. The stochastic regularity of the subordinated fields is studied in Sect. 6. There, we derive conditions which ensure the existence of pointwise moments. In the last section we present some numerical experiments on the theoretical results presented in this paper intended to help fitting random fields to data. 2 Preliminaries In this section we give a short introduction to Lévy processes and Gaussian random fields as basis for the construction of subordinated Gaussian random fields. Throughout the paper, let (Ω,F,ℙ) be a complete probability space. 2.1 Lévy Processes Let T ⊆ ℝ+ ∶= [0,+∞) be an arbitrary time domain. A stochastic process X = (X(t), t ∈ T) on T is a family of random variables on the probability space (Ω,F,ℙ) . A stochastic process l on T = [0,+∞) is said to be a Lévy process if l(0) = 0 ℙ − a.s. , l has independent and stationary increments and l is stochastically continuous. A very impor- tant characterization property of Lévy processes is given by the so called Lévy-Khinchin formula. Theorem  2.1 (Lévy-Khinchin formula, see (Applebaum  2009, Th. 1.3.3 and p. 29)) Let l be a real-valued Lévy process on T = ℝ+ ∶= [0,+∞) . There exist parameters b ∈ ℝ, �2 N ∈ ℝ+ and a measure � on (ℝ,B(ℝ)) such that the pointwise characteristic func- tion �l(t) , for t ∈ T , admits the representation for � ∈ ℝ . The measure � on (ℝ,B(ℝ)) satisfies and a measure with this property is called Lévy measure. It follows from the Lévy-Khinchin formula that every Lévy process is fully character- ized by the so called Lévy triplet (b, �2 N , �). Within the class of Lévy processes there exists a subclass which is given by the so called subordinators: A Lévy subordinator on T is a Lévy process that is non-decreasing ℙ -almost surely. The pointwise characteristic function of a Lévy subordinator l(t), for t ∈ T , admits the form �l(t)(�) ∶= 𝔼(exp(i�l(t))) = exp ( t ( ib� − �2 N 2 �2 + � ℝ⧵{0} ei�y − 1 − i�y1{|y|≤1}(y) �(dy) )) , ∫ ℝ min(y2, 1) 𝜈(dy) < ∞, 2663Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 (see (Applebaum 2009, Theorem 1.3.15)). Here, � is the Lévy measure and � is called drift parameter of l. The Lévy measure � on (ℝ,B(ℝ)) of a Lévy subordinator satisfies Since any Lévy subordinator l is a Lévy process, the Lévy-Khinchin formula holds and we obtain �2 N = 0 and b = � + ∫ 1 0 y �(dy) in Theorem 2.1 for the subordinator l. In the follow- ing, we always mean the triplet (� , 0, �) corresponding to representation (1) if we refer to the characteristic triplet of a Lévy subordinator. 2.2 Gaussian Random Fields Let D ⊂ ℝ d be a spatial domain. A random field R = (R(x), x ∈ D) is a family of random variables on the probability space (Ω,F,ℙ) . In our approach to extend Lévy processes on higher-dimensional parameter domains, one important component is given by the Gaussian random field. Definition 2.2 (see (Adler and Taylor 2007, Sc. 1.2)) A random field W ∶ Ω ×D → ℝ on a d-dimensional domain D ⊂ ℝ d is said to be a Gaussian random field (GRF) if, for any x(1),… , x(n) ∈ D with n ∈ ℕ , the n-dimensional random variable (W(x(1)),… ,W(x(n))) is multivariate Gaussian distributed. For a GRF W and arbitrary points x(1), x(2) ∈ D , we define the mean function by �W (x (1)) ∶= �(W(x(1))) and the covariance function by The GRF W is called centered, if �W (x (1)) = 0 for all x(1) ∈ D. Note that every Gaussian random field is determined uniquely by its mean and covariance function. We denote by Q ∶ L 2(D) → L 2(D) the covariance operator of W which is, for � ∈ L 2(D) , defined by Here, L2(D) denotes the set of all square integrable functions over D . Further, if D is compact, there exists a decreasing sequence (�i, i ∈ ℕ) of real eigenvalues of Q with cor- responding eigenfunctions (ei, i ∈ ℕ) ⊂ L 2(D) which form an orthonormal basis of L2(D) (see (Adler and Taylor 2007, Section 3.2) and (Werner 2011, Theorem VI.3.2 and Chap- ter II.3)). The GRF W is called stationary if the mean function �W is constant and the covariance function qW (x(1), x(2)) only depends on the difference x(1) − x(2) of the values x(1), x(2) ∈ D . Further, the stationary GRF W is called isotropic if the covariance function qW (x (1), x(2)) only depends on the Euclidean length ‖x(1) − x(2)‖2 of the difference of the val- ues x(1), x(2) ∈ D (see Adler and Taylor (2007), p. 102 and p. 115). (1)�l(t)(�) = 𝔼(exp(i�l(t))) = exp ( t ( i�� + ∫ ∞ 0 ei�y − 1 �(dy) )) , for � ∈ ℝ, 𝜈(−∞, 0) = 0 and ∫ ∞ 0 min(y, 1) 𝜈(dy) < ∞. qW (x (1), x(2)) ∶= � ( (W(x(1)) − �W (x (1)))(W(x(2)) − �W (x (2))) ) . Q(�)(x) = ∫ D qW (x, y)�(y)dy, for x ∈ D. 2664 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 3 The Subordinated Gaussian Random Field Throughout the rest of this paper, let d ∈ ℕ be a natural number with d ≥ 2 and T1,… , Td > 0 be positive values. We define the horizon vector � ∶= (T1,… , Td) and consider the spatial domain [0, 𝕋 ]d ∶= [0, T1] ×⋯ × [0,Td] ⊂ ℝ d . In the following, we will also use the notation (0, � ]d ∶= (0, T1] ×⋯ × (0,Td] . After a short motivation we define next the subordinated field and show that it is indeed measurable. 3.1 Motivation: The Subordinated Brownian Motion In order to motivate the novel subordination approach for the extension of Lévy pro- cesses, we shortly repeat the main ideas of the subordinated Brownian motion which is defined as a Lévy-time-changed Brownian motion: Let B = (B(t), t ∈ ℝ+) be a Brown- ian motion and l = (l(t), t ∈ ℝ+) be a subordinator. The subordinated Brownian motion is then defined to be the process It follows from (Applebaum 2009, Theorem 1.3.25) that the process L is again a Lévy pro- cess. Note that the class of subordinated Brownian motions is a rich class of processes with great distributional flexibility. For example, the well known Generalized Hyperbolic Lévy process can be represented as a NIG-subordinated Brownian motion (see Barth and Stein (2018b) and especially Lemma 4.1 therein). 3.2 The Definition of the Subordinated Gaussian Random Field Let W = (W(x), x = (x1,… , xd) ∈ ℝ d + ) be a GRF such that W is F⊗B(ℝd + ) −B(ℝ) -measurable. We denote by �W ∶ ℝ d + → ℝ the mean function and by qW ∶ ℝ d + ×ℝ d + → ℝ the covariance function of W. Let lk = (lk(x), x ∈ [0, Tk]) be independent Lévy subor- dinators with triplets (�k, 0, �k) , for k ∈ {1,… , d} , corresponding to representation (1). Further, we assume that the Lévy subordinators are stochastically independent of the GRF W. We consider the random field and call it subordinated Gaussian random field (subordinated GRF). Remark 3.1 Note that assuming that W has continuous paths is sufficient to ensure that W is a jointly measurable function since W is a Carathéodory function in this case (see (Aliprantis and Border 2006, Lemma 4.51)). A sufficient condition for the pathwise continuity of GRFs is given, for example, by (Adler and Taylor 2007, Theorem 1.4.1) (see also the discussion in (Adler and Taylor 2007, Section 1.3, p. 13)). A specific example for a class of GRFs with at least continuous samples is given by the Matérn GRFs: for a given smoothness parameter 𝜈 > 1 2 , a correlation parameter r > 0 and a variance parameter 𝜎2 > 0 the Matérn-� covari- ance function on ℝd + ×ℝ d + is given by qM W (x, y) = �M(‖x − y‖2) with L(t) ∶= B(l(t)), t ∈ ℝ+. L ∶ Ω × [0, 𝕋 ]d → ℝ with L(x1,… , xd) ∶= W(l1(x1),… , ld(xd)), for (x1,… , xd) ∈ [0, 𝕋 ]d, 2665Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 where Γ(⋅) is the Gamma function and K�(⋅) is the modified Bessel function of the second kind (see (Graham et  al.  2015,  Section  2.2 and Proposition 1)). Here, ‖ ⋅ ‖2 denotes the Euclidean norm on ℝd . A Matérn-� GRF is a centered GRF with covariance function qM W . The subordinated GRF constructed above is one possible way to extend the concept of the subordinated Brownian motion to higher-dimensional parameter domains. How- ever, construction of a random field by subordination in each spatial variable is not con- fined to this approach. For example, the construction itself is not limited to the case that W is a GRF and l is a Lévy subordinator. One could consider more general random fields (R(x), x ∈ ℝ d + ) subordinated by d positive valued stochastic processes. However, in general it might be difficult or impossible to investigate theoretical properties of the resulting ran- dom field. In contrast, the subordinated GRFs inherits several properties from the GRF and the Lévy subordinators investigated in the following sections. 3.3 Measurability In Subsection 3.2 we introduced the subordinated GRF L as a random field. Strictly speak- ing, we therefore have to verify that point evaluations of the field L are random variables, meaning that we have to ensure measurability of these objects. Note that this is not triv- ial, since - due to the construction of L - the Lévy subordinators induce an additional � -dependence in the spatial direction of the GRF W. The following lemma proves joint measurability of L. Lemma 3.2 Let L be a subordinated GRF on the spatial domain [0, � ]d as constructed in Subsection 3.2, where we use the notation x = (x1,… , xd) ∈ [0, � ]d . The mapping is F⊗B([0, 𝕋 ]d) −B(ℝ)-measurable. Proof For any k ∈ {1,… , d} , the Lévy process lk has càdlàg paths and, hence, the mapping lk ∶ Ω × [0,Tk] → ℝ+, is F⊗B([0, Tk]) −B(ℝ+)-measurable (see (Protter 2004, Chapter 1, Theorem 30) and (Sato 2013, Section 30)). We consider domain-extended versions of the processes: for any k ∈ {1… , d} , we define the mapping l̃k(𝜔, x) ∶= lk(𝜔, xk), for (𝜔, x) ∈ Ω × [0, � ]d , which is F⊗B([0, 𝕋 ]d) −B(ℝ+) measurable by (Aliprantis and Border (2006), Lemma 4.51). An application of (Aliprantis and Border 2006, Lemma 4.49) yields the F⊗B([0, 𝕋d]) −B(ℝd + ) -measurability of the mapping Further, the mapping (�, x) ↦ � is F⊗B([0, � ]d) −F-measurable and, hence, (Aliprantis and Border 2006, Lemma 4.49) yields that the mapping �M(s) = �2 2 1−� Γ(�) �2s√� r �� K� �2s√� r � , for s ≥ 0, L ∶ Ω × [0, 𝕋 ]d → ℝ, (�, x) ↦ W(�, l1(�, x1)… , ld(�, xd)), Ω × [0, 𝕋 ]d → ℝ d + , (𝜔, x) ↦ (l̃1(𝜔, x),… , l̃d(𝜔, x)) = (l1(𝜔, x1),… , ld(𝜔, xd)). Ω × [0, 𝕋 ]d → Ω ×ℝ d + , (𝜔, x) ↦ (𝜔, (l̃1(𝜔, x),… , l̃d(𝜔, x))) = (𝜔, (l1(𝜔, x1),… , ld(𝜔, xd))), 2666 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 is F⊗B([0, 𝕋d]) −F⊗B(ℝd + )-measurable. By assumption, the GRF W is F⊗B(ℝd + ) −B(ℝ)-measurable and, therefore, the mapping is F⊗B([0, 𝕋d]) −B(ℝ)-measurable as composition of measurable functions. ◻ 4 The Pointwise Distribution of the Subordinated GRF and the Lévy‑Khinchin‑formula In this section we prove a Lévy-Khinchin-type formula for the subordinated GRF in order to have access to the pointwise distribution. This is important, for example, in view of statistical fitting and other applications. In order to be able to do so we need the follow- ing technical lemma about the expectation of the composition of independent random variables, which is a generalization of the corresponding assertion given in the proof of (Sato 2013, Theorem 30.1). Lemma 4.1 Let W ∶ Ω ×ℝ d + → ℝ be a ℙ − a.s. continuous random field and let Z ∶ Ω → ℝ d + be a ℝd + -valued random variable which is independent of the random field W. Further, let g ∶ ℝ → ℝ be a deterministic, continuous function. It holds where m(z) ∶= �(g(W(z)) for deterministic z ∈ ℝ d + . Proof Step 1: Assume that g is globally bounded. We denote by Cb(ℝ d + ) the space of con- tinuous, bounded functions on ℝd + equipped with the usual supremum norm. We define the function which is continuous and, hence, Borel-measurable. For a fixed threshold A > 0 , we define the cut function �A(x) ∶= min(x,A) , for x ∈ ℝ and consider the random field WA(�, x) ∶= W(�,�A(x1),… ,�A(xd)) , for � ∈ Ω and x = (x1,… , xd) ∈ ℝ d + . Since W has continuous paths and [0,A]d is compact, WA has paths in Cb(ℝ d + ) and we have the pathwise identity g(WA(Z)) = F(WA, Z) . Using the independence of W and Z together with (Da Prato and Zabczyk 2014, Proposition 1.12) yields where mA(z) ∶= �(g(WA(z))) for z ∈ ℝ d + . Further, since g is continuous and bounded and W has continuous paths, we obtain the pathwise convergence g(WA(Z)) → g(W(Z)) and mA(Z) → m(Z) , for A → ∞ . Using again the boundedness of g and the dominated conver- gence theorem, we obtain �(g(W(Z))) = �(m(Z)). Step 2: In this step we assume that g(x) ≥ 0 on ℝ but g does not necessarily have to be bounded. It follows that m is also non-negative on ℝd + . Since g and m are non-negative we obtain the ℙ − a.s. monotone convergence of �A(g(W(Z)))) → g(W(Z))) for A → +∞ . We define mA(z) ∶= �(�A(g(W(z))) , for z ∈ ℝ d + , and obtain by the monotone convergence theorem L ∶ Ω × [0, 𝕋 ]d → ℝ, (�, x) ↦ W(�, (l1(�, x1)… , ld(�, xd))), �(g(W(Z)) = �(m(Z)), F ∶ Cb(ℝ d + ) ×ℝ d + → ℝ, (f , x) ↦ g(f (x)), �(g(WA(Z)) = �(F(WA, Z)) = �(�(F(WA, Z) | �(Z))) = �(mA(Z)), 2667Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Using Step 1 and the monotone convergence theorem we obtain: Step 3: Finally, we consider an arbitrary continuous function g ∶ ℝ → ℝ . We write g+ = max{0, g}, g− = −min{0, g} as well as m̃+(z) = �(g+(W(z))), m̃−(z) = �(g−(W(z))) for z ∈ ℝ d + and obtain the additive decomposition g(x) = g+(x) − g−(x) for x ∈ ℝ and m(z) = m̃+(z) − m̃−(z) for z ∈ ℝ d + by the additivity of the integral with respect to the inte- gration domain. We apply Step 2 to optain which proves the assertion. ◻ Remark 4.2 We emphasize that the assumptions on the random field W and the random variable Z in Lemma 4.1 are very mild. In particular, we do not assume the existence of continuous densities of the random field W or the random vector Z. Further, we mention that the assertion obviously may be extended to deterministic, bounded and continuous functions g ∶ ℝ → ℂ which are complex-valued. A GRF W is pointwise normally distributed with parameters specified by the mean �W and covariance function qW . Using this together with Lemma 4.1 and Remark 4.2 we obtain the following semi-explicit formula for the pointwise characteristic function of a subordinated GRF. Corollary 4.3 Let W be a ℙ − a.s. continuous GRF on ℝd + with mean function �W ∶ ℝ d + → ℝ and covariance function qW ∶ ℝ d + ×ℝ d + → ℝ . Further, let lk = (lk(t), t ∈ [0, Tk]) , for k = 1,… , d , be independent Lévy subordinators which are independent of W. The pointwise characteristic function of the subordinated GRF defined by L(x) ∶= W(l1(x1),… , ld(xd)) , for x = (x1,… , xd) ∈ [0, � ]d , admits the formula for � ∈ ℝ and any fixed point x = (x1,… , xd) ∈ [0, � ]d . Here, the variance function �2 W ∶ ℝ d + → ℝ+ is given by �2 W (x) ∶= qW (x, x) for x ∈ ℝ d + . In the one-dimensional case, the Lévy-Khinchin formula gives an explicit representation of the pointwise characteristic function of a Lévy process. This representation also applies to the subordinated Brownian motion, since it is itself a Lévy process (see Subsection 3.1). Note that in the construction of the subordinated Brownian motion one cannot replace the Brownian motion by a general one-parameter GRF on ℝ+ without losing the validity of the Lévy-Khinchin formula. Hence, in the case of a subordinated GRF on a higher-dimensional parameter space, it is natural that we have to restrict the class of admissible GRFs in order to obtain a Lévy-Khinchin-type formula which is the d-dimensional analogon of Theorem 2.1. We recap that the pointwise characteristic function of a standard Brownian motion B is given by mA(Z) → m(Z) ℙ − a.s. for A → +∞. �(g(W(Z))) = lim A→+∞ �(�A(g(W(Z)))) = lim A→+∞ �(mA(Z)) = �(m(Z)). �(g(W(Z))) = �(g+(W(Z))) − �(g−(W(Z))) = �(m̃+(Z)) − �(m̃−(Z)) = �(m(Z)), �(exp(i�L(x))) = �(exp(i�W(l1(x1),… , ld(xd)))) = � ( exp ( i�W (l1(x1),… , ld(xd)) − 1 2 �2�2 W (l1(x1),… , ld(xd)) )) , 2668 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 for t ≥ 0 . Note that the Brownian motion ist not characterized by this property, i.e. not every zero-mean GRF on ℝ+ with the above pointwise characteristic function is a Brownian motion, since this specific characteristic function can be attained by different covariance functions, whereas the covariance function of the Browian motion is given uniquely by qBM(s, t) = Cov(B(s),B(t)) = min{s, t} for s, t ≥ 0 (see for example (Schoutens 2003, Sec- tion 3.2.2)). Motivated by this, we impose the following assumptions on the GRF on ℝd + . Assumption 4.4 Let W = (W(x), x ∈ ℝ d + ) be a zero-mean continuous GRF. We assume that there exists a constant 𝜎 > 0 such that the pointwise characteristic function of W is given by for � ∈ ℝ and every x = (x1,… , xd) ∈ ℝ d + . Remark 4.5 Note that for a zero-mean, continuous and stationary GRF W̃ = (W̃(x), x ∈ ℝ d + ) , the GRF W defined by satisfies Assumption 4.4. We are now able to derive the Lévy-Khinchin-type formula for the subordinated GRF. Theorem  4.6  (Lévy‑Khinchin‑type formula) Let Assumption 4.4 hold. We assume inde- pendent Lévy subordinators lk = (lk(x), x ∈ [0, Tk]) , with Lévy triplets (�k, 0, �k) , for k = 1,… , d , are given corresponding to representation (1). Further, we assume that these processes are independent of the GRF W. We consider the subordinated GRF defined by L ∶ Ω × [0, 𝕋 ]d → ℝ with L(x) ∶= W(l1(x1),… , ld(xd)) for x = (x1,… , xd) ∈ [0, 𝕋 ]d . The pointwise characteristic function of the random field L is, for any x = (x1,… , xd) ∈ [0, � ]d , given by for � ∈ ℝ . Here, the jump measure �ext is defined through for a, b ∈ ℝ where the Lévy measure �# k , for k = 1,… , d and a, b ∈ ℝ , is given by �B(t)(�) = 𝔼(exp(i�B(t))) = exp ( − 1 2 t�2 ) , for � ∈ ℝ, �W(x)(�) = �(exp(i�W(x))) = exp ( − 1 2 �2�2(x1 +⋯ + xd) ) , W(x) ∶= √ x1 +⋯ + xdW̃(x), for x = (x1,… , xd) ∈ ℝ d + , �L(x)(�) = 𝔼(exp(i�W(l1(x1),… , ld(xd))) = exp ( − (x1,… , xd) ⋅ (�2�2 2 (�1,… , �d) t + � ℝ⧵{0} 1 − ei�z + i�z1{|z|≤1}(z)�ext(dz) )) , �ext([a, b]) ∶= ⎛⎜⎜⎝ �# 1 ([a, b]) ⋮ �# d ([a, b]) ⎞⎟⎟⎠ , �# k ([a, b]) ∶= ∫ ∞ 0 ∫ b a 1√ 2��2t exp � − x2 2�2t � dx �k(dt). 2669Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Proof It follows by (Sato 2013, Theorem 30.1 and Lemma 30.3) that the measures �# k are Lévy measures for k = 1, ..., d . For notational simplicity we prove the assertion for d = 2 . For general d ∈ ℕ the assertion follows by the same arguments. Claim 1: For a Lévy measure � on (ℝ+,B(ℝ+)) it holds for every � ∈ ℝ: where the measure 𝜈♯ is defined by 𝜈♯(I) = ∫ ∞ 0 ∫ b a 1√ 2𝜋t exp(− x2 2t )dx𝜈(dt) , for I = [a, b] with a, b ∈ ℝ . We use the notation fs(x) ∶= 1√ 2�s exp(− x2 2s ) for s > 0 and x ∈ ℝ and derive this equation by a direct calculation using the definition of the measure 𝜈♯: In the last step we used that the characteristic function of a N(0, s)-distributed random variable is given by �(�) = exp(− s�2 2 ) for � ∈ ℝ and s > 0 . Further, we used the fact that f � s (x) = −x∕sfs(x) to see that Claim 2: (See (Applebaum 2009), P. 53) For a Lévy subordinator l with triplet (� , 0, �) it holds for t ≥ 0 and 𝜉 > 0. With these two assertions at hand we can now prove the Lévy-Khinchin-type for- mula. The case � = 0 is trivial since both sides equal 1 in this case. Let (x, y) ∈ [0, � ]2 and 0 ≠ � ∈ ℝ be fixed. Using Lemma 4.1 and Remark 4.2 with g(⋅) = exp(i�⋅) and Z = (l1(x), l2(y)) we calculate where where we used Assumption 4.4. Therefore, using the independence of the processes l1 and l2 together with Claim 2 we obtain � ∞ 0 exp(− 𝜉2 2 y) − 1𝜈(dy) = � ℝ⧵{0} exp(i𝜉x) − 1 − i𝜉x1{|x|≤1}(x)𝜈♯(dx), � ℝ⧵{0} exp(i𝜉x) − 1 − i𝜉x1{|x|≤1}(x)𝜈♯(dx) = � ℝ⧵{0} (exp(i𝜉x) − 1 − i𝜉x1{|x|≤1}(x))� ∞ 0 fs(x)𝜈(ds)dx = � ∞ 0 � ℝ⧵{0} exp(i𝜉x)fs(x)dx − 1 − i𝜉 � 1 −1 xfs(x)dx𝜈(ds) = � ∞ 0 exp(− s𝜉2 2 ) − 1𝜈(ds). ∫ 1 −1 xfs(x) = −s(fs(1) − fs(−1)) = 0. �(exp(−�l(t))) = exp(−t(�� + ∫ ∞ 0 (1 − exp(−�y))�(dy))), �(exp(i�W(l1(x), l2(y)))) = �(m(l1(x), l2(y))) m(x�, y�) ∶= 𝔼(exp(i�W(x�, y�))) = exp(− 1 2 �2�2(x� + y�)) for (x�, y�) ∈ ℝ 2 + , 2670 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 where we define the (Lévy-)measures �̂�1 and �̂�2 by �̂�k([a, b]) = 𝜈k([a∕𝜎 2, b∕𝜎2]) for a, b ∈ ℝ+ and k = 1, 2 . Now, using Claim 1 we calculate where the measures �̂�♯ k for k = 1, 2 are given by: for a, b ∈ ℝ . This finishes the proof. ◻ Using Theorem  4.6 together with the convolution theorem (see for example (Klenke 2013, Lemma 15.11 (iv))) we immediately obtain the following corollary. Corollary 4.7 Let Assumption 4.4 hold. We assume d independent Lévy sub- ordinators lk = (lk(x), x ∈ [0, Tk]) are given for k = 1,… , d , which are inde- pendent of W and the corresponding Lévy triplets are given by (�k, 0, �k) for k = 1,… , d . We consider the subordinated GRF L ∶ Ω × [0, 𝕋 ]d → ℝ defined by L(x) ∶= W(l1(x1),… , lk(xd)), for x = (x1,… , xd) ∈ [0, � ]d . Further, we assume that inde- pendent Lévy processes l̃k on [0,Tk] are given with triplets (0, �2�k, � # k ) for k = 1,… , d in the sense of the one-dimensional Lévy-Khinchin formula, see Theorem 2.1. Here, the Lévy measure �# k is defined by for k = 1,… , d and a, b ∈ ℝ . The pointwise marginal distribution of the subordinated GRF satisfies for every x = (x1,… , xd) ∈ [0, � ]d. 𝜙L(x,y)(𝜉) = �(exp(− 1 2 𝜎2𝜉2l1(x)))�(exp(− 1 2 𝜎2𝜉2l2(x))) = exp(−x(𝛾1 𝜎2𝜉2 2 + ∫ ∞ 0 (1 − exp(− 𝜉2 2 y))�̂�1(dy))) ⋅ exp(−y(𝛾2 𝜎2𝜉2 2 + ∫ ∞ 0 (1 − exp(− 𝜉2 2 y))�̂�2(dy))), 𝜙L(x,y)(𝜉) = exp ( − x(𝛾1 𝜎2𝜉2 2 − � ℝ⧵{0} exp(i𝜉x) − 1 − i𝜉x1{|x|≤1}(x)�̂�♯1(dx))) − y(𝛾2 𝜎2𝜉2 2 − � ℝ⧵{0} exp(i𝜉x) − 1 − i𝜉x1{|x|≤1}(x)�̂�♯2(dx)) ) , �̂�♯ k ([a, b]) = ∫ ∞ 0 ∫ b a 1√ 2𝜋t exp(− x2 2t )dx�̂�k(dt) = ∫ ∞ 0 ∫ b a 1√ 2𝜋𝜎2t exp(− x2 2𝜎2t )dx𝜈k(dt), �# k ([a, b]) ∶= ∫ ∞ 0 ∫ b a 1√ 2��2t exp � − x2 2�2t � dx �k(dt), L(x) D =l̃1(x1) +⋯ + l̃d(xd), 2671Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 We point out that the case of stationary GRFs is excluded by Assumption 4.4. There- fore, we consider this situation in the following remark where we again assume d = 2 for notational simplicity. Remark 4.8 Let W be a stationary, centered GRF with covariance function qW ((x, y), (x �, y�)) = q̃W ((x − x�, y − y�)) , for (x, y), (x�, y�) ∈ ℝ 2 + , and pointwise variance 𝜎2 ∶= q̃W ((0, 0)) > 0 . Let l1 and l2 be independent Lévy subordinators, which are also inde- pendent of W. We obtain by Lemma 4.1 the following representation for the pointwise char- acteristic function of the subordinated random field defined by L(x, y) ∶= W(l1(x), l2(y)) , for (x, y) ∈ [0, � ]2: where which is a constant function in (x�, y�) . Therefore we obtain for (x, y) ∈ [0, � ]2 . Hence, in case of a stationary GRF, the subordinated GRF is pointwise normally distributed with variance �2. We conclude this subsection with a remark on the given Lévy-Khinchin formula. Remark 4.9 With the approach of subordinating GRFs on a higher-dimensional domain, we obtain a discontinuous Lévy-type random field and a Lévy-Khinchin formula which allows access to the pointwise distribution of the random field. Further we obtain a simi- lar parametrization of the class of subordinated random fields, as it is the case for Lévy processes on a one-dimensional parameter space: Under the assumptions of Theorem 4.6, every subordinated GRF can be characterized by the tuple (�2, �1,… , �d, �ext, qW ) , where qW ∶ ℝ d + ×ℝ d + → ℝ is the covariance function of the GRF. Further, the class of sub- ordinated GRFs is linear in the sense that for the sum of two independent subordinated GRFs one can construct a single subordinated GRF with the same pointwise characteristic function. 5 Covariance Function One advantage of the subordinated GRF is that the correlation between spatial points is accessible. The correlation structure is hereby determined by the covariance function of the underlying GRF and the specific choice of the subordinators. For statistical applica- tions it is often important to image or enforce a specific correlation structure in view of fitting random fields to physical phenomena. In this context the question arises whether one can find analytically explicit formulas for the covariance function of a subordinated Gaussian random field. This will be explored in the following section. �L(x,y)(�) = �(exp(i�W(l1(x), l2(y))) = �(m(l1(x), l2(y))), m(x�, y�) = �(exp(i�W(x�, y�))) = exp ( − 1 2 �2�2 ) , �L(x,y)(�) = exp ( − 1 2 �2�2 ) , 2672 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 For notational simplicity we restrict the dimension to be d = 2 in this section but we point out that analogous results apply for dimensions d ≥ 3 . A direct application of Lemma 4.1 yields the following corollary. Corollary 5.1 Let W be a continuous, zero-mean GRF on ℝ2 + . Further, let l1 and l2 be two independent Lévy subordinators which are independent of W. Then the subordinated GRF L defined by L(x, y) ∶= W(l1(x), l2(y)) , for (x, y) ∈ ℝ 2 + , is zero-mean with covariance function for (x, y), (x�, y�) ∈ ℝ 2 + , where qW ∶ ℝ 2 + ×ℝ 2 + → ℝ denotes the covariance function of the GRF W. Proof For (x, y) ∈ [0, � ]2 , we use Lemma 4.1 and the fact that the GRF W is centered to deduce �(L(x, y)) = �(W(l1(x), l2(y))) = 0 . Let (x, y), (x�, y�) ∈ [0, � ]2 be fixed. Another application of Lemma 4.1 with W̃(x1, y1, x2, y2) ∶= W(x1, y1) ⋅W(x2, y2) , g = id ℝ and Z ∶= (l1(x), l2(y), l1(x �), l2(y �)) yields the desired formula. ◻ 5.1 The Isotropic Case We use Corollary 5.1 to derive a semi-explicit formula for the covariance function of the sub- ordinated GRF, where the underlying GRF is isotropic. Lemma 5.2 Let W ∶ Ω ×ℝ 2 + → ℝ be a zero-mean, continuous and isotropic GRF with covariance function qW ((x, y), (x�, y�)) = q̃W (|x − x�|, |y − y�|) . Further, suppose that l1 and l2 are independent Lévy subordinators on [0,T1] (resp. [0,T2] ) with density functions f1 and f2 , i.e. f x 1 (⋅) (resp. ly 2 (⋅) ) is the density function of l1(x) (resp. l2(y) ) for (x, y) ∈ (0, � ]2 . The covariance function of the subordinated GRF L with L(x, y) ∶= W(l1(x), l2(y)) , for (x, y) ∈ [0, � ]2 , admits the representation for (x, y), (x�, y�) ∈ [0, � ]2 with x ≠ x′ and y ≠ y′. For x = x� and y ≠ y′ it holds for x ≠ x′ and y = y� one obtains and for (x, y) = (x�, y�) the pointwise variance is given by qL((x, y), (x �, y�)) ∶= �(L(x, y)L(x�, y�)) = � ( qW ((l1(x), l2(y)), (l1(x �), l2(y �))) ) , qL((x, y), (x �, y�)) = ∫ ℝ+ ∫ ℝ+ q̃W (s, t)f |x−x�| 1 (s)f |y−y�| 2 (t)dsdt, qL((x, y), (x, y �)) = ∫ ℝ+ q̃W (0, t)f |y−y�| 2 (t)dt, qL((x, y), (x �, y)) = ∫ ℝ+ q̃W (s, 0)f |x−x�| 1 (s)ds, Var(L(x, y)) = qL((x, y), (x, y)) = q̃W (0, 0). 2673Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Proof The assertion follows immediately by Corollary 5.1 together with the independence of the processes l1 and l2 and the fact that |lk(x) − lk(x �)|D=lk(|x − x�|) for x, x� ∈ [0, Tk] and k = 1, 2 by the definition of a Lévy process. ◻ 5.2 The Non‑isotropic Case In this subsection, we derive a formula for the covariance function of the subordinated GRF for the case that the underlying GRF is not isotropic. In the following, we use the notation x ∧ y ∶= min(x, y) and x ∨ y ∶= max(x, y) for real numbers x, y ∈ ℝ . The next lemma will be useful in the proof of the covariance representation. Lemma 5.3 Let l = (l(x), x ∈ [0, T]) be a general Lévy process with density func- tion f ∶ (0,T] ×ℝ → ℝ , i.e. the probability density function of the random varia- ble l(x) is given by f x(⋅) , for x ∈ (0,T] . In this case, the joint probability density func- tion of the random vector Z ∶= (l(x ∧ x�), l(x ∨ x�)) , with x ≠ x� ∈ (0,T] , is given by fZ(s, t) = fmin(x,x�)(s) ⋅ f |x�−x|(t − s) for t, s ∈ ℝ. Proof Let x, x� ∈ (0, T] with x < x′ and x1, x2 ∈ ℝ be fixed. The increment l(x�) − l(x) is stochastically independent of the random variable l(x), which yields For the case that x′ < x the same argument yields which finishes the proof. ◻ Remark 5.4 Note that Lemma 5.3 immediately implies that the joint density fZ(s, t) of the two-dimensional random vector Z = (l(x ∧ x�), l(x ∨ x�)) for a Lévy subordinator l with x ≠ x� ∈ (0,T] is given by With this lemma at hand we are able to derive a formula for the covariance func- tion of the subordinated (non-isotropic) GRF. Without loss of generality we consider points (x, y), (x�, y�) with x ≤ x′ and y ≤ y′ in the following Lemma. Formulas for the other cases follow by the same arguments with Lemma 5.3 and Remark 5.4. ℙ(l(x) ≤ x1 ∧ l(x�) ≤ x2) = 𝔼(1{l(x)≤x1}1{l(x�)≤x2}) = 𝔼(1{l(x)≤x1}1{l(x�)−l(x)≤x2−l(x)}) = � ℝ � ℝ 1{s≤x1}1{t≤x2−s}f x(s)f x �−x(t)dtds = � x1 −∞ � x2−s −∞ f x(s)f x �−x(t)dtds = � x1 −∞ � x2 −∞ f x(s)f x �−x(t − s)dtds. ℙ(l(x) ≤ x1 ∧ l(x�) ≤ x2) = � x1 −∞ � x2 −∞ f x � (s)f x−x � (t − s)dsdt, fZ(s, t) = fmin(x,x�)(s) ⋅ f |x�−x|(t − s), for s, t ∈ ℝ+. 2674 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Lemma 5.5 Let W ∶ Ω ×ℝ 2 + → ℝ be a zero-mean, continuous and non-isotropic GRF with covariance function qW . Further, suppose that l1 and l2 are independent Lévy subordinators on [0,T1] (resp. [0,T2] ) with density functions f1 and f2 , i.e. f x 1 (⋅) (resp. ly 2 (⋅) ) is the density function of l1(x) (resp. l2(y) ) for (x, y) ∈ (0, � ]2 . The covariance function of the subordin- ared GRF L with L(x, y) ∶= W(l1(x), l2(y)) , for (x, y) ∈ [0, � ]2 , admits the representation for (x, y), (x�, y�) ∈ (0, � ]2 with x < x′ and y < y′. For x = x� and y < y′ , it holds and for x < x′ and y = y� it holds For (x, y) = (x�, y�) one obtains for the pointwise variance of the field Proof Using Corollary 5.1, the independence of the processes l1 and l2 , Lemma 5.3 and Remark 5.4 we calculate for (x, y), (x�, y�) ∈ (0, � ]2 with x < x′ and y < y′: The remaining cases follow by the same argument. ◻ 5.3 Statistical Fitting of the Covariance Function The parametrization property of the subordinated GRF (see Remark 4.9) motivates a direct approach of covariance fitting: For a natural number N ∈ ℕ , we assume that discrete points {(xi, yi), i = 1,… ,N} are given with corresponding empirical covariance function data Cemp = {C emp i,j , i, j = 1,… ,N} , where Cemp i,j represents the empirical covariance of the field evaluated at the points (xi, yi) and (xj, yj) . We search for the solution to the problem qL((x, y), (x �, y�)) = ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ qW ((x1, x2), (x3, x4))f x 1 (x1)f y 2 (x2) × f x �−x 1 (x3 − x1)f y�−y 2 (x4 − x2)dx1 dx2 dx3 dx4, qL((x, y), (x, y �)) = ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ qW ((x1, x2), (x1, x4))f x 1 (x1)f y 2 (x2) × f y�−y 2 (x4 − x2)dx1 dx2 dx4, qL((x, y), (x �, y)) = ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ qW ((x1, x2), (x3, x2))f x 1 (x1)f y 2 (x2) × f x �−x 1 (x3 − x1)dx1 dx2 dx3. Var(L(x, y)) = qL((x, y), (x, y)) = ∫ ℝ+ ∫ ℝ+ qW (x1, x2, x1, x2)f x 1 (x1)f y 2 (x2)dx1dx2. qL((x, y), (x �, y�)) = ∫ ℝ 4 + qW ((x1, x2), (x3, x4))dℙ(l1(x),l2(y),l1(x �),l2(y �))(x1, x2, x3, x4) = ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ ∫ ℝ+ qW ((x1, x2), (x3, x4))f x 1 (x1)f y 2 (x2) × f x �−x 1 (x3 − x1)f y�−y 2 (x4 − x2)dx1 dx2 dx3 dx4. 2675Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 where we use the notation q̃L ∶= {qL(xi, yi), i, j = 1,… ,N} and ‖ ⋅ ‖∗ is an appropriate norm on ℝN , e.g. the euclidian norm. In order to solve this type of problem, the formulas for the covariance function given by Lemma 5.2 and Lemma 5.5 can be used, but still accessing the solution will be challenging due to the complexity of the set of admissible parameters. 6 Stochastic Regularity ‑ Pointwise Moments In this section we consider pointwise moments of a subordinated GRF L. In particular, we derive conditions which ensure the existence of pointwise p-th moments of the sub- ordinated GRF L defined by L(x) ∶= W(l1(x1),… , ld(xd)) , for x = (x1,… , xd) ∈ [0, � ]d. Obviously, in order to guarantee the existence of moments of the random variable L(x) , we have to impose conditions on the GRF W and the subordinators l1,… , ld . The following theorem gives a better insight into the interaction between the underlying GRF and the stochastic regularity of the subordinators and presents coupled regularity conditions on the tail behaviour of both components of the random field. Theorem  6.1 We assume that W is a centered and continuous GRF on ℝd + with covari- ance function qW ∶ ℝ d + ×ℝ d + → ℝ . Further, we assume that there exist a positive num- ber N ∈ ℕ , coefficients {cj, j = 1,… ,N} ⊂ [0,+∞) and d-dimensional exponents {𝛼(j), j = 1,… ,N} ⊂ ℝ d + such that the pointwise variance function �2 W of W satisfies Here, we use the notation z� = z �1 1 ⋅ ⋯ ⋅ z �d d for z = (z1,… , zd) ∈ ℝ d + and � = (� 1 ,… , � d ) ∈ ℝ d + . We consider a fixed point x ∈ [0, � ]d and assume that the densities f x1 1 ,… , f xd d of the evaluated processes l1(x1),… , ld(xd) fulfill with positive decay rates {�i, i = 1,… , d} . Here, the constants C and K are independent of z but may depend on the evaluation point x = (x1,… , xd) and �i may depend on xi , for i = 1,… , d . We define the number Then, the random variable L(x) admits a p-th moment for p ∈ [1, a) , i.e. L(x) ∈ L p(Ω;ℝ) for p ∈ [1, a). Proof Let Z ∼ N(0, �2) be a real-valued, centered, normally distributed random variable with variance 𝜎2 > 0 . It follows by Eq. (18) in Winkelbauer (2012) that the p-th absolute moment of Z admits the form �(|Z|p) = Cp� p, for all p > −1 , with a constant Cp depending on p. Let p ≥ 1 be a fixed number. We use Lemma 4.1 to calculate argmin � ‖q̃L − Cemp‖∗ ��� admissible tuples (𝜎2, 𝛾1, 𝛾2, 𝜈ext, qW ) � (2)�W (z) = qW (z, z) 1∕2 ≤ N∑ j=1 cjz �(j) , for z1,… , zd ≥ 0. (3)f xi i (z) ≤ C|z|−�i , for z ≥ K and i = 1,… , d, a ∶= min{(�i − 1)∕�(j) i || i = 1,… , d, j = 1,… ,N, �(j) i ≠ 0}. 2676 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 with for (x� 1 ,… , x� d ) ∈ ℝ d + . Hence, we obtain Next, we use the tail estmations (2) and (3), Hölder’s inequality and the independence of the subordinators to calculate It remains to show that all the integrals Ij i are finite. For i ∈ {1,… , d} and j ∈ {1,… ,N} with �(j) i = 0 we have Ij i = 1 . If �(j) i ≠ 0 it holds where the integral in the last step is finite since p𝛼(j) i − 𝜂i < −1 for all i ∈ {1,… , d} and j ∈ {1,… ,N} with �(j) i ≠ 0 . ◻ We close this section with three remarks on the assumptions and possible extensions of Theorem 6.1. Remark 6.2 The assumption given by Eq. (2) is, for example, fulfilled for the d-dimensional Brownian sheet with N = 1 , c1 = 1 and �(1) = (1∕2,… , 1∕2) ∈ ℝ d + . Condition (2) also accomodates the GRFs we considered in the Lévy-Khinchin formula (see Theorem 4.6 and Assumption 4.4) with N = d , c1 = ⋯ = cd = 1 and 𝛼(j) = 1∕2 ⋅ êj for j = 1,… , d , where êj is the j-th unit vector in ℝd . Further, this assumption is fulfilled for any stationary GRF W. Indeed, in case of a stationary GRF the assumption is satisfied for �(1) = (�, 0,… , 0) for any 𝜀 > 0 and, hence, Theorem 6.1 yields that every moment of the corresponding evalu- ated subordinated GRF exists, independently of the specific choice of the subordinators. This is consistent with Remark 4.8. The assumption on the Lévy subordinators in Eq. (3) is natural and can be verified easily in many cases, see also (Barth and Stein 2018b, Assump- tion 3.7 and Remark 3.8). For example, if for some non-negative integer n ∈ ℕ , the n-th derivative of the characteristic function �l(xi) (⋅) is integrable over ℝ , then Eq. (3) holds with �i = n , K = 0 and C = 1 2� ∫ +∞ −∞ | dn dtn �l(xi) (t)| dt (cf. (Hughett 1998, Lemma 12)). �(|L(x)|p) = �(|W(l1(x1),… , ld(xd))|p) = �(m(l1(x1),… , ld(xd))), m(x� 1 ,… , x� d ) ∶= �(|W(x� 1 ,… , x� d )|p) = Cp� p W (x� 1 ,… , x� d ), �(|L(x)|p) = Cp� ( �p W (l1(x1),… , ld(xd)) ) . 𝔼(|L(x)|p) = Cp𝔼 ( �p W (l1(x1),… , ld(xd) ) ≤ Cp � ℝ d + ( N∑ j=1 cjz �(j) )p f x1 1 (z1)… f xd d (zd)d(z1,… , zd) ≤ C(N, p) N∑ j=1 c p j d∏ i=1 � +∞ 0 z p�(j) i i f xi i (zi)dzi ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ =∶I j i . I j i = ( � K 0 +� +∞ K ) z p𝛼(j) i i f xi i (zi)dzi ≤ Kp𝛼(j) i + C � +∞ K z p𝛼(j) i −𝜂i i dzi < +∞, 2677Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Remark 6.3 We point out that the statement of Theorem 6.1 remains valid if we consider Lévy distributions with discrete probability distribution which satisfy a discrete version of (3): If the GRF W satisfies (2) and the evaluated discrete subordinators l1(x1),… , ld(xd) satisfy then we obtain that �(|L(x1,… , xd)|p) < ∞ for p ∈ [1, a) with the real number a defined in Theorem 6.1. Remark 6.4 For the pointwise existence of moments given by Theorem 6.1, it is not neces- sary to restrict the subordinating processes to the class of Lévy subordinators. More gener- ally, one could consider a GRF W satisfying (2) and general Lévy processes l1,… , ld satis- fying (3) for |z| ≥ K . In this case, Theorem 6.1 still holds for the random field L defined by L(x) ∶= W(|l1(x1)|,… , |ld(xd)|) , for x = (x1,… , xd) ∈ [0, � ]d. 7 Numerical Examples In the following, we present numerical experiments on the theoretical results given in this paper. The goal of this section is to use the knowledge on theoretical properties of the subordi- nated GRF to investigate existing numerical methods for the approximation of pointwise dis- tributions (Subsection 7.1) as well as methods to verify or disprove the existence of moments of a random variable (Subsection 7.2). The numerical methods may also be useful for a fitting of random fields to existing data in applications. All our numerical experiments have been performed with MATLAB. 7.1 Experiments on the Lévy‑Khinchin Formula The Lévy-Khinchin-type formula (Theorem 4.6) allows access to the pointwise distribution of a subordinated GRF which motivates the investigation of numerical methods to approximate the pointwise distribution. To be more precise, we use Corollary 4.7 to obtain a pointwise distributional representation of a subordinated GRF as the sum of one-dimensional Lévy pro- cesses with transformed Lévy triplets. We use this representation to investigate the perfor- mance of different methods to approximate the distribution of Lévy processes. Assume L = (W(l1(x), l2(y)), (x, y) ∈ [0, 1]2) is a subordinated GRF where the GRF W satisfies Assumption 4.4 and the two subordinators l1 and l2 are characterized by the Lévy triplets (�k, 0, �k) for k = 1, 2 . It follows by Corollary 4.7 that L admits the pointwise distribu- tional representation for (x, y) ∈ [0, 1]2 . Here, the processes l̃k on [0, 1] are independent Lévy processes with tri- plets (0, �2�k, � # k ) , for k = 1, 2 , in the sense of the one-dimensional Lévy-Khinchin formula (see Theorem 2.1) and the Lévy measure �# k is defined by (4)f xi i (k) = ℙ(li(xi) = k) ≤ C|k|−�i , for k ≥ K and i ∈ {1,… , d}, (5)L(x, y) D =l̃1(x) + l̃2(y), �# k ([a, b]) ∶= ∫ ∞ 0 ∫ b a 1√ 2��2t exp � − x2 2�2t � dx �k(dt), 2678 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 a, b ∈ ℝ and k = 1, 2 . We choose specific spatial points and use two different methods to approximate the distribution of the Lévy processes on the right hand side of (5): the com- pound Poisson approximation (CPA) (see (Schoutens 2003, Section 8.2.1)) and the Fourier inversion method for Lévy processes (see Gil-Pelaez (1951) and Barth and Stein (2018b)) which allows for a direct approximation of the density of the right hand side of (5). In order to investigate the performance of these two approaches, the corresponding results are then compared with samples of the subordinated GRF on the left hand side of Eq. (5). 7.1.1 Compound Poisson Approximation We recall that a Gamma(aG, bG) process lG has independent Gamma-distributed incre- ments and lG(t) follows a Gamma(aG ⋅ t, bG) distribution. In our first example we choose Gamma(4, 12) processes to subordinate the GRF W defined by W(x, y) = √ x + y W̃(x, y) , for (x, y) ∈ ℝ 2 + , where W̃ is a Matérn-1.5-GRF with pointwise standard deviation � = 2 (see Remark 4.5). We fix the evaluation point (x, y) = (1, 1) and use the CPA method to obtain samples of the Lévy process on the right hand side of (5) which can then be com- pared with samples of the subordinated GRF. Figure 2 (left and middle) shows the corre- sponding histograms for 10.000 samples of each distribution. We observe an accurate fit of the samples generated by the different approaches: the first histogram, corresponding to the exact sampling of the subordinated GRF, displays the same characteristics as the histogram generated by CPA, which shows that the CPA method is appropriate to simulate the distribution of the right hand side of Eq. (5). 7.1.2 Fourier Inversion Method The second approach is to approximate the density function of the right hand side of (5) by the Fourier inversion (FI) method (see Gil-Pelaez (1951) and Barth and Stein (2018b)) and compare it with samples of the subordinated GRF. Figure 3 illustrates the results for this approach where we used the evaluation point (x, y) = (1, 1) , the same GRF as in Subsec- tion 7.1.1, Gamma(4, 12) subordinators and 100.000 samples of the subordinated GRF. As one can see in Fig. 3, the approximated density of the right hand side of (5) per- fectly matches the pointwise distribution of the subordinated GRF. We want to confirm this observation by a Kolmogorov-Smirnov-Test (see for example (Pestman 2009, Section VII.4)). Figure 4 illustrates how the empirical CDF, obtained by sampling the subordinated GRF, converges to the target CDF which is approximated by the Fourier inversion method using Eq. (5). A Kolmogorov-Smirnov-test with 10.000 samples and a level of significance of 5% is passed. Fig. 2 Samples of the subordinated GRF W(l 1 (1), l 2 (1)) (left), the sum of the corresponding transformed Lévy processes l̃ 1 (1) + l̃ 2 (1) generated by the CPA method (middle) and both histograms in one plot (right) 2679Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 In the next experiment we use a modified subordinator, which results in a less smooth pointwise density of the subordinated GRF. We repeat the experiment with Gamma(0.5, 10) subordinators where the GRF, the evaluation point and the sample size remain unchanged. Figure 5 shows 100.000 samples of the subordinated GRF and the density of the process given by the right hand side of Eq. (5) approximated via the Fourier Inversion method. As in the first experiment, the results given by Fig.  5 indicate that the approximated density of the right hand side of (5) matches the pointwise distribution of the subordinated GRF. Figure  6 illustrates how the empirical CDF, obtained by sampling of the subordi- nated GRF, converges to the approximated target CDF of the right hand side of Eq. (5), which is computed by the Fourier inversion method. A Kolmogorov-Smirnov-test with a level of significance of 5% is passed. 7.2 Pointwise Moments Theorem 6.1 guarantees the existence of pointwise moments of the subordinated GRF if the GRF and the corresponding subordinators satisfy certain conditions. In the following numerical experiments, we investigate the results of different statistical methods to inves- tigate numerically the existence of pointwise moments of a certain order in the specific situation of the subordinated Gaussian random field. We set d = 2 and assume W to be a Fig. 3 Samples of Gamma(4, 12)-subordinated GRF and approximated density (FI) Fig. 4 Approximated target CDF (FI) vs. empirical CDF using 100 (left), 1.000 (middle) and 10.000 (right) samples of the subordinated GRF with Gamma(4, 12) subordinators 2680 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Brownian sheet on ℝ2 + . Further, we use Lévy processes with different stochastic regularity - in terms of the existence of moments - to subordinate the GRF W. 7.2.1 Statistical Methods to Test the Existence of Moments of a Random Variable The existence of moments of a specific distribution is one of the most frequently formu- lated assumptions in statistical applications. For example, already the strong law of large numbers assumes finiteness of the first moment of the corresponding random variable. Nevertheless, in the literature only few statistical methods exist to verify or disprove the existence of moments, given a specific sample of random variables (see e.g. Mandelbrot (2012); Hill (1975); Ng and Yau (2018); Fedotenkov (2013a, 2014, 2013b)). One of the earlier methods to verify the existence of moments of a distribution was proposed in 1963 by Mandelbrot (see Mandelbrot (2012) and Cont (2001)). It is based on the simple observa- tion that the estimated (sample-)moments will converge to a certain value for an increasing sample size if the theoretical moment exists. On the other side, if the theoretical moment does not exist, the estimated moment will diverge or behave unstable when the sample size increases. However, this quite intuitive method is rather heuristic and depends highly on the experience of the researcher (see also Fedotenkov (2013b)). Another popular direct way to investigate the existence of moments of a certain distribution is the sample-based Fig. 5 Samples of Gamma(0.5, 10)-subordinated GRF and approximated density (FI) Fig. 6 Approximated target CDF (FI) vs. empirical CDF using 100 (left), 1.000 (middle) and 10.000 (right) samples of the subordinated GRF with Gamma(0.5, 10) subordinators 2681Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 estimation of a decay rate � for the corresponding density function proposed by Hill in Hill (1975). However, the Hill-estimator requires a parameter k > 0 which specifies the sample values which are considered as the tail of the distribution and it turned out that the Hill- estimator is very sensitive to the choice of this parameter k. Further, the method makes the quite restrictive assumption that the underlying distribution is of Pareto-type (see Ng and Yau (2018); Fedotenkov (2013b, 2014, 2013a)). In 2013, Fedotenkov proposed a bootstrap test for the existence of moments of a given distribution (see Fedotenkov (2013b)). The test performs well for specific distributions, however, its accuracy deteriorates fast when moments of higher order are considered (see also Fedotenkov (2014)). Recently, Ng and Yau proposed another sample-based bootstrap test for the existence of moments which outperforms the previously mentioned methods for many distributions (see Ng and Yau (2018)). The test is based on a result from bootstrap asymptotic theory which states that the m out of n bootstrap sample mean (see Bickel et al. (1997)) converges weakly to a normal distribution. For a detailed description of the test statistic and further theoretical investiga- tions we refer the interested reader to Ng and Yau (2018). Based on these observations, we investigate the results of direct moment estimation via Monte Carlo (MC) and the bootstrap test proposed by Ng and Yau to analyze the existence of (pointwise) moments of the subordinated GRF. For our numerical examples we choose three different Lévy distributions to subordi- nate the Brownian sheet W: a Poisson distribution, a Gamma distribution and a Student-t distribution. Therefore, we use a discrete and a continuous distribution where all moments are finite and a continuous distribution, which only admits a limited number of moments. Hence, we consider three fundamentally different situations. In all three experiments, we consider the evaluation point x = (x1, x2) = (1, 1) ∈ ℝ 2 + for the subordinated GRF L. Note that the two-dimensional Brownian sheet satisfies Eq. (2) in Theorem  6.1 with N = 1 , c1 = 1 and �(1) = (1∕2, 1∕2). 7.2.2 Poisson‑subordinated Brownian sheet In this example, we use Poisson(3) processes to subordinate the two-dimensional Brownian sheet. It is easy to verify that condition (4) is satisfied for any 𝜂i > 0 , i = 1, 2 , since point evaluations of a Poisson process are Poisson distributed. Theorem 6.1 implies the existence of the p-th moment of the evaluated field L(1, 1) for any p < ∞ (see Remark 6.3). We esti- mate the p-th moment for p ∈ {4, 6, 8} by a MC-estimation using M samples of the evalu- ated GRF L(1, 1) for different values of M ∈ ℕ , i.e. where (L(i)(1, 1), i ≥ 1) are i.i.d. samples of the evaluated field L(1,  1). As explained in Subsection 7.2.1, the MC-estimator EM(|L(1, 1)|p) is expected to converge for M → ∞ if the p-th moment exists and one expects unstable behaviour if this is not the case. Figure 7 shows the development of the MC-estimator EM(|L(1, 1)|p) for the p-th moment as a func- tion of the number of samples M. For every moment, we take 5 independent MC-runs to validate that they converge to the same value. �(|L(1, 1)|p) ≈ EM(|L(1, 1)|p) = 1 M M∑ i=1 |L(i)(1, 1)|p, 2682 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 As expected, Fig. 7 shows a stable convergence of the MC-estimator for a growing num- ber of samples for every considered moment. Further, the different independent MC-runs converge to the same value - the theoretical p-th moment for p ∈ {4, 6, 8}. In the next step, we perform the bootstrap test (see Subsection  7.2.1 and Ng and Yau (2018)). We test the existence of the p-th moment for p ∈ {1, 2, 3, 4, 5, 6, 7, 8} using M = 107 samples of the subordinated evaluated GRF L(1, 1). Hence, the null and alter- native hypothesis are given by for the different values of p. We choose the significance level �s = 1% and perform 100 independent test runs. Figure 8 shows the proportion of acceptance of the null hypothesis in the 100 test runs as a function of the considered moment p ∈ {1, 2, 3, 4, 5, 6, 7, 8}. As we see in Fig. 8, the bootstrap test accepts the null hypothesis H0 in almost every test run for every considered moment p ∈ {1, 2, 3, 4, 5, 6, 7, 8} which is in line with our expectatations since all these moments exist. We conclude that both approaches, the MC moment estimation and the bootstrap test, perform as expected in this experiment. H0 ∶ �(|L(1, 1)|p) < +∞ vs. H1 ∶ �(|L(1, 1)|p) = +∞, Fig. 7 Five independent realizations of the MC-estimator EM(|L(1, 1)|p) ≈ �(|L(1, 1)|p) as a function of the sample numbers M with a Poisson(3)-subordinated Brownian sheet; p = 4 (left), p = 6 (middle), p = 8 (right) Fig. 8 Results for 100 independ- ent runs of the bootstrap test for the existence of the p-th moment using Poisson(3) processes to subordinate the Brownian sheet 2683Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 7.2.3 Gamma‑subordinated Brownian sheet In our second numerical example we consider Gamma processes to subordinate the Brown- ian sheet. We recall that, for aG, bG > 0 , a Gamma(aG, bG)-distributed random variable admits the density function where Γ(⋅) denotes the Gamma function. A Gamma process (l(t))t≥0 has independent Gamma distributed increments and l(t) follows a Gamma(aG ⋅ t, bG)-distribution for t > 0 . Therefore, condition (3) holds for any 𝜂i > 0 , for i = 1, 2 and, hence, Theorem 6.1 again implies the existence of every moment, i.e. �(|L(1, 1)|p) < ∞ for any p ≥ 1 . We choose aG = 4 , bG = 10 and estimate the p-th moment of L(1, 1) with p ∈ {4, 6, 8} by a MC-esti- mation using a growing number of samples M ∈ ℕ . Figure 9 shows the development of the MC-estimator EM(|L(1, 1)|p) for the p-th moment as a function of the number of samples M. As in the first experiment, we take 5 independent MC-runs to validate the convergence to a unique value. In line with our expectations, the results show a stable convergence of the MC-estimations for the different moments of this subordinated GRF. In this experiment we again perform the bootstrap test for the existence of the p-th moment for p ∈ {1, 2, 3, 4, 5, 6, 7, 8} using M = 107 samples of the subordinated evaluated GRF L(1, 1). Hence, the null and alternative hypothesis are given by for the different values of p. We choose the significance level �s = 1% and perform 100 independent test runs. Figure 10 shows the proportion of acceptance of the null hypothesis in the 100 test runs as a function of the considered moment p ∈ {1, 2, 3, 4, 5, 6, 7, 8}. As in the first experiment, the test results meet our expectations, since almost every test run accepts the null hypothesis for any moment p ∈ {1, 2, 3, 4, 5, 6, 7, 8}. 7.2.4 Student t‑subordinated Browinan Sheet In our last experiment we consider a Lévy process where the pointwise distribution only admits a finite number of moments. The Student’s t-distribution with three degress of free- dom admits the density function x ↦ b aG G Γ(aG) xaG−1 exp(−xbG), for x > 0, H0 ∶ �(|L(1, 1)|p) < +∞ vs. H1 ∶ �(|L(1, 1)|p) = +∞, Fig. 9 Five independent realizations of the MC-estimator EM(|L(1, 1)|p) ≈ �(|L(1, 1)|p) as a function of the sample numbers M with a Gamma(4, 10)-subordinated Brownian sheet; p = 4 (left), p = 6 (middle), p = 8 (right) 2684 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 It follows by (Kelker (1971), Theorem 3) that a Student-t distributed random variable with three degrees of freedom is infinitely divisible. Hence, we can define Lévy processes lj , for j = 1, 2 , such that lj(1) follows a Student-t distribution with three degrees of freedom (see Sato (2013, Theorem 7.10)). Using these processes and the Brownian sheet W, we consider the subordinated GRF L(x1, x2) ∶= W(|l1(x1)|, |l2(x2)|) for (x1, x2) ∈ [0,T1] × [0, T2] (see Remark 6.4). For our numerical experiment we again evaluate the field at (x1, x2) = (1, 1) . Using (6) we obtain Therefore, condition (3) is satisfied for �i = 4 , for i = 1, 2 , and it is violated for any 𝜂i > 4 (see also Remark 6.4). Since the Brownian sheet satisfies condition (2) with N = 1 , c1 = 1 and �(1) = (1∕2, 1∕2) , Theorem 6.1 yields that �(|L(1, 1)|p) < ∞ for p < 6 and we expect that this boundary is sharp, i.e. we expect that �(|L(1, 1)|p) = ∞ for p ≥ 6. We estimate the p-th moment for p ∈ {5, 6, 8} with the MC-estimator EM(|L(1, 1)|p) with growing sample number M ∈ ℕ . In Fig. 11 we see the development of the MC-estimator (6)ft(x) = Γ(2)√ 3�Γ(3∕2) � 1 + x2 3 �−2 , for x ∈ ℝ. ft(x) ≤ C|x|−4, for x ∈ ℝ. Fig. 10 Results for 100 inde- pendent runs of the bootstrap test for the existence of the p-th moment using Gamma(4,10) processes to subordinate the Brownian sheet Fig. 11 Five independent realizations of the MC-estimator EM(|L(1, 1)|p) ≈ �(|L(1, 1)|p) as a function of the sample numbers M with a Student-t-subordinated Brownian sheet; p = 5 (left), p = 6 (middle), p = 8 (right) 2685Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 for the p-th moment as a function of the number of samples. For every moment, we take 5 independent MC-runs. The results indicate a convergence of the MC-estimations of the p-th moment for p = 5 : in this case the estimation stabilizes with growning sample size and all 5 independent MC-estimations seem to converge to a unique value. However, for the higher moments p = 6 and p = 8 , we see upward breakouts and instable behaviour of the corresponding MC-estimator for increasing sample sizes. Further, the 5 independent MC-runs do not indicate a convergence to a unique value. For all the considered moments p ∈ {5, 6, 8} , these results are in line with our expectations, since the p-th moment of the evaluated subordinated GRF L(1, 1) admits a p-th moment for p < 6 and this boundary is sharp (see Theorem 6.1). We perform the bootstrap test for the existence of the p-th moment for p ∈ {1, 2, 3, 4, 4.5, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.5, 7, 8} using M = 107 samples of the subordi- nated GRF L(1, 1). Hence, the null and alternative hypothesis are again given by for the different values of p. We choose the significance level �s = 1% and perform 100 independent test runs. Figure 12 shows the proportion of acceptance of the null hypothesis in the 100 test runs as a function of the considered moment p and the test statistic values for the 100 test runs. In all of the 100 test runs the null hypothesis is accepted for p ∈ {1, 2, 3, 4, 4.5, 5} . Further, in almost all of the 100 test runs H0 is rejected for the cases p ∈ {6, 6.5, 7, 8} , which is absolutely in line with the theoretical results for this specific choice of the subordinated GRF. Only for p ∈ (5, 6) , the test rejects the null hypothesis in some of the test runs although the theoretical moment exists. Overall, the test results for the exist- ence of moments of the Student-t-subordinated GRF match our expectations based on Theorem 6.1. Acknowledgements The authors would like to thank two anonymous referees for their remarks, which led to a significant improvement of the manuscript. Funding Open Access funding enabled and organized by Projekt DEAL. Funded by Deutsche Forschun- gsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2075 - 390740016.  H0 ∶ �(|L(1, 1)|p) < +∞ vs. H1 ∶ �(|L(1, 1)|p) = +∞, Fig. 12 Results for 100 independent runs of the bootstrap test for the existence of the p-th moment using Student-t-distributed random variables as subordinators: share acceptance of H 0 (left), test statistic values for the test runs (right) 2686 Methodology and Computing in Applied Probability (2022) 24:2661–2688 1 3 Data Availability All data used in our experiments have been produced with MATLAB random number generators and no external datasets have been used. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request. Declarations Conflicts of Interest The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com- mons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. References Adler RJ, Taylor JE (2007) Random fields and geometry. Springer Monographs in Mathematics. Springer, New York ISBN 978-0-387-48112-8 Aliprantis CD, Border KC (2006) Infinite dimensional analysis, 3rd edn. Springer, Berlin. ISBN 978-3-540- 32696-0; 3-540-32696-0. A hitchhiker’s guide Applebaum D (2009) Lévy processes and stochastic calculus, Cambridge studies in advanced mathematics, vol 116, 2nd edn. Cambridge University Press, Cambridge. ISBN 978-0-521-73865-1. https:// doi. org/ 10. 1017/ CBO97 80511 809781 Barndorff-Nielsen OE, Pedersen J, Sato K (2001) Multivariate subordination, self-decomposability and sta- bility. Adv in Appl Probab 33(1), 160–187 Barth A, Merkle R (2020) Subordinated Gaussian random fields in elliptic partial differential equations. ArXiv e-prints, arXiv: 20110 9311 [mathPR] To appear in stochastics and partial differential equations: analysis and computations Barth A, Stein A (2018a) A study of elliptic partial differential equations with jump diffusion coefficients. SIAM/ASA J Uncertain Quantif 6(4):1707–1743. https:// doi. org/ 10. 1137/ 17M11 48888 Barth A, Stein A (2018b) Approximation and simulation of infinite-dimensional Lévy processes. Stoch Partial Differ Equ Anal Comput 6(2):286–334. ISSN 2194-0401. https:// doi. org/ 10. 1007/ s40072- 017- 0109-2 Bastian P (2014) A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure. Comput Geosci 18(5):779–796. ISSN 1420-0597. https:// doi. org/ 10. 1007/ s10596- 014- 9426-y Bickel PJ, Götze F, van Zwet WR (1997) Resampling fewer than n observations: gains, losses, and remedies for losses. vol 7. Empirical Bayes, sequential analysis and related topics in statistics and probability (New Brunswick, NJ, 1995), pp 1–31 Buchmann B, Kaehler B, Maller R, Szimayer A (2017) Multivariate subordination using generalised gamma convolutions with applications to variance gamma processes and option pricing. Stochastic Process Appl 127(7):2208–2242. ISSN 0304-4149. https:// doi. org/ 10. 1016/j. spa. 2016. 10. 008 Buchmann B, Lu KW, Madan DB (2019) Weak subordination of multivariate Lévy processes and variance generalised gamma convolutions. Bernoulli 25(1):742–770. ISSN 1350-7265. https:// doi. org/ 10. 3150/ 17- bej10 04 Buchmann B, Lu KW, Madan DB (2020) Self-decomposability of weak variance generalised gamma convo- lutions. Stochastic Process Appl 130(2):630–655. ISSN0304-4149. https:// doi. org/ 10. 1016/j. spa. 2019. 02. 012 Cont R (2001) Empirical properties of asset returns: stylized facts and statistical issues. Quantitative Finance 1:223–236 Da Prato J G Zabczyk (2014) Stochastic equations in infinite dimensions, encyclopedia of mathematics and its applications, vol 152, 2nd edn. Cambridge University Press, Cambridge. ISBN 978-1-107-05584-1. https:// doi. org/ 10. 1017/ CBO97 81107 295513 2687Methodology and Computing in Applied Probability (2022) 24:2661–2688 http://creativecommons.org/licenses/by/4.0/ https://doi.org/10.1017/CBO9780511809781 https://doi.org/10.1017/CBO9780511809781 http://arxiv.org/abs/201109311 https://doi.org/10.1137/17M1148888 https://doi.org/10.1007/s40072-017-0109-2 https://doi.org/10.1007/s10596-014-9426-y https://doi.org/10.1007/s10596-014-9426-y https://doi.org/10.1016/j.spa.2016.10.008 https://doi.org/10.3150/17-bej1004 https://doi.org/10.3150/17-bej1004 https://doi.org/10.1016/j.spa.2019.02.012 https://doi.org/10.1016/j.spa.2019.02.012 https://doi.org/10.1017/CBO9781107295513 1 3 Dobrushin RL (1979) Gaussian and their subordinated self-similar random generalized fields. Ann Probab 7(1):1–28.  ISSN 0091-1798. http:// links. jstor. org/ sici? sici= 0091- 1798(197902) 7: 1<1: GATSS R>2. 0. CO;2- Y& origin= MSN Fedotenkov I (2013a) A bootstrap method to test for the existence of finite moments. J Nonparametr Stat 25(2):315–322. ISSN 1048-5252. https:// doi. org/ 10. 1080/ 10485 252. 2012. 752487 Fedotenkov I (2013b) A simple nonparametric test for the existence of finite moments. Econometrics: Econometric and Statistical Methods - General eJournal Fedotenkov I (2014) A note on the bootstrap method for testing the existence of finite moments. Statistica 74(4):447–453. https:// ideas. repec. org/a/ bot/ rivsta/ v74y2 014i4 p447- 453. html Gil-Pelaez J (1951) Note on the inversion theorem. Biometrika 38:481–482. ISSN 0006-3444. https:// doi. org/ 10. 1093/ biomet/ 38.3- 4. 481 Graham IG, Kuo FY, Nichols JA, Scheichl R, Schwab C, Sloan IH (2015) Quasi-Monte Carlo finite element methods for elliptic PDEs with lognormal random coefficients. Numer Math 131(2):329–368.  ISSN 0029-599X. https:// doi. org/ 10. 1007/ s00211- 014- 0689-y Hill BM (1975) A simple general approach to inference about the tail of a distribution. Ann Statist 3(5):1163– 1174. ISSN 0090-5364. http:// links. jstor. org/ sici? sici= 0090- 5364(197509) 3: 5< 1163: ASGAT I>2. 0. CO;2- Y& origin= MSN Hughett P (1998) Error bounds for numerical inversion of a probability characteristic function. SIAM J Numer Anal 35(4):1368–1392. ISSN 0036-1429. https:// doi. org/ 10. 1137/ S0036 14299 63108 5X Kelker D (1971) Infinite divisibility and variance mixtures of the normal distribution. Ann Math Statist 42:802–808. ISSN 0003-4851. https:// doi. org/ 10. 1214/ aoms/ 11776 93436 Klenke A (2013) Wahrscheinlichkeitstheorie. Springer Berlin Heidelberg. https:// doi. org/ 10. 1007/ 978-3- 642- 36018-3 Leonenko NN, Ruiz-Medina MD, Taqqu MS (2017) Rosenblatt distribution subordinated to Gaussian ran- dom fields with long-range dependence. Stoch Anal Appl 35(1):144–177. ISSN 0736-2994. https:// doi. org/ 10. 1080/ 07362 994. 2016. 12307 23 Makogin V, Spodarev E (2021) Limit theorems for excursion sets of subordinated Gaussian random fields with long-range dependence. Stochastics 0(0):1–32. https:// doi. org/ 10. 1080/ 17442 508. 2021. 19146 20 Mandelbrot B (2012) The variation of certain speculative prices [reprint of J. Bus. 36 (1963), no. 4, 394– 419]. In: Financial risk measurement and management, Internat. Lib. Crit. Writ. Econ., vol 267. Edward Elgar, Cheltenham, pp 230–255 Ng WL, Yau CY (2018) Test for the existence of finite moments via bootstrap. J Nonparametr Stat 30(1):28– 48. ISSN 1048-5252. https:// doi. org/ 10. 1080/ 10485 252. 2017. 14028 96 Pestman WR (2009) Mathematical statistics, 2nd edn. De Gruyter Textbook, Walter de Gruyter & Co., Ber- lin. ISBN 978-3-11-020852-8. https:// doi. org/ 10. 1515/ 97831 10208 535 Protter PE (2004) Stochastic Integration and differential equations, applications of mathematics (New York), vol 21, 2nd edn. Springer-Verlag, Berlin, stochastic Modelling and Applied Probability. ISBN 3-540-00313-4 Sato K (2013) Lévy processes and infinitely divisible distributions, Cambridge studies in advanced math- ematics, vol 68. Cambridge University Press, Cambridge, translated from the 1990 Japanese original, Revised edition of the 1999 English translation. ISBN 978-1-107-65649-9 Schoutens W (2003) Lévy processes in finance: pricing financial derivatives. Wiley Series in Probability and Statistics, Wiley. ISBN 9780470851562. https:// books. google. de/ books? id= Be4Nq Iz7h- kC Werner D (2011) Funktionalanalysis. Springer-Lehrbuch, Springer Berlin Heidelberg.  ISBN 9783642210174. https:// books. google. de/ books? id= 580Nc HHrAX sC Winkelbauer A (2012) Moments and absolute moments of the normal distribution. ArXiv e-prints, arXiv: 12094 340v2 [mathST] Zhang D, Kang Q (2004) Pore scale simulation of solute transport in fractured porous media. Geophysical Research Letters 31. https:// doi. org/ 10. 1029/ 2004G L0198 86 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 2688 Methodology and Computing in Applied Probability (2022) 24:2661–2688 https://www.jstor.org/stable/2242835 https://www.jstor.org/stable/2242835 https://doi.org/10.1080/10485252.2012.752487 https://ideas.repec.org/a/bot/rivsta/v74y2014i4p447-453.html https://doi.org/10.1093/biomet/38.3-4.481 https://doi.org/10.1093/biomet/38.3-4.481 https://doi.org/10.1007/s00211-014-0689-y https://www.jstor.org/stable/2958370 https://www.jstor.org/stable/2958370 https://doi.org/10.1137/S003614299631085X https://doi.org/10.1214/aoms/1177693436 https://doi.org/10.1007/978-3-642-36018-3 https://doi.org/10.1007/978-3-642-36018-3 https://doi.org/10.1080/07362994.2016.1230723 https://doi.org/10.1080/07362994.2016.1230723 https://doi.org/10.1080/17442508.2021.1914620 https://doi.org/10.1080/10485252.2017.1402896 https://doi.org/10.1515/9783110208535 https://books.google.de/books?id=Be4NqIz7h-kC https://books.google.de/books?id=580NcHHrAXsC http://arxiv.org/abs/12094340v2 http://arxiv.org/abs/12094340v2 https://doi.org/10.1029/2004GL019886 On Some Distributional Properties of Subordinated Gaussian Random Fields Abstract 1 Introduction 2 Preliminaries 2.1 Lévy Processes 2.2 Gaussian Random Fields 3 The Subordinated Gaussian Random Field 3.1 Motivation: The Subordinated Brownian Motion 3.2 The Definition of the Subordinated Gaussian Random Field 3.3 Measurability 4 The Pointwise Distribution of the Subordinated GRF and the Lévy-Khinchin-formula 5 Covariance Function 5.1 The Isotropic Case 5.2 The Non-isotropic Case 5.3 Statistical Fitting of the Covariance Function 6 Stochastic Regularity - Pointwise Moments 7 Numerical Examples 7.1 Experiments on the Lévy-Khinchin Formula 7.1.1 Compound Poisson Approximation 7.1.2 Fourier Inversion Method 7.2 Pointwise Moments 7.2.1 Statistical Methods to Test the Existence of Moments of a Random Variable 7.2.2 Poisson-subordinated Brownian sheet 7.2.3 Gamma-subordinated Brownian sheet 7.2.4 Student t-subordinated Browinan Sheet Acknowledgements References