& \propto\frac{1}{\sigma_{1}^{2}},\\
\pi^{R}({\boldsymbol \beta}_{t},\gamma_{t-1},\sigma_{t}^{2}) & \propto\frac{1}{\sigma_{t}^{2}},\,t=2,\ldots,s.
\end{split}\end{aligned}$$

The following lemma gives the predictive distribution of ${\mathbf{y}}({\mathbf{x}}_0):=(y_1({\mathbf{x}}_0), \ldots, y_s({\mathbf{x}}_0))^\top$ given ${\mathbf{y}}$ and ${\boldsymbol \phi}$.

\[lem: conditional predictive dists\] According to the cokriging model , we have $$\begin{aligned}
\label{eqn: joint predictive dist under noninformative priors}
\begin{split}
\pi({\mathbf{y}}({\mathbf{x}}_0)\mid {\mathbf{y}},{\boldsymbol \phi}) &=\pi(y_1({\mathbf{x}}_0)\mid {\mathbf{y}}_{1},{\boldsymbol \phi}_{1})\prod_{t=2}^{s-1}\pi(y_t({\mathbf{x}}_0)\mid {{\mathbf{y}}}_{t-1}, y_{t-1}({\mathbf{x}}_0), {\mathbf{y}}_{t},{\boldsymbol \phi}_{t})\\
&\times \pi(y_s({\mathbf{x}}_0) \mid y_{s-1}({\mathbf{x}}_0), {\mathbf{y}}_s, {\boldsymbol \phi}_s),
\end{split}\end{aligned}$$ where conditional distributions on the right-hand side are Student $t$-distributions $t_{n_t-q_t}(\mu_t({\mathbf{x}}_0), $ $ \Sigma_t({\mathbf{x}}_0))$ given by $$\begin{aligned}
\begin{split}
\mu_{t}({\mathbf{x}}_0) & :={\mathbf{X}}_t^\top({\mathbf{x}}_0)\hat{{\mathbf{b}}}_t+{\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t),\\
\Sigma_{t}({\mathbf{x}}_0) & :=\hat{\sigma}^2_t c_t^*,
\end{split}\end{aligned}$$ with $$\begin{aligned}
\begin{split}
\hat{\sigma}^2_t & :=({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t)^{\top}{\mathbf{R}}_t^{-1}({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t) / (n_t-q_t),\\
c_t^* & :=r({\mathbf{x}}_0,{\mathbf{x}}_0|{\boldsymbol \phi}_{t})-{\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}{\mathbf{r}}_t({\mathbf{x}}_0) \\
&+ [{\mathbf{X}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)]^{\top}({\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{X}}_{t})^{-1}[{\mathbf{X}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)],
\end{split}\end{aligned}$$ where $\hat{{\mathbf{b}}}_t:=(\hat{{\boldsymbol \beta}}^\top, \hat{\gamma}_{t-1})^\top = ({\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t)^{-1}{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{y}}_t$, ${\mathbf{r}}_t({\mathbf{x}}_0):=r({{{\cal X}}}_t, {\mathbf{x}}_0| {\boldsymbol \phi}_t)$, ${\mathbf{X}}_1:={\mathbf{H}}_1$, and ${\mathbf{X}}_t:=[{{\mathbf{H}}}_{t}, y_{t-1}({{\cal X}}_t)]$ for $t>1$. $q_t$ is the number of columns in ${\mathbf{X}}_t$.

See Appendix \[app: conditional predictive dist\].

This result shows that one can generate a random sample from the predictive distribution $\pi({\mathbf{y}}({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \phi})$ by sequentially sampling from a collection of conditional distributions. Notice that samples are obtained across all the code levels. The computation associated with each conditional distribution only requires $O(n_t^3)$ flops for $t=1, \ldots, s$. Before stating the next theorem that provides 