the predictive mean and predictive variance in  can be computed exactly with computational cost $O(\sum_{t=1}^s n_t^3)$. As a byproduct, predictions at code levels from $t=1$ to $t=s-1$ are obtained automatically. For $t=1$, the predictive mean and predictive variance in the cokriging model are exactly the universal kriging predictor and universal kriging variance in a kriging model. For $t>1$, the cokriging predictor is a sum of a kriging predictor and an additional constant, i.e., $\hat{y}_t({\mathbf{x}}_0)={\mathbf{h}}_t^\top({\mathbf{x}}_0)\hat{{\boldsymbol \beta}}_t + {\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}({\mathbf{y}}_t - {\mathbf{H}}_t\hat{{\boldsymbol \beta}}_t) + [\hat{y}_{t-1}({\mathbf{x}}_0) + {\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}W_{t-1}]\hat{\gamma}_{t-1}$. The computational cost for both parameter estimation and prediction in an $s$-level cokriging model is equivalent to the one in $s$ independent kriging models. This predictive distribution allows us to explicitly integrate out models parameters ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$ except the range parameters ${\boldsymbol \phi}$, and hence carries several advantages over the predictive distribution given in [@Gratiet2013]. Specifically, the predictive distribution in [@Gratiet2013] is a normal distribution conditioned on all model parameters $\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi}\}$, and intensive Monte Carlo approximations are required to take into account the uncertainty due to estimation of $\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2\}$ in order to derive the predictive distribution $\pi(y_s({\mathbf{x}}_0)\mid {\mathbf{y}}, {\boldsymbol \phi})$.

