{'fraction_non_alphanumeric': 0.05631174533479693, 'fraction_numerical': 0.025649469447493595, 'mean_word_length': 3.929743867243867, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Nonparametric Analysis of Multivariate Time Series: A Matrix Gamma Process Approach

While there is an increasing amount of literature about Bayesian time series analysis, only a few Bayesian nonparametric approaches to multivariate time series exist. Most methods rely on Whittle's Likelihood, involving the second order structure of a stationary time series by means of its spectral density matrix. This is often modeled in terms of the Cholesky decomposition to ensure positive definiteness. However, asymptotic properties such as posterior consistency or posterior contraction rates are not known. A different idea is to model the spectral density matrix by means of random measures. This is in line with existing approaches for the univariate case, where the normalized spectral density is modeled similar to a probability density, e.g. with a Dirichlet process mixture of Beta densities. In this work, we present a related approach for multivariate time series, with matrix-valued mixture weights induced by a Hermitian positive definite Gamma process. The proposed procedure is shown to perform well for both simulated and real data. Posterior consistency and contraction rates are also established.


Introduction
With the surge of research in Bayesian nonparametrics over the last two decades, several nonparametric Bayesian approaches to analysing univariate stationary time series have been developed, such as Carter and Kohn (1997), Gangopadhyay et al. (1999), Liseo et al. (2001), Choudhuri et al. (2004a), Hermansen (2008), Chopin et al. (2013), Cadonna et al. (2017) and Edwards et al. (2018). These are generally based on Whittle's likelihood (Whittle, 1957), essentially the product of the asymptotically independent discrete Fourier transformed observations with variances equal to the spectral density at the corresponding Fourier frequencies, providing an approximation to the true likelihood. The Whittle likelihood avoids matrix inversion (for Gaussian time series) and significantly speeds up computation. Even for non-Gaussian time series, it provides a reasonable approximation and asymptotically correct inference in many situations (Shao and Wu, 2007). The nonparametric prior on the spectral density used in Choudhuri et al. (2004a) is based on the Bernstein polynomial approximation first employed in Petrone (1999) for nonparametric density estimation. Whereas no asymptotic results were provided for other nonparametric priors such as those based on state-space methodology by Carter and Kohn (1997) or fitting local polynomials to the log spectral density in Gangopadhyay et al. (1999), Choudhuri et al. (2004a) proved consistency of the posterior distribution based on the Bernstein polynomial mixtures and the Whittle likelihood for Gaussian time series.
The Whittle likelihood for multivariate time series is a function of the spectral density matrix, a Hermitian positive definite (Hpd) matrix-valued function. Several frequentist and Bayesian extensions of nonparametric spectral density estimation to multivariate time series have been suggested in the recent literature (Dai and Guo (2004), Rosen and Stoffer (2007) and Li and Krafty (2018)). To ensure positive definiteness, these approaches are based on smoothing the components of the Cholesky decomposition of the periodogram matrix. Rosen and Stoffer (2007) proposed a Bayesian nonparametric approach based on fitting smoothing splines to each real and imaginary component of the Cholesky decomposition of the periodogram matrix. Similarly, Zhang (2016) proposed an extension to nonstationary multivariate time series by an adaptive estimation of time-varying spectra and cross-spectra based on the Cholesky decomposition of the inverse of the spectral density matrix. Zhang (2018) described an extension to quantile-based spectra based on the Cholesky decomposition of the copula spectral density kernels. However, using the Cholesky decomposition of the spectral density matrix will make it difficult to elicit priors. While simulations studies and applications to real data generally demonstrate good performance of these extensions, no asymptotic results for any Bayesian nonparametric approaches to multivariate time series are known. Therefore, with the aim of proving posterior consistency of a nonparametric Bayesian approach to multivariate time series, we extend the Bernstein-Dirichlet process prior of Choudhuri et al. (2004a) from the space of positive real-valued functions to the space of Hpd matrix-valued functions. We employ the equivalence of the Dirichlet process and the normalized Gamma process used to define the mixture weights of the Bernstein polynomial mixture and generalize the completely random measure (CRM) on X = [0, π] with independent Gamma increments to a completely random Hpd matrix-valued measure with independent matrix-Gamma increments. To this end, we first define an infinitely divisible Hpd Gamma distribution in terms of the Lévy-Khinchine representation of its Laplace transform, following Pérez-Abreu and Stelzer (2014). Then, in analogy to the Gamma process, an Hpd Gamma process is defined using the Kingman construction (Kingman, 1992) of a CRM with Hpd Gamma increments. A software implementation of the methodology is included in the R package beyondWhittle, which is available on the Comprehensive R Archive Network (CRAN), see Meier et al. (2018).
The paper is organized as follows: Section 2 first introduces an infinitely divisible Hpd Gamma distribution defined via the Lévy-Khinchine representation of its Laplace transform. It then details the construction of the Hpd Gamma process using the Kingman construction based on a Poisson process with a suitably defined intensity measure corresponding to the Lévy measure of a Hpd Gamma distribution. It is proven that this defines a completely random Hpd measure with independent Hpd Gamma increments. Furthermore, we derive an almost surely convergent series representation of the Hpd Gamma process. Section 3 provides a hierarchical prior for the spectral density matrix of a multivariate stationary time series using the mixture of Bernstein polynomial prior with matrix-valued weights induced by a Hpd Gamma process. The posterior distribution is obtained by combining the prior with the Whittle likelihood which is sampled by an efficient Inverse Lévy measure algorithm as in Wolpert and Ickstadt (1998). In Section 4, we essentially follow the steps in Choudhuri et al. (2004a) and first show mutual contiguity of Whittle likelihood and full Gaussian likelihood under certain conditions on the Gaussian stationary time series. We then prove posterior consistency in the L 1 norm and strengthen the result to uniform consistency under the additional assumption of uniform boundedness of the Lipschitz constants of the spectral density. Section 4 also establishes posterior contraction rates in the Hellinger distance. The performance of the proposed nonparametric Bayesian approach is illustrated in Section 5 in a simulation study where it is compared to a parametric vector-autoregressive model. It is then applied to the analysis of the bivariate monthly time series of the Southern Oscillation Index and fish recruitment from 1950-1987, previously analyzed by Rosen and Stoffer (2007) and Shumway and Stoffer (2010). Section 6 concludes with an outlook to the analysis of non-Gaussian time series, frequentist coverage probabilities of Bayesian credible sets and alternative applications of the Hpd Gamma process.

Hpd Gamma Process
In this section, we construct an Hermitian positive definite (Hpd) Gamma process that will be used in the Bayesian procedure described in Section 3. We start our considerations by revisiting the Hpd Gamma distribution from Pérez-Abreu and Stelzer (2014). Let us first introduce some notation. For a complex d×d matrix A, let |A| := | det A|. Denote the trace of A by tr A, let etr A := exp(tr A) and denote the Hermitian conjugate by A * . The real-and imaginary part of A will be denoted by A and A.

Infinitely divisible Hpd Gamma distribution
Let α be a finite measure onS + d and let β :S + d → (0, ∞) be a measurable function, such that the integral S + d log(1 + 1/β(U ))α(dU ) is finite. The Hpd Gamma distribution Ga d×d (α, β) with parameters α, β is defined in terms of the Lévy-Khinchine representation of the Laplace transform as follows: holds for all Θ ∈S + d . Under the above assumptions on α and β, the Ga d×d (α, β) distribution is well-defined and it holds P (X ∈S + d ) = 1. Observe that for d = 1, the Hpd Gamma distribution coincides with the Gamma distribution. The parameters α, β are generalizations of the scale and rate parameter. Since the Ga d×d (α, β) distribution is defined in terms of the Lévy-Khinchine representation of its Laplace transform, it is necessarily infinitely divisible (see Theorem 1 and Remark 2 in Pérez-Abreu and Rosiński (2007)). This constitutes the key property for the upcoming Hpd Gamma process construction. We emphasize that e.g. the complex Wishart distribution (being another well-known generalization of the Gamma distribution to the Hpd matrix case) lacks this property (Lévy, 1948). Observe that the Lévy measure of the Ga d×d (α, β) dis- (2.1) Remark 2.1. A special case of the Ga d×d (α, β) distribution is the so-called AΓ distribution, that has been considered in Pérez-Abreu and Stelzer (2014) and generalized to the Hpd setting in Meier (2018), Section 2.4. To elaborate, the AΓ(η, ω, Σ) distribution is defined with the parameters η > d−1, ω > 0 and Σ ∈ S + d as the Ga d×d (α η,Σ , β Σ ) distribution, with where Γ denotes the Gamma function andΓ d the complex multivariate Gamma function (see Mathai and Provost (2005)), and β Σ (U ) = tr(Σ −1 U ). It has the advantage that for X ∼ AΓ(η, ω, Σ), the formulas for mean and covariance structure are explicitly known: where H = d i,j=1 H i,j ⊗ H j,i and H i,j being the matrix having a one at (i, j) and zeros elsewhere, see Meier (2018), Lemma 2.8. Thus the AΓ-distribution is particularly well suited for Bayesian prior modeling if the prior knowledge is given in terms of mean and covariance structure.

Process construction
We can now utilize the infinitely divisible Hpd Gamma distribution to define an Hpd Gamma process, i.e. a stochastic process with independent Ga d×d (α, β) distributed increments. The process construction is based on Poisson processes and generalizes the famous Kingman construction of the Gamma process (see Section 8.2 in Kingman (1992)) to the Hpd matrix case. Let us briefly recall the notion of Poisson processes. For a Borel space Y, a Poisson process Π on Y is a countable subset such that for all m > 0 and all disjoint subsets A 1 , . . . , A m , it holds that #{Π ∩ A 1 }, . . . , #{Π ∩ A m } are independent random variables, distributed as Poi(ν(A j )) for j = 1, . . . , m. Here, ν is a measure on Y, which is called the mean measure of Π and we write Π ∼ PP(ν). For a rigorous treatment of Poisson processes, the reader is referred to Kingman (1992).
Let X be a Polish space, equipped with a locally compact, σ-finite and nontrivial Borel measure. To define a Ga d×d (α, β) process on X , we allow the distributional parameters α, β to vary among X . To elaborate, denote by B For x ∈ X , the measure ν(dU , dr|x) := 1 r exp(−β(x, U )r)drα(x, dU ) onS + d corresponds to the Lévy measure of the Ga d×d (α(x, ·), β(x, ·)) distribution, see (2.1).
In what follows, we will make the following assumption on ν from (2.2): This property ensures that ν is a feasible Poisson process mean measure (see Section 2.5 in Kingman (1992)). Let Π ∼ PP(ν) and define the process Φ : The following result shows that Φ is well-defined and an independent increment process, with Ga d×d distributed increments.
Theorem 2.2. Let Φ be defined as in (2.4), with Π ∼ PP(ν) and ν from (2.2) fulfilling assumption (2.3). Then it holds: In the situation of Theorem 2.2, we write Φ ∼ GP d×d (α, β) and call Φ an Hpd Gamma process. If the process parameter β does not vary among X , i.e. β(x, ·) = β(·) for a functionβ :S + d → (0, ∞), then Φ is called homogeneous. In this case, part (a) of Theorem 2.2 reveals that Φ . Part (b) shows that Φ generalizes the notion of a completely random measure (i.e. a random measure Φ such that for any finite disjoint collection A 1 , . . . , A m ⊂ X , the random variables Φ(A 1 ), . . . , Φ(A m ) are independent, see Section 8 in Kingman (1992)). This is why we will call Φ a completely random Hpd measure.
The following result shows that Φ obeys a convenient almost surely convergent series representation involving iid components. This is related to the famous stick-breaking representation of the Dirichlet and the Gamma process (Sethuraman, 1994;Roychowdhury and Kulis, 2015) and will be of great usefulness for later practical applications, in particular for the implementation of MCMC algorithms.
Theorem 2.3. Let the assumptions of Theorem 2.2 be fulfilled and assume additionally that C α := X α(x,S + d )dx is finite. With α * := α/C α , assume that β(x, U ) ≥ β 0 holds for α * -almost all (x, U ) and some constant β 0 > 0. Then The proof is analogous to Section 3(B) in Rosiński (2001) and based on the Interval Theorem and the Marking Theorem for Poisson processes (see Kingman (1992)). It may be noted that the assumptions of Theorem 2.3 can be generalized slightly to hold outside a nullset in X . A detailed version of the proof can be found in Meier (2018), proof of Lemma 3.13.
For later proofs in the time series applications, further distributional properties of the Hpd Gamma process are needed, such as full support and lower probability bounds. These results may be of independent interest for different applications as well and can be found in Lemma A.1 and Lemma A.2 in the Appendix.

Spectral Density Inference
We will now illustrate how the Hpd Gamma process can be incorporated in a nonparametric prior model for the spectral density matrix f of a stationary multivariate time series. To elaborate, let I j,k := ((j − 1)π/k, jπ/k] for k > 0 and j = 1, . . . , k. Denote by b(x|j, k − j + 1) = Γ(k + 1)Γ(j) −1 Γ(k − j + 1) −1 x j−1 (1−x) k−j for 0 ≤ x ≤ 1 the density of the Beta(j, k−j+1) distribution. We will consider the following Bernstein-Hpd-Gamma prior for f : where α is a (nonrandom) measure onS + d and β :S + d → (0, ∞), such that (2.3) is fulfilled. Under this prior, the spectral density is modeled as a Bernstein polynomial mixture of (random) degree k, where the polynomial mixture weights are Hermitian positive semidefinite matrices induced by the Hpd Gamma process Φ. The latter is defined on X = [0, π] and the prior probability mass function p(k) of the polynomial degree is fully supported on N. A common choice (which is motivated from asymptotic considerations, as outlined in the upcoming Section 4), is p(k) ∝ exp(−ck log k) for k ∈ N, where c is a positive constant. This Bernstein polynomial mixture approach is inspired by existing methods for univariate spectral density inference (Choudhuri et al., 2004a) and density estimation (Petrone, 1999) and relies on the approximation properties of Bernstein polynomials (see e.g. Section 1.6 in Lorentz (2012) or Appendix E in Ghosal and van der Vaart (2017)).
The model specification is completed by employing Whittle's likelihood P n W for the spectral density. It is defined in terms of the Fourier coefficients of the data Z 1 , . . . , Z n . Asymptotically, the Fourier coefficients are independent and normally distributed, a result that is well known and holds true in quite great generality, beyond Gaussian time series, see Hannan (1970), p.224. Motivated from this asymptotic result, Whittle's likelihood is a pseudo likelihood that mirrors the asymptotic distribution of the Fourier coefficients.
With N := n/2 − 1, it is defined in terms of the Lebesgue density The coefficients corresponding to the Fourier frequencies ω = 0 and ω = π (the latter occurring for n even) are omitted in the likelihood, since they represent the mean and alternating mean resp., which are both typically treated separately from the autocovariance structure. By Bayes' Theorem, this leads to the posterior distribution P n W (f |Z 1 , . . . , Z n ) ∝ P n W (Z 1 , . . . , Z n |f )P (f ).
Since the posterior is not tractable analytically, we employ a Markov Chain Monte Carlo (MCMC) algorithm to draw random samples from it. To this end, we approximate the infinite series representation of Φ from Theorem 2.3 by a finite sum as Φ ≈ L j=1 δ xj r j U j for some large integer L, serving as a truncation parameter. In practice, the choice of L should depend on the sample size n. In chapter 5 below we choose L = max{20, n 1/3 }, which is the same value that has also been used by Choudhuri et al. (2004a) and Meier (2018), Section 3.4.3 and Section 5.2. With this approximative representation of Φ, the spectral density f is parametrized by the finite-dimensional vector Θ f = (k, x 1 , . . . , x L , U 1 , . . . , U L , r 1 , . . . , r L ). (3.4) Posterior samples can be drawn with a Metropolis-within-Gibbs sampler, which is discussed in more detail in Appendix B.

Asymptotic Properties
We will establish L 1 -consistency (Theorem 4.5 in Section 4.2) as well as uniform consistency (Theorem 4.6 in Section 4.3) and Hellinger contraction rates (Theorem 4.8 in Section 4.4) for the posterior of f . As a first important observation, we will derive that Whittle's likelihood and the full Gaussian likelihood are mutually contiguous (Theorem 4.2 in Section 4.1). Our considerations are based on the following assumption on the true spectral density f 0 .
Assumption 4.1. Let {Z t : t ∈ Z} be a Gaussian stationary time series in R d with spectral density f 0 fulfilling: (a) There exist constants b 0 , b 1 > 0 such that where λ 1 , λ d denote the smallest and largest eigenvalue respectively.
Part (a) of Assumption 4.1 has a statistical interpretation in terms of decay of linear dependence coefficients, at least for the univariate case (see Theorem 1.6 in Bradley (2002)). Part (b) can be seen as a regularity condition, since it implies f to be continuously differentiable with derivative being Hölder of order a − 1 > 0.

Contiguity
We will consider the following versionP n W =P n W (·|f ) of Whittle's likelihood, where the coefficients corresponding to the Fourier frequencies ω = 0 and ω = π are taken into account: with p n W from (3.3) and p 0,n (·|f ) being the density of the N d (0, 2πf (0)) distribution and p n/2,n (·|f ) the density of the N d (0, 2πf (π)) distribution, the latter being included if and only if n is even. The joint distribution ofZ 0 , . . . ,Z n/2 will be denoted byP n =P n (·|f ). The following result generalizes the findings from Choudhuri et al. (2004b) to the multivariate case. Recall that two sequences (P n ) and (Q n ) of measures on measurable spaces X n are called mutually contiguous, if for every sequence (A n ) of measurable sets it holds: Theorem 4.2. Let Assumption 4.1 be fulfilled. ThenP n W andP n are mutually contiguous.
The result carries over to the version P n W of Whittle's likelihood from (3.3), as formulated in the following Corollary.
Corollay 4.3. Let Assumption 4.1 be fulfilled and denote by P n = P n (·|f ) the joint distribution ofZ 1 , . . . ,Z N . Then P n W and P n are mutually contiguous. The proof of contiguity relies on the Gaussianity assumption and it can be shown that the result does not extend beyond Gaussianity in general (see the discussion in Section 4.3 in Meier (2018). It is conjectured that posterior consistency and contraction rates may be valid for non-Gaussian time series as well, however, since the proof relies on the contiguity result, it is unclear how this can be shown.

L 1 -consistency
Consistency for the spectral density of a univariate time series under a Bernstein-Dirichlet prior in the L 1 -topology has been derived in Choudhuri et al. (2004a). The authors relied on the assumption that the normalizing constant of the spectral density is known a priori. Our method of proof is slightly different and does not rely on such an assumption, but requires (the eigenvalues of) the spectral density to be uniformly bounded a priori. To elaborate, denote by C the set of continuousS + d -valued functions on [0, π]. We endow C with the σalgebra induced by the maximum Frobenius norm f F,∞ := max 0≤ω≤π f (ω) , where we will also consider the Frobenius L p norms f F,p := ( f (ω) p dω) 1/p for p > 0. Let 0 ≤ τ 0 < τ 1 ≤ ∞ and consider the set (4.1) Denoting by P the Bernstein-Hpd-Gamma prior (3.1) on f ∈ C, let P τ0,τ1 be the restriction of P to C τ0,τ1 . The pseudo-posterior when updating P τ0,τ1 with Whittle's likelihood P n W from (3.3) will be denoted by P n W ;τ0,τ1 (f |Z 1 , . . . , Z n ). We will make the following assumptions on Φ ∼ GP d×d (α, β) and on the prior probability mass function p(k) of the Bernstein polynomial degree k.
holds for all k ∈ N.
Part (a) and (b) of Assumption 4.4 are needed to ensure that prior probability mass is allocated in arbitrarily small neighborhoods of the true spectral density f 0 . In fact, it ensures that the Hpd Gamma process Φ assigns positive mass to neighborhoods of the true spectral measure F 0 := f 0 (ω)dω. Part (c) of Assumption 4.4 states full prior support for k and poses a condition on the decay of the tail of the prior. This is similar to frequentist tuning parameters that control the smoothness of the estimator, and is the same condition as in the univariate case, see Choudhuri et al. (2004a). Now we can formulate the first main result of this section, stating L 1 consistency of the posterior for Gaussian time series.
Theorem 4.5. Let {Z t } be a stationary Gaussian time series in R d with true spectral density f 0 fulfilling Assumption 4.1. Let τ ∈ (b 1 , ∞). Consider the prior P 0,τ on f , with Bernstein-Hpd-Gamma prior fulfilling Assumption 4.4. Then for all ε > 0 it holds with P n (·, f 0 ) denoting the joint distribution ofZ 1 , . . . ,Z N from (3.2).

Uniform consistency
The posterior consistency result can be strengthened from the L 1 topology to the uniform topology, under strengthened assumptions of the prior. The idea is to restrict the prior not only in the · F,∞ topology, but also introduce a uniform bound on the Lipschitz constants. To elaborate, consider the following Lipschitz norm on the space C 1 of continuously differentiable Hpd matrix valued functions on [0, 2π]: For τ > 0, letC and denote byP τ the restriction of the Bernstein-Hpd-Gamma prior (3.1) toC τ . Denote byP n W,τ the pseudo-posterior distribution for f when employingP τ in conjunction with Whittle's likelihood P n W from (3.3) Theorem 4.6. Let {Z t } be a stationary Gaussian time series in R d with true spectral density f 0 fulfilling Assumption 4.1 such that . Consider the priorP τ on f , with Bernstein-Hpd-Gamma prior fulfilling Assumption 4.4. Then for all ε > 0 it holds

Posterior contraction rates
Consider the following average squared Hellinger distance d 2 n,H (f 0 , f ) between two spectral density matrices f 0 , f : dx denotes the squared Hellinger distance between two probability densities p, q and denotes the density of the complex multivariate normal CN d (0, 2πf (ω j )) distribution. Observe that p j,N corresponds to the distribution of the Fourier co-efficientZ j under Whittle's likelihood (3.3). It can readily be seen that d n,H is a semimetric (i.e. a symmetric nonnegative function satisfying the triangle inequality, but possibly lacking the identity of indiscernibles) on C. We will need the following assumptions on the prior.
Part (a) and (b) of Assumption 4.7 are stronger than part (a) and (b) of Assumption 4.4. They ensure not only positive prior mass around f 0 , but also enable the derivation of lower bounds for the prior probability mass of small neighborhoods. As an example, the prior choice in Section 5 fulfills these conditions. Further examples are discussed in Meier (2018), Lemma 3.9. Part (c) is a condition on the prior tail of k. In contrast to part (c) of Assumption 4.7, the decay is bounded not only from below but also from above. The following theorem establishes contraction rates for the posterior of f . The rates coincide with the rates that are known for the univariate case, see Ghosal and Van Der Vaart (2007) and Example 9.19 in Ghosal and van der Vaart (2017).
Theorem 4.8. Let {Z t } be a stationary Gaussian time series in R d with true spectral density f 0 fulfilling Assumption 4.1 with 1 < a ≤ 2. Let 0 < τ 0 < b 0 < b 1 < τ 1 < ∞ and let the prior on f be given by P τ0,τ1 , with Bernstein-Hpd-Gamma prior fulfilling Assumption 4.7.
Then with ε n = n −a/(2+2a) (log n) (1+2a)/(2+2a) it holds in P n (·, f 0 probability, for every positive sequence M n with M n → ∞. Remark 4.9. In Theorem 4.8, the root average squared Hellinger distance d n,H is used for technical reasons, as the metric in consideration needs to yield exponentially powerful and uniformly exponentially consistent tests for this proof to be valid. It is known that such tests always exist for the Hellinger topology, while deriving such tests for different topologies may constitute a difficult task. Under additional prior assumptions, the Hellinger contraction rates from Theorem 4.8 can be used to establish rates in different topologies. As an example, if the prior is restricted to a Lipschitz class with uniformly bounded Lipschitz constant, then Lemma A.19 reveals that ε n is also a contraction rate with respect to the · F,1 norm and that ε 2/3 n is a contraction rate in the · F,∞ norm. It is, however, not clear whether these rates could be improved (or the prior assumptions relaxed) by a different proof technique.

Illustration
We will compare the performance of our proposed method with a parametric VAR model, where we will consider both simulated (Section 5.1) and a real data example (Section 5.2). An implementation of all procedures presented below is included in the R package beyondWhittle, which is available on CRAN, see Meier et al. (2018). As Bayes estimates, we will consider the pointwise posterior median functionf 0 , withf 0 (ω) consisting of the posterior median of real-and imaginary parts of the components of f (ω) for 0 ≤ ω ≤ π. We will compute pointwise 90% posterior credibility regions as the area between the pointwise 0.05 and 0.95 quantiles of the real and imaginary parts of the components of f . We will also consider uniform credibility regions. To elaborate, assume that we have a posterior sample f (1) , . . . , f (M ) at hand, as e.g. generated with the algorithm described in Appendix B. Denote by H : S + d → R d 2 the transformation that maps each Hermitian positive definite matrix A = (a rs ) d r,s=1 to a vector HA consisting of the logarithmized diagonal elements log a 11 , . . . , log a dd and the (non-logarithmized) real and imaginary parts of the entries {a rs : r < s} above the diagonal. For the transformed versions h Hf (j) for j = 1, . . . , M , denote the pointwise sample median function byĥ := (ĥ 1 , . . . ,ĥ d 2 ). Letσ := (σ 1 , . . . ,σ d 2 ) withσ r (ω) being the median absolute deviation of {h is called a uniform 90% credibility region. By construction, it holds that f ∈ C uni (·|0.9) with (empirical) posterior probability of at least 90%.

Prior choice
The parameters of the Hpd Gamma process in the Bernstein-Hpd-Gamma prior (3.1) are chosen such that the mean is large and the covariance is large, and proportional to the identity matrix. This is done in order to achieve a prior for f that is vague, homogeneous and isotropic (i.e. not preferring any particular directions in the function space). To elaborate, we choose the process parameters α(x, dU ) = α 0 (dU ) and β(x, dU ) ≡ β 0 of GP d×d (α, β) as α 0 (dU ) = 2dU and β 0 = 10 4 . As mentioned in Remark 2.1, we have that for Z ∼ Ga d×d (α 0 , β 0 ) it holds EZ = 10 4 I d with CovZ := E[Z ⊗Z]−(EZ ⊗EZ) being component-wise proportional to 10 8 I d .
To achieve a more stable mixture behavior at the left and right boundary of [0, π], the beta densities b(·|j, k −j +1) in (3.1) are replaced by their truncated and dilated counterparts. These are defined, for 0 < ξ l < ξ r < 1, as and have the advantage that the polynomial mixture at x = 0 (and x = 1 resp.) is not only determined by b(·|1, k) (and b(·|k, 1) resp.), but other basis functions are also incorporated. See Section 5.1 in Meier (2018) for a more detailed discussion of this matter. We will employ the truncated Bernstein polynomial basis from (5.1) with ξ l = 0.1 and ξ r = 0.9. The prior of the polynomial degree k is p(k) ∝ exp(−0.01k log k). The values of k are thresholded at k max = 500, which was found to be large enough in preliminary pilot runs, for the sake of computational speed-up (see Appendix B). Posterior samples are obtained with the MCMC algorithm from Appendix B. Each Markov chain is run for a total of 80,000 iterations, where the first 30,000 iterations are discarded as burn-in period and the remaining 50,000 iterations are thinned by a factor of 5, yielding a posterior sample size of 10,000. This procedure will be referred to as the NP procedure in the following, where NP stands for nonparametric method.

Simulated data
We consider simulated data drawn from the following bivariate VAR(2) model which has been considered in Rosen and Stoffer (2007): where WN d (µ, Σ) denotes d-dimensional White Noise with mean µ ∈ R d and covariance matrix Σ ∈ S + d , which is generated from a standard White Noise {ẽ t } withẽ t = (ẽ t,1 , . . . ,ẽ t,d ) andẽ t,1 , . . . ,ẽ t,d iid ∼ WN(0, 1) as e t := Σ 1/2ẽ t + µ. Furthermore, we draw from the following bivariate VMA(1) model: We will consider normally distributed innovations (i.e. e t,j iid ∼ N (0, 1)). Furthermore, we will also consider Student-t distributed innovations with ν = 4 degrees of freedom and centered exponential innovations (i.e. e t,j +1 iid ∼ Exp (1)). For N = n/2 −1, we compare the L 1 -error f 0 −f 0 1 := N −1 N j=1 f 0 (ω j )−f 0 (ω j ) of the pointwise posterior median functionf 0 . This is close to f 0 (ω)−f 0 (ω) dω for N large. The L 2 -error is defined analogously. As a parametric comparison model, we employ a Gaussian Vector Autoregession (VAR) with Normal-Inverse-Wishart prior Σ ∼ Wish −1 d×d (10 −4 , 10 −4 I d ) (see Section 3.4 in Gupta and Nagar (1999)) and vec(B 1 , . . . , B p ) ∼ N pd 2 (0, 10 4 I pd 2 ), where vec denotes the vectorization operator that stacks all columns of a matrix below each other. The order p is determined in a preliminary model selection step based on Akaike's Information Criterion (Akaike, 1974). To draw posterior samples, the Gibbs sampling algorithm from Section 2.2.3 in Koop and Korobilis (2010) is employed. The Markov Chain lengths and posterior sample size are chosen as for the NP procedure. We will refer to this method as the VAR procedure in the following. We consider M = 500 independent realizations of models (5.2) as well as (5.3) for each length n = 256, 512, 1024.
The results are shown in Table 1. As for the L 1 -and L 2 -error, it can be seen that the VAR procedure outperforms the NP procedure for VAR(2) data, which can be expected since in this case the parametric model is well-specified. For   VMA(1) data however, the NP procedure yields better results illustrating the benefit of employing a nonparametric procedure in comparison to a parametric method by being much less susceptible to misspecification. The empirical coverage of uniform 90% credibility regions is smaller for the NP procedure than for the VAR procedure in all examples (e.g. 0.37 vs. 0.90 for VAR(2) data or 0.44 vs. 0.98 for VMA(1) data with n = 512). With view on the discussion in Szabó et al. (2015), we conjecture that this property is due to the usage of Bernstein polynomials, which are known to have suboptimal approximation rates and tend to produce over-smoothed results (see Appendix E in Ghosal and van der Vaart (2017) and the simulation study in Edwards et al. (2018)).
Possible alternative approaches include the usage of a different polynomial basis (e.g. B-splines as in Edwards et al. (2018)) or to employ a parametric working model that increases the prior flexibility, as considered in Kirch et al. (2018) for the univariate case.

Southern Oscillation Index
We analyze the Southern Oscillation Index (SOI) and Recruitment series that have been analyzed in Rosen and Stoffer (2007) and Shumway and Stoffer (2010). Both time series are available as datasets soi and rec in the R package astsa (Stoffer, 2017). They consist of monthly data for 452 months. The SOI is defined as the normalized difference in air pressure between Tahiti (French Polynesia) and Darwin (Northern Territory, Australia). It constitutes a key indicator for warming or cooling effects of the central and eastern Pacific ocean known as El Niño and La Niña, see Bureau of Meteorology of the Australian Government (2018a,b). The Recruitment Series consists of the number of new spawned fish in a population in the Pacific Ocean. Since the fish are known to spawn better in colder waters (Rosen and Stoffer, 2007), it can be expected that there exists a cross-correlation between SOI and Recruitment. The rescaled and centered version of the data that we analyze is shown in Figure 1 (a).
The spectral inference results of the NP procedure are visualized in Figure 2 (a). It can be seen that a spectral peak at the frequency ω yearly := 2/12 ≈ 0.52 corresponding to a temporal distance of 12 months is estimated in the individual spectrum (f 11 ) of the SOI time series. This is in line with the previous findings from Rosen and Stoffer (2007). We would like to emphasize that the pointwise and uniform posterior credibility regions have to be interpreted with care, since we saw in the case of simulated data that they are in general far from being honest frequentist confidence sets.
To get a deeper insight into the dependency structure beyond the annual peak, we also investigate the data differenced at lag 12, i.e. the time series Y t := Z t − Z t−12 . The transformed data is visualized in Figure 1 (b) and the corresponding results of the NP procedure in Figure 2 (b). It can be seen that the annual peak at ω yearly is not present anymore. However, since the lag 12 differenced model does not constitute a perfect fit, the harmonics of the difference operator are visible in terms of six wavy bumps, particularly prominent in the individual spectra. On the other hand, the peak in the co-spectrum in the region of ω ≈ ω yearly /4 which is already visible in the original data, is even more prominent in the differenced version. This effect has also been observed by the authors in Rosen and Stoffer (2007), who explained it by a possible El Niño effect.

Conclusion
In this work, we presented a new nonparametric prior for the spectral density matrix of a stationary multivariate time series. Based on a mutual contiguity result, we established L 1 -consistency and Hellinger contraction rates for Gaussian time series under the Hpd Gamma process based prior and Whittle's likelihood. Regarding future work, it will be of interest to gain a better understanding of frequentist coverage properties of Bayesian credible sets, in particular if honest uncertainty quantification can be achieved, i.e. credible regions that asymptotically match confidence regions. Furthermore, it will be interesting to see if the asymptotic properties (at least consistency) carry over to different time series models beyond Gaussianity. The main difficulty in this setting will be the lack of a mutual contiguity result. Even for Gaussian data, it may be investigated whether the prior truncation of the eigenvalues can be dropped or at least relaxed towards bounds that are asymptotically growing to infinity (and shrinking to 0 respectively). In this work, we focused our attention on spectral density inference. However, the Hpd Gamma process based approach can be applied to other models as well.
As an example, one may consider semiparametric models, where in a parametric model of interest (e.g. linear regression, change point model), the time series constitutes the noise parameter, which can be modeled nonparametrically. From a theoretical perspective, it will be of interest whether a Bernstein-von-Mises result can be established for the parameter of interest. First considerations in this direction have been outlined in Section 9 in Meier (2018).

Acknowledgements
The foundation of this work was laid when Alexander Meier was visiting the University of Auckland. The visit has been supported by a travel scholarship (Kurzstipendien für Doktoranden) of the German Academic Exchange Service (DAAD). This work was supported by DFG grant AZ KI 1443/3-1. Renate Meyer was also supported by the James Cook Research Fellowships from Government funding, administered by the Royal Society Te Apārangi.

Appendix A: Proofs
We will make extensive use of matrix algebra. For a comprehensive overview on the most important calculation rules and results, the reader may confer Appendix B.1 in Meier (2018).

A.1. Proofs of Section 2
Detailed versions of the proofs in this section can be found in Meier (2018), Sections 3.1 and 3.2.
Let Π 1 , . . . , Π m be independent with Π j ∼ PP(ν j ). By the Superposition Theorem (see Section 2.2 in Kingman (1992)), it follows that Π is equal in distribution to ∪ k j=1 Π j . Since the result follows.
The following result concerns full support of the Hpd Gamma process. It will be used later in Lemma A.6 to show full support of the Bernstein-Hpd-Gamma prior (as needed for posterior consistency). For a probability measure µ, we denote by supp(µ) the support of µ and similarly, for a random variable X, let supp(X) denote the support of the distribution of X.
The results from Lemma A.1 can be refined, quantifying the probability mass that increments of the Hpd Gamma process puts in small neighborhoods under slightly stronger assumptions. This is done in the following Lemma A.2. It will be used later in Lemma A.14 and Lemma A.15 to derive lower probability bounds for the Bernstein-Hpd-Gamma prior (as needed for posterior contraction rates).

A.2. Proofs of Section 4.1
For a detailed proof of all the results in this section, the reader is referred to Section 4.2 and Section 4.3 in Meier (2018). To establish Theorem 4.2 with methods from probability theory, we translate the complex-valued version of Whittle's likelihood into a real-valued version, which is obtained by representing the real and imaginary parts separately. Indeed, it can be shown thatP n W can equivalently be written as being the block diagonal with the blocks f (0), Bf (ω 1 ), . . . , Bf (ω N ) (and f (π) for n even) and the algebra isomorphism see p.224 in Hannan (1970). Similarly, the full Gaussian likelihoodP n can be written in the frequency domain in terms ofZ as where R nd is a permutation matrix fulfilling R nd (X 1 , X 2 , . . . , X nd ) = (X 1 , X d+1 , . . . , X (n−1)d+1 , X 2 , X d+2 , . . . , X nd ) and F n is the univariate Discrete Fourier Transform (DFT) matrix (see Section 10.1 in Brockwell and Davis (1991)). We need the following auxiliary results. Lemma A.3 is a classic result that links the determinant of the time-domain Block Toeplitz covariance matrix with the integrated log spectral density. Lemma A.4 asymptotically bounds the norm of I nd − Γ nd Γ − nd and Lemma A.5 establishes rates when approximating the Gaussian frequency domain covariance matrix by a block circulant matrix (see e.g. Gray (2006) for an overview of this topic), which is closely related to the covariance matrix under Whittle's likelihood. Proof. The proof for the case d = 1 can be found in Lemma A1.4 in Dzhaparidze and Kotz (2012) and the proof for the case d > 1 follows along the lines.
We will employ the following notational convention: For a matrix A ∈ R nd×nd , denote the (disjoint) d × d blocks of A by A(i, j) ∈ R d×d for i, j = 1, . . . , n. Furthermore, for two positive sequences (a n ), (b n ) we write a n b n if there exists a constant c > 0 such that a n ≤ cb n for all n.
Lemma A.5. Let f fulfill Assumption 4.1 and H nd := F nd Γ nd F T nd − D nd with D nd from (A.7). Then there exists a constant C > 0 such that H nd (i, j) ≤ Cn −1 holds for i, j = 1, . . . , n and all n.
Proof. Consider the symmetric block circulant matrix Γ • nd ∈ R nd×nd defined as Using the representation of the block com- by Assumption 4.1. This shows G nd (i, j) n −1 uniformly in i, j as n → ∞. Now consider f n (ω) = 1 2π |h|≤ n/2 Γ(h) exp(−ihω) and the corresponding block diagonal matrix D nd [f n ] ∈ R nd×nd of f n as in (A.7). Using Proposition 4.5.1 in Brockwell and Davis (1991), a few elementary calculations show for 0 ≤ ω ≤ π by Assumption 4.1 and hence D nd [f n ](i, j)−D nd [f ](i, j) n −1 uniformly in i, j as n → ∞, the claim follows. Now we can present the proof of the contiguity result.
Proof of Theorem 4.2. It suffices to show thatP n W andP n from (A.6) and (A.8) are mutually contiguous. To this end, following the arguments from the proof of the univariate case in Choudhuri et al. (2004b), it suffices to show that the sequence of random variables has uniformly bounded mean and variance under bothP n W andP n . Using the result from Lemma A.3, the boundedness of log |Γ nd | − log |D nd | as n → ∞ follows with the same argument as for the univariate case, see Choudhuri et al. (2004b). LettingΛ n :=Z T (Γ −1 nd − D −1 nd )Z withΓ nd = F nd Γ nd F T nd , it remains to show that each of the following sequences is bounded: Let H nd be defined as in Lemma A.5. Using | tr(AB)| ≤ A 2 B with A 2 denoting the largest singular value, we get by Lemma A.5 and Assumption 4.1. A similar calculation that follows the arguments in Choudhuri et al. (2004b) shows VarP nΛn = O(1). For the mean and variance under Whittle's likelihood, consider a time series with spectral density matrix f −1 . From Lemma 13.3.2 in Gröchenig (2013), we get that f −1 also fulfills Assumption 4.1. Using Γ −1 nd 2 ≤ max 0≤ω≤π f (ω) −1 2 = O(1) from Lemma 2.1 in Hannan and Wahlberg (1989)

A.3. Proofs of Section 4.2
In this Section we will present the proof of Theorem 4.5. It relies on the contiguity result from Corollary 4.3. The proof technique consists of two main parts: prior positivity of neighborhoods and existence of exponentially powerful tests, both which are discussed in the following. Detailed proofs of all results in this section can be found in Meier (2018), Section 7.1.
Proof of Theorem 4.5. We will apply a general consistency theorem for noniid observations from Choudhuri et al. (2004a), see Theorem A.1 in Choudhuri et al. (2004a). The prior KL support and testability assumptions are verified in Lemma A.7 and Lemma A.11 below. It thus remains to bound the prior mass of the sieve complement: P 0,τ (Θ c n ) k>kn p(k), which is bounded from above by exp(−cn) for a constant c > 0 by Assumption 4.4 concluding the proof.

A.3.1. Prior positivity of neighborhoods
We start our considerations with prior posivity of uniform neighborhoods. The maximum Frobenius norm is defined as f F,∞ := max 0≤ω≤π f (ω) .
Lemma A.6. Let the assumptions of Theorem 4.5 be fulfilled. Then with B ε : Proof. The proof is analogously to the proof of (B.1) in Choudhuri et al. (2004a), using the insight from Lemma A.1, observing that for ε small enough it holds B ε ⊂ C 0,τ .
The result from Lemma A.6 can be used to show prior positivity of Kullback Leibler (KL) neighborhoods, as summarized in the following Lemma. Recall the definition of the Fourier coefficientsZ 1 , . . . ,Z N with N = n/2 − 1 from (3.2) and their respective probability densities p j,N from (4.4) under Whittle's likelihood P n W . Consider the KL terms where E f0 and Var f0 denote mean and variance under P n W (·|f 0 ). Lemma A.7. Under the assumptions of Theorem 4.5, it holds Proof. Using λ 1 (A+B) ≥ λ 1 (A)− B for A ∈ S + d and B ∈ S d we get λ 1 (f (ω)) ≥ λ 1 (f 0 (ω)) − f − f 0 , which is larger or equal to b 0 /2 for all f ∈ B ε and all ε small enough. Consider Q(ω) := f (ω) −1/2 f 0 (ω)f (ω) −1/2 , with A 1/2 denoting the (unique) Hpd matrix square root of A ∈ S + d .

A.3.2. Existence of tests
. Let us consider the following sieve (Θ n ) of the parameter space Θ = C 0,τ : where δ > 0 will be specified later. A few calculations using the equivalence of matrix norms reveal 14) The following result quantifies the metric entropy of Θ n in terms of ε-covering numbers (see Appendix C in Ghosal and van der Vaart (2017)). It is similar in spirit to Lemma B.4 in Choudhuri et al. (2004a) and will be needed later for the construction of uniformly exponentially powerful tests from tests against fixed alternatives, where the alternative set is covered by small balls and the number of balls is controlled to ensure test consistency.
Lemma A.8. The ε-covering number of Θ n with respect to · F,∞ is bounded by log N (ε, Θ n , · F,∞ ) k n log k n + log τ ε .
The following result is related to Lemma B.2 in Choudhuri et al. (2004a). It translates the integral condition f (ω) − f 0 (ω) dω > ε from the consistency complement sets to a suitable pointwise condition that holds on sufficiently many Fourier frequencies. This allows for testing at these frequencies. Denote by Ω n := {ω 0 , . . . , ω N , ω n/2 } the set of Fourier frequencies, with ω n/2 being included if and only if n is even.
Lemma A.9. Let τ > 0 and f 0 ∈ C 0,τ . Then there exists k 0 ∈ N such that for every k ≥ k 0 and every polynomial f ∈ C 0,τ of degree k with and bound the cardinality of A ∩ Ω n from below. Using A 1/2 BA 1/2 2 ≥ λ 1 (A) tr(BAB) = λ 1 (A) tr(ABB) ≥ λ 1 (A) 2 B 2 for A ∈ S + d and B ∈ S d (see Marshall et al. (2011), p.341), it follows that dτε is a subset of A and we will continue to bound the cardinality of B ∩Ω n from below. For k ≥ k 0 , there exists a polynomialf 0 ∈ C 0,τ of degree k with max ) 2 ) being a polynomial of degree at most 2k, it follows with the same argument as in the proof of Lemma B2 in Choudhuri et al. (2004a) that the cardinality of C ∩ Ω n is at least nε/(8πτ ) − 4k.
Based on these findings, we can now construct exponentially powerful tests. We start with tests against fixed alternatives.
A standard large deviation argument (as e.g. in Lemma A2 in Kirch et al. (2018) and another large deviation argument concludes the proof.
The previous result can now be used to construct uniform tests.
Lemma A.11. Let the assumptions of Theorem 4.5 be fulfilled. Let ε > 0 and U ε := {f : Then there exists a test ϕ n fulfilling for a constant c > 0, with E f denoting the mean under P n W (·|f ). Proof. Let f ∈ U c ε ∩ Θ n . We will show that there exists a test ϕ n,f fulfilling with positive constants c 0 , c 1 not depending on f . The test ϕ n as claimed can then be constructed by the standard approach (see e.g. Section B.2 in Choudhuri et al. (2004a) or Section 6.4 in Ghosal and van der Vaart (2017)) of covering the alternative set U c ε ∩ Θ n with small balls in conjunction with the bound for the covering number from Lemma A.8.

A.4. Proofs of Section 4.3
Proof of Theorem 4.6. The proof is analogous to the proof of Theorem 4.5 in Section A.3, replacing Lemma A.6 with Lemma A.12 below and Lemma A.9 with Lemma A.13 below.
Lemma A.12. Let the assumptions of Theorem 4.6 be fulfilled. For ε > 0, Proof. It suffices to show P (B ε ) > 0 for ε > 0 small enough, where P denotes the Bernstein-Hpd-Gamma prior (without restriction). Consider the Bernstein polynomial approximation of f 0 of degree k. Since f 0 is continuously differentiable by Assumption 4.1, it follows that the derivative of f k converges uniformly to the one of f 0 as k → ∞ (see Sections 1.4 and 1.8 in Lorentz (2012)). Hence, there exists k 0 such that it holds f 0 − f k L < ε/2 for all k ≥ k 0 . This yields Using b j,k L ≤ 2k(k − 1) + k ≤ 3k 2 (this follows from basic properties of Bernstein polynomials, see Section 1.4 in Lorentz (2012)) this leads to with ε k = ε/(6k 3 ). The right hand side of (A.19) can be shown to be positive with the same arguments as in the proof of (B.1) in Choudhuri et al. (2004a).
Lemma A.13. Let f 0 ∈C b1 be continuously differentiable. Then there exists k 0 ∈ N such that for every k ≥ k 0 and every f ∈C τ for τ < b 1 with components being polynomials of degree k, it holds true that and Ω n denotes the set of Fourier frequencies.

A.5. Proofs of Section 4.4
The proof of Theorem 4.8 relies on a general contraction rate result from Ghosal and Van Der Vaart (2007) and is again split in two parts: Prior mass of neighborhoods (Lemma A.15) and and sieve entropy (Lemma A.17). For detailed proofs of all results in this section, see Section 7.2 in Meier (2018).
Proof of Theorem 4.8. Following the arguments of Section 9.5.2 in Ghosal and van der Vaart (2017), we use Lemma A.15 below to find P τ0,τ1 (Θ c n ) P τ0,τ1 (B N,2 (f 0 , ε n )) = o(exp(−2nε 2 n )). (A.20) The result now follows by an application of Theorem 1 and Lemma 1 in Ghosal and Van Der Vaart (2007), noting that the assumptions thereof are fulfilled by Lemma A.15 and Lemma A.17 below, Lemma 2 in Ghosal and Van Der Vaart (2007) and (A.20).

Prior mass of neighborhoods
Denote by F 0 (A) := A f 0 (ω)dω ∈S + d for A ⊂ [0, π] the spectral measure corresponding to spectral density f 0 . Recall that {I j,k : j = 1, . . . , k} denotes the equidistant partition of [0, π] of size k. In the setting of Theorem 4.8, we will consider the following sieve on C τ0,τ1 : where ρ > 0 will be specified later and S +k d as defined before (A.13). The following result quantifies the probability mass that Φ puts in neighborhoods of F 0 .
Proof. Using 1/k n ≥ ε n for n large enough, we find where Theorem 2.2 (b) was used in the last step. The assertion now follows from Lemma A.2, noting that with κ α (I j,kn ) = exp(−cC α (I j,kn )) k n , C α (I j,kn ) = I j,kn α(x,S + d )dx, for j = 1, . . . , k n it holds kn j=1 C α,j,kn = π 0 α(x,S + d )dx < ∞ as well as kn j=1 κ α (I j,kn ) ≥ c/k kn n ≥ c exp(k n log ε n ). The previous result can be used to quantify the probability mass that the Bernstein-Hpd-Gamma prior puts in balls around f 0 . Recall the definition of the KL terms from (A.9) and consider the KL neighborhoods B N,2 (f 0 , ε) := {f ∈ Θ n : K N (f 0 , f ) < ε 2 , V N (f 0 , f ) < ε 2 } with Θ n from (A.21).
Lemma A.15. Let the assumptions of Theorem 4.8 be fulfilled. Then there exists a constant c > 0 such that P τ0,τ1 (B N,2 (f 0 , ε n )) ≥ exp(−cnε 2 n ). Proof. The proof is analogous to Section 9.5.2 in Ghosal and van der Vaart (2017), using the result from Lemma A.14.

Sieve entropy
The next aim is to derive an upper bound for the ε-covering number of Θ n in the Hellinger topology. We will use the bound in the maximum Frobenius topology from Lemma A.8 for this purpose. The following result established the needed link between the Hellinger and Frobenius topology.
The previous insight can be used to bound the root average squared Hellinger distance from (4.3) for spectral densities with uniformly bounded Lipschitz constants from below.

Appendix B: MCMC Algorithm for Posterior Computation
To draw random samples from the posterior distribution of f , we employ a Gibbs sampler (see Section 6.3.2 in Christensen et al. (2011)) based on the approximate parametrization from (3.4) for f . We assume that the Hpd Gamma process prior on Φ fulfills the assumptions of Theorem 2.3. Additionally, we assume that the measure α on X ×S + d has a Lebesgue density, which will be denoted by g. The corresponding probability density will be denotes by g * = g/C α . From Theorem 2.3, we obtain that the prior for Θ f from (3.4) is given by where J T denotes the Jacobian of the transformation r := (r 1 , . . . , r L ) → (v 1 , . . . , v L ) =: v as in Theorem 2.3. It can be shown that log |J T (r)| = L log C α − L l=1 (log r l + β(x l , U l )r l ) , see Section 5.2.1 in Meier (2018) for the details. The matrices U 1 , . . . , U L ∈ S + d are parametrized in terms of hyperspherical coordinates as U l =T ϕ l for l = 1, . . . , L, whereT is a one-to-one transformation and ϕ l is a (d 2 −1)-dimensional vector with components in (0, π/2) and (0, π). The Jacobian determinant ofT is available analytically, see (27)-(31) in Mittelbach et al. (2012) or Section 3.4.1 in Meier (2018).
Since the full conditionals are not available in a closed form, we employ Metropolis-Hastings (MH) steps (see Section 6.3.3 in Christensen et al. (2011)) for the components of Θ f , where the starting values and proposal distributions will be discussed in the following.
For the sake of computational speed-up, the values of k can be restricted to a finite set {1, . . . , k max }, where k max is some large integer. This has the computational advantage that the Bernstein polynomial basis functions can be precomputed and stored for usage in the iterations of the sampling algorithm. A feasible choice of k max depends on the data at hand and can be determined by preliminary pilot runs. For all examples considered in this work, a value of k max = 500 was large enough. To draw samples from the full conditional of k, we choose some large integer as starting value. Although samples from the full conditional could be drawn from the finite set {1, . . . , k max }, we instead use a MH step to avoid the computationally expensive task of k max conditional posterior evaluations. The proposal is chosen as a random walk scheme with discretized Cauchy distribution increments.
The radial parts r 1 , . . . , r L are updated one at a time by individual MH steps with normal random walk proposals on the log scale. To elaborate, given the value r The x l 's in (0, π) are also updated by individual MH steps, with random walk proposals, where the increments are Unif([−δ l , δ l ]) distributed with δ l = πl/(l+2 √ n). The values of δ l are as suggested in Choudhuri et al. (2004a) for the location parameters of a Dirichlet process. We also employ the author's strategy of circular proposals whenever a proposal lands outside the interval [0, π]. The starting values are drawn iid from the Unif([0, π]) distribution.
The matrices U l are sampled one at a time with a MH step for ϕ l . The components of ϕ l = (ϕ l,1 , . . . , ϕ l,i 2 −1 ) are sampled blockwise, where the proposals are drawn independently from a random walk scheme with Unif([−a jδl , a jδl ]) increments, where a j denotes the length of the interval of definition of ϕ l,j andδ l > 0 is a scaling parameter that is determined adaptively durin burn-in (again with target acceptance rate of 0.44).