exponential penalty to avoid diagonal correlation matrix. $n$ here refers to the number of model runs; $x_i^{max}, x_i^{min}$ refer to the maximum and minimum of input parameter $x_i$. [@Gu2019] recommends default settings for these parameters: $a_0=0.5-d$ and $b_0=1$. However, it was pointed out in [@Gu2019] that the choice of $a_0$ is an open problem and is problem-specific. In the following numerical study, we fix $b_0$ at 1, and tune the parameter $a_0$ to achieve comparable results. For the autoregressive cokriging model, independent jointly robust priors are assumed for correlation parameters across different levels of fidelity. In the following numerical comparison, the proposed cokriging predictors and cokriging variances explicitly take into account the uncertainties in estimating ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$, while the closed-form predictive formulas in [@Gratiet2013] do not take into account the uncertainties in estimating ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$, and ${\boldsymbol \phi}$. As numerical examples in [@Gratiet2013] indicate that the approach in [@Gratiet2013] performs better than the approaches in [@Kennedy2000; @Qian2008], which is computationally much more expensive, it is more interesting to compare the proposed inference approach and the approach in [@Gratiet2013]. In addition, the proposed inference approach is very similar to the approach in [@Gratiet2013] 