uncertainty due to estimation of model parameters ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi}$, Monte Carlo approximation is required. In [@Gratiet2013], a closed-form predictive distribution is also only given for $\pi(y_s({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi})$, and Monte Carlo approximation is used to account for uncertainty due to estimation of other model parameters. [@Gratiet2013] also develops an iterative formula to invert the $(\sum_{t=1}^s n_t)\times (\sum_{t=1}^s n_t)$ correlation matrix of code output at all levels, which reduces computation cost to $O(\sum_{t=1}^s n_t^3)$.

In what follows, we give a new way to derive closed-form predictive distributions for $\pi(y_t({\mathbf{s}}_0) \mid {\mathbf{y}}, {\boldsymbol \phi}), t=1, \ldots, s$ that not only explicitly account for the uncertainty due to estimation of ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$ but also has $O(\sum_{t=1}^s n_t^3)$ computational cost. The formula for these predictive distributions is derived based upon the idea that the new input ${\mathbf{x}}_0$ is added to each ${{{\cal X}}}_t$ such that a hierarchically nested design can be obtained.

To deal with these unknown parameters ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$, the following standard reference priors are used for the location-scale parameters: ${\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2}$: $$\begin{aligned}
 \label{eqn: location-scale prior}
\begin{split}
\pi^{R}({\boldsymbol \beta}_{1},\sigma_{1}^{2}) 