of generality, we focus on the simple form, i.e., ${\mathbf{k}}_{t-1}$ is assumed to be 1 and ${\boldsymbol \zeta}_{t-1}$ is assumed to be a scalar parameter.

To account for uncertainties in the unknown functions $y_{1}(\cdot)$ and $\delta_{t}(\cdot)$, Gaussian process priors can be assigned: $$\begin{aligned}
\begin{split}y_{1}(\cdot)\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1} & \sim\mathcal{GP}({\mathbf{h}}_{1}(\cdot)^{\top}{\boldsymbol \beta}_{1},\,\sigma_{1}^{2}r(\cdot,\cdot|{\boldsymbol \phi}_{1})),\\
\delta_{t}(\cdot) & \sim\mathcal{GP}({\mathbf{h}}_{t}(\cdot)^{\top}{\boldsymbol \beta}_{t},\,\sigma_{t}^{2}r(\cdot,\cdot|{\boldsymbol \phi}_{t})),
\end{split}
\label{eqn: co-kriging model}\end{aligned}$$ for $t=2,\ldots,s$, $i=1,\ldots,n_{t}$, where $r(\cdot,\cdot|{\boldsymbol \phi}_{t})$ is a correlation function with correlation parameters ${\boldsymbol \phi}_{t}$. A popular form is to choose the product correlations with the power exponential family and the Matérn family. ${\mathbf{h}}_{t}(\cdot)$ is a vector of (fixed) basis functions and ${\boldsymbol \beta}_{t}$ is a vector of unknown coefficients at code level $t$. $\sigma^2_t$ is the variance parameter.

The cokriging model in  and  has been used to model computer model output at different fidelity levels in previous works [@Kennedy2000; @Gratiet2013] with hierarchically nested design, i.e., ${{\cal X}}_{t}\subset{{\cal X}}_{t-1}$. In [@Qian2008], a measurement-error process is incorporated in to link observations with computer model outputs at two fidelity levels. The assumption of hierarchically nested design in these works allows for closed-form likelihood-based inference in discussions that follow. Let ${{\mathbf{y}}}_{t}$ be a vector of output values at all inputs in ${{{\cal X}}}_{t}$ at 