the posterior predictive distribution $\pi({\mathbf{y}}({\mathbf{x}}_0)\mid {\mathbf{y}}, \hat{{\boldsymbol \phi}})$. Notice that there is no need to estimate other model parameters $\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2\}$. If desired, these model parameters can be estimated based on the posterior distribution $\pi({\boldsymbol \beta}, {\boldsymbol \gamma}\mid {\mathbf{y}}, \hat{{\boldsymbol \phi}})$ and $\pi({\boldsymbol \sigma}^2\mid {\mathbf{y}}, \hat{{\boldsymbol \phi}})$. The detailed procedures to estimate these parameters are given in Appendix \[app: parameter estimation\].

Numerical Illustration {#sec: numerical demonstration}
======================

The main goal of the numerical illustration is to demonstrate the predictive performance of the autoregressive cokriging model with objective priors developed in previous sections. In addition, we also include the jointly robust prior [@Gu2019] in the comparison, since the jointly robust prior mimics the behavior of reference priors for Gaussian process models and it is a proper prior that allows fast computation. The form of the jointly robust prior is $$\begin{aligned}
\pi^{JR}(B_1, \ldots, B_d) = C\left( \sum_{i=1}^d C_i B_i  \right)^{a_0} \exp\left\{ - b_0 \left( \sum_{i=1}^d C_i B_i  \right)  \right\},\end{aligned}$$ where $B_i$’s are inverse range parameters; $C$ is a normalizing constant; $a_0>-(d+1)$, $b_0>0$ and $C_i=n^{-1/d}|x_i^{max} - x_i^{min} |$are hyperparameters. $a_0$ is a parameter controlling the polynomial penalty to avoid singular correlation matrix and $b_0$ is a parameter controlling the 