How to create a multivariate Brownian Bridge? It is known, that a standard multivariate Brownian bridge $ y(\mathbf u) $ is a centered Gaussian process with covariance function
$$ \mathbb E(y(\mathbf u) y(\mathbf v)) = \prod_{j=1}^d (u_j \wedge v_j) - \prod_{j=1}^d u_j v_j $$
I am not sure about how to constuct such a multivariate Brownian bridge.
My first thought was to start somehow with a univariate Brownian bridge. I have found information about that and even a package in R that can do this, but only for the univariate Brownian bridge.
I found this, but as I understand it, what has been done there is not a standard multivariate Brownian bridge as defined above or e.g. in this paper.
I would appreciate any hints and support. 
 A: As you already pointed out in the comments, the question reduces to simulating a Brownian sheet. This can be done by generalizing simulation of Brownian motion in a straightforward way.
To simulating the Brownian motion, one can take an i.i.d. mean-0 variance-1 time series $W_i$, $i = 1, 2, \cdots$, and construct the normalized partial sum process
$$
X_n(t) = \frac{1}{\sqrt{n}} \sum_{i = 1}^{[nt]} W_i.
$$ 
As $n \rightarrow \infty$, $X_n$ convergence weakly (in the sense of Borel probability measures on a metric space) to the standard Brownian $B$ on the Skorohod space $D[0,1]$.
The i.i.d. with finite second moment case is the simplest way to simulate. The mathematical result (Functional Central Limit Theorem/Donsker's Theorem/Invariance Principle) holds in much greater generality.
Now to simulating the (say, two-dimensional) Brownian sheet, takes i.i.d. mean-0 variance-1 array $W_{ij}$, $i,j = 1, 2, \cdots$, and construct the normalized partial sum process
$$
X_n(t_1, t_2) = \frac{1}{ n } \sum_{1 \leq i \leq [nt_1] , 1 \leq j \leq [nt_2]} W_{ij}.
$$ 
As $n \rightarrow \infty$, $X_n$ convergence weakly to the standard Brownian sheet on the Skorohod space $D([0,1]^2)$ on the unit square.
(The proof is a standard weak convergence argument:


*

*Convergence of finite dimensional distribution follows from the Levy-Lindeberg CLT.

*Tightness on $D([0,1]^2)$ follows from a sufficient moment condition that holds trivially in the i.i.d. finite second moment case---see, e.g. Bickel and Wichura (1971). 
)
Then, by the continuous mapping theorem
$$
X_n(t_1, t_2) - \prod_{j=1}^2 t_j X_n(t_1, t_2)
$$
converges weakly to the two-dimensional Brownian bridge.
