code level $t$. Let ${\boldsymbol \beta}:=({\boldsymbol \beta}_1^\top, \ldots, {\boldsymbol \beta}_s^\top)^\top$, ${\boldsymbol \gamma}:=(\gamma_1, \ldots, \gamma_{s-1})^\top$, ${\boldsymbol \sigma}^2:=(\sigma^2_1, \ldots, \sigma^2_s)^\top$, ${\boldsymbol \phi}:=({\boldsymbol \phi}_1^\top, \ldots, {\boldsymbol \phi}_s^\top)^\top$, and ${\mathbf{y}}=({\mathbf{y}}_{1}^{\top},\ldots,{\mathbf{y}}_{s}^{\top})^{\top}$. Then the marginal likelihood is $$\begin{aligned}
 \label{eqn: sampling dist}
\begin{split}L({{\mathbf{y}}}\mid{\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2},{\boldsymbol \phi}) & =\pi({{\mathbf{y}}}_{1}\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1})\prod_{t=2}^{s}\pi({{\mathbf{y}}}_{t}\mid{{\mathbf{y}}}_{t-1},\gamma_{t-1},{\boldsymbol \beta}_{t},\sigma_{t}^{2},{\boldsymbol \phi}_{t})\end{split},\end{aligned}$$ where $$\begin{aligned}
\begin{split}\pi({{\mathbf{y}}}_{1}\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1}) & =\mathcal{N}({{\mathbf{H}}}_{1}{\boldsymbol \beta}_{1},\,\sigma_{1}^{2}{{\mathbf{R}}}_1),\\
\pi({{\mathbf{y}}}_{t}\mid {{\mathbf{y}}}_{t-1},\gamma_{t-1},{\boldsymbol \beta}_{t},\sigma_{t}^{2},{\boldsymbol \phi}_{t}) & =\mathcal{N}({{\mathbf{H}}}_{t}{\boldsymbol \beta}_{t} + {W}_{t-1}\gamma_{t-1},\,\sigma_{t}^{2} {{\mathbf{R}}}_t),
\end{split}
\label{eqn: conditional sampling distributions}\end{aligned}$$ with ${{\mathbf{H}}}_t:={\mathbf{h}}_t({{{\cal X}}}_t)$ and ${W}_{t-1}:=y_{t-1}({{{\cal X}}}_t)$, where $y_{t-1}(A):=[y_{t-1}(x),x\in A]$ is a vector of output values over inputs in $A$. This sampling distribution provides a convenient form to perform closed-form likelihood-based inference.

The Cokriging Predictor and Cokriging Variance {#sec: prediction in univariate model}
==============================================

For any new input ${\mathbf{x}}_{0}\in {{\cal X}}$, the goal is to make prediction for $y_s({\mathbf{x}}_0)$ based upon the code output ${\mathbf{y}}$. In [@Kennedy2000], a closed-form predictive distribution is derived for $\pi(y_s({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi})$, which only accounts for the uncertainty due to estimation of ${\boldsymbol \beta}$ and has $O((\sum_{t=1} n_t)^3)$ computational cost. In [@Qian2008], a closed-form predictive distribution is only given for $\pi(y_s({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi})$, which also has $O((\sum_{t=1}^s n_t)^3)$ computational cost. To account for 