The resulting cokriging model is often referred to as an autoregressive cokriging model. Several extensions of this autoregressive cokriging model have been proposed with increased model flexibility and Bayesian inference approaches [@Qian2008; @Gratiet2013; @Gratiet2014].

A common feature found in these works is that the predictive distribution for the high-fidelity code given both low and high fidelity output at a set of inputs as well as correlation parameters (or range parameters) requires numerical integration, which leads to intensive computations. Indeed, the predictive distribution can be available in a closed-form when other model parameters such as regression parameters and variance parameters are conditioned upon, but this leads to under-estimation of uncertainties associated with predictions. The uncertainty analysis about the quantity of interest, often a transformation of predictors, will hence suffer severely from this artifact. To avoid this, a recursive predictive formula is derived such that predictive distributions are only conditioned upon the code output and correlation parameters. Closed-form predictive means and predictive variances are also derived at each code level that can be computed without the need of Monte Carlo approximations as in [@Qian2008; @Gratiet2013]. This formula explicitly accounts for the uncertainty due to the estimation of regression parameters, scale discrepancy parameter 