---
abstract: 'Autoregressive cokriging models have been widely used to emulate multiple computer models with different levels of fidelity. The dependence structures are modeled via Gaussian processes at each level of fidelity, where covariance structures are often parameterized up to a few parameters. The predictive distributions typically require intensive Monte Carlo approximations in previous works. This article derives new closed-form formulas to compute the means and variances of predictive distributions in autoregressive cokriging models that only depend on correlation parameters. For parameter estimation, we consider objective Bayesian analysis of such autoregressive cokriging models. We show that common choices of prior distributions, such as the constant prior and inverse correlation prior, typically lead to improper posteriors. We also develop several objective priors such as the independent reference prior and the independent Jeffreys prior that are shown to yield proper posterior distributions. This development is illustrated with a borehole function in an eight-dimensional input space and applied to an engineering application in a six-dimensional input space. codes are available in the Supplementary Material to reproduce the numerical results.'
author:
- 'Pulong Ma[^1]'
bibliography:
- 'references.bib'
title: Objective Bayesian Analysis of a Cokriging Model for Hierarchical Multifidelity Codes
---

Introduction
============

Complex computer codes have been widely used to solve mathematical models that represent real-world processes in virtually every field of science and engineering. They are often referred to as *simulators* in Uncertainty Quantification (UQ) and computer experiments [@OHagan2006; @Santner2018]. In practice, computer codes can be too timing-consuming to be used for adequately addressing UQ tasks. To overcome this bottleneck, Gaussian processes have been widely used as surrogate models to approximate simulators due to its computational advantages and attractive theoretical properties [@Sacks1989].

In real applications, computer codes can be run at different levels of accuracy due to sophistication of physics incorporated in mathematical models, accuracy of numerical solvers and resolutions of meshes; see [@Peherstorfer2018] for formal definition of multifidelity models. Several works have been proposed to combine output from computer codes at different fidelity levels based on a well-known geostatistical method called *cokriging*; see Chapter 3 of [@Cressie1993]. The approach of cokriging to synthesizing multiple computer model outputs is originated in [@Kennedy2000], which is developed based upon a first order Markov assumption that given output from a low-fidelity code run at an input, no more information can be learnt about the high-fidelity code with output from the low-fidelity code at any other input. 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 and variance parameters.

Inference about model parameters has been approached in several different ways in autoregressive cokriging models with focus on empirical Bayesian approaches. A particular challenge in autoregressive cokriging models is to estimate correlation parameters. To tackle this issue, [@Kennedy2000] assume independent noninformative priors for all the parameters, and then carry out numerical maximization for the marginal likelihood functions after integrating out regression parameters with respect to scale discrepancy parameters, variance parameters and correlation parameters [@Kennedy2000]. In [@Qian2008], conjugate priors are assumed for regression parameters, scale discrepancy parameters and variance parameters. For correlation parameters, proper gamma priors are assumed. The correlation parameters are then estimated by maximizing the corresponding marginal posteriors. For the code at the first level, standard nonlinear optimization is performed; while for the code at the second level, the corresponding posterior does not have a closed form and its evaluation requires numerical integration. This optimization procedure is then recast into a stochastic programming problem. To alleviate computational difficulties in [@Kennedy2000; @Qian2008], [@Gratiet2013] develop an efficient joint Bayesian estimation approach with either non-information priors or informative priors for all the model parameters except the correlation parameters. Without further assuming prior distributions for correlation parameters, [@Gratiet2013] maximizes a concentrated restricted likelihood to obtain estimates for correlation parameters at each code level. However, whether the choices of priors in [@Qian2008; @Gratiet2013] will lead to good estimates is not discussed. We show that vague proper priors for correlation parameters in [@Qian2008] will lead to an improper posterior. Thus, the usage of such vague priors will not solve but hide the problem, see [@Berger2006] for detailed discussion and references therein. We also show that the concentrated restricted likelihood in [@Gratiet2013] with noninformative priors and informative priors (when chosen to be vague) can have nonrobust estimates in autoregressive cokriging models, where nonrobustness is defined to be the situation where the correlation matrix becomes either singular or identity when correlation parameters go to zero or infinity; see [@Gu2018] for detailed discussions.

This article has two primary objectives that are of interest from computational and theoretical perspectives. The first objective is the derivation of new formulas for predictive distributions of the code output at any level over a new input given code output and correlation parameters. Realizations can be simulated from predictive distributions based upon a set of conditional distributions and the predictive means and predictive variances can be computed exactly at any fidelity level in a computationally efficient way. The new closed-form predictive formulas will take into account uncertainties due to the estimation of the location and scale parameters. The second objective is the development of objective priors. The objective priors can be used as default priors when elicitation of prior information is challenging. It also enables more accurate uncertainty estimation in the predictive distribution than the typical maximum-likelihood based approaches with commonly-used noninformative priors.

Section \[sec: univariate model\] gives the general assumptions and reviews the autoregressive cokriging models. In Section \[sec: prediction in univariate model\], new closed-form expressions are derived for predictive distributions at all code levels conditioned on all code output and correlation parameters. The $s$-level cokriging model turns out to have the same computational cost as $s$ independent kriging models in both parameter estimation and prediction. Section \[sec: objective Bayes\] begins with discussions on commonly-chosen noninformative priors and proves that the resulting posteriors are improper with such noninformative priors. The objective priors including independent reference priors and independent Jeffreys priors are then developed and are shown to yield proper posteriors. Section \[sec: numerical demonstration\] gives several numerical examples. Section \[sec: discussion\] is concluded with further discussions.

The Autoregressive Cokriging Model {#sec: univariate model}
==================================

Suppose that we have $s$ levels of code $y_{1}(\cdot),\ldots,$ $y_{s}(\cdot)$, where the code $y_{t}(\cdot)$ is assumed to be more accurate than the code $y_{t-1}(\cdot)$ for $t=2,\ldots,s$. Let ${{\cal X}}$ be a compact subset of $\mathbb{R}^{d}$, which is assumed to be the input space of computer code. Further assume that the code $y_{t}(\cdot)$ is run at a set of input values denoted by ${{\cal X}}_{t}\subset{{\cal X}}$ for $t=1,\ldots,s$, where ${{\cal X}}_{t}$ is assumed to contain $n_{t}$ input values. Consider the following autoregressive model as in [@Kennedy2000; @Gratiet2013]: $$\begin{aligned}
 \label{eqn: AR}
y_{t}({\mathbf{x}})=\gamma_{t-1}y_{t-1}({\mathbf{x}})+\delta_{t}({\mathbf{x}}), \, {\mathbf{x}}\in {{\cal X}},\end{aligned}$$ for $t=2,\ldots,s$, where $y_{t-1}(\cdot)$ is an unknown function of input. $\delta_{t}(\cdot)$ is the unknown location discrepancy function representing the local adjustment from level $t-1$ to level $t$. $\gamma_{t-1}$ is the scale discrepancy representing the scale change from level $t-1$ to level $t$. Notice that currently $\gamma_{t-1}$ does not depend on input. A more general assumption is to take $\gamma_{t-1}(\cdot)$ to be a basis-function representation, i.e., $\gamma_{t-1}(\cdot)={\mathbf{k}}_{t-1}(\cdot)^{\top}{\boldsymbol \zeta}_{t-1}$ for $t=2,...,s$, where ${\mathbf{k}}_{t-1}(\cdot)$ is a vector of basis functions and ${\boldsymbol \zeta}_{t-1}$ is a vector of unknown coefficients with dimension $q_{\zeta}$. The development in this article is true for this general parameterization. Without loss of generality, we focus on the simple form, i.e., ${\mathbf{k}}_{t-1}$ is assumed to be 1 and ${\boldsymbol \zeta}_{t-1}$ is assumed to be a scalar parameter.

To account for uncertainties in the unknown functions $y_{1}(\cdot)$ and $\delta_{t}(\cdot)$, Gaussian process priors can be assigned: $$\begin{aligned}
\begin{split}y_{1}(\cdot)\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1} & \sim\mathcal{GP}({\mathbf{h}}_{1}(\cdot)^{\top}{\boldsymbol \beta}_{1},\,\sigma_{1}^{2}r(\cdot,\cdot|{\boldsymbol \phi}_{1})),\\
\delta_{t}(\cdot) & \sim\mathcal{GP}({\mathbf{h}}_{t}(\cdot)^{\top}{\boldsymbol \beta}_{t},\,\sigma_{t}^{2}r(\cdot,\cdot|{\boldsymbol \phi}_{t})),
\end{split}
\label{eqn: co-kriging model}\end{aligned}$$ for $t=2,\ldots,s$, $i=1,\ldots,n_{t}$, where $r(\cdot,\cdot|{\boldsymbol \phi}_{t})$ is a correlation function with correlation parameters ${\boldsymbol \phi}_{t}$. A popular form is to choose the product correlations with the power exponential family and the Matérn family. ${\mathbf{h}}_{t}(\cdot)$ is a vector of (fixed) basis functions and ${\boldsymbol \beta}_{t}$ is a vector of unknown coefficients at code level $t$. $\sigma^2_t$ is the variance parameter.

The cokriging model in  and  has been used to model computer model output at different fidelity levels in previous works [@Kennedy2000; @Gratiet2013] with hierarchically nested design, i.e., ${{\cal X}}_{t}\subset{{\cal X}}_{t-1}$. In [@Qian2008], a measurement-error process is incorporated in to link observations with computer model outputs at two fidelity levels. The assumption of hierarchically nested design in these works allows for closed-form likelihood-based inference in discussions that follow. Let ${{\mathbf{y}}}_{t}$ be a vector of output values at all inputs in ${{{\cal X}}}_{t}$ at code level $t$. Let ${\boldsymbol \beta}:=({\boldsymbol \beta}_1^\top, \ldots, {\boldsymbol \beta}_s^\top)^\top$, ${\boldsymbol \gamma}:=(\gamma_1, \ldots, \gamma_{s-1})^\top$, ${\boldsymbol \sigma}^2:=(\sigma^2_1, \ldots, \sigma^2_s)^\top$, ${\boldsymbol \phi}:=({\boldsymbol \phi}_1^\top, \ldots, {\boldsymbol \phi}_s^\top)^\top$, and ${\mathbf{y}}=({\mathbf{y}}_{1}^{\top},\ldots,{\mathbf{y}}_{s}^{\top})^{\top}$. Then the marginal likelihood is $$\begin{aligned}
 \label{eqn: sampling dist}
\begin{split}L({{\mathbf{y}}}\mid{\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2},{\boldsymbol \phi}) & =\pi({{\mathbf{y}}}_{1}\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1})\prod_{t=2}^{s}\pi({{\mathbf{y}}}_{t}\mid{{\mathbf{y}}}_{t-1},\gamma_{t-1},{\boldsymbol \beta}_{t},\sigma_{t}^{2},{\boldsymbol \phi}_{t})\end{split},\end{aligned}$$ where $$\begin{aligned}
\begin{split}\pi({{\mathbf{y}}}_{1}\mid{\boldsymbol \beta}_{1},\sigma_{1}^{2},{\boldsymbol \phi}_{1}) & =\mathcal{N}({{\mathbf{H}}}_{1}{\boldsymbol \beta}_{1},\,\sigma_{1}^{2}{{\mathbf{R}}}_1),\\
\pi({{\mathbf{y}}}_{t}\mid {{\mathbf{y}}}_{t-1},\gamma_{t-1},{\boldsymbol \beta}_{t},\sigma_{t}^{2},{\boldsymbol \phi}_{t}) & =\mathcal{N}({{\mathbf{H}}}_{t}{\boldsymbol \beta}_{t} + {W}_{t-1}\gamma_{t-1},\,\sigma_{t}^{2} {{\mathbf{R}}}_t),
\end{split}
\label{eqn: conditional sampling distributions}\end{aligned}$$ with ${{\mathbf{H}}}_t:={\mathbf{h}}_t({{{\cal X}}}_t)$ and ${W}_{t-1}:=y_{t-1}({{{\cal X}}}_t)$, where $y_{t-1}(A):=[y_{t-1}(x),x\in A]$ is a vector of output values over inputs in $A$. This sampling distribution provides a convenient form to perform closed-form likelihood-based inference.

The Cokriging Predictor and Cokriging Variance {#sec: prediction in univariate model}
==============================================

For any new input ${\mathbf{x}}_{0}\in {{\cal X}}$, the goal is to make prediction for $y_s({\mathbf{x}}_0)$ based upon the code output ${\mathbf{y}}$. In [@Kennedy2000], a closed-form predictive distribution is derived for $\pi(y_s({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi})$, which only accounts for the uncertainty due to estimation of ${\boldsymbol \beta}$ and has $O((\sum_{t=1} n_t)^3)$ computational cost. In [@Qian2008], a closed-form predictive distribution is only given for $\pi(y_s({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi})$, which also has $O((\sum_{t=1}^s n_t)^3)$ computational cost. To account for 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}) & \propto\frac{1}{\sigma_{1}^{2}},\\
\pi^{R}({\boldsymbol \beta}_{t},\gamma_{t-1},\sigma_{t}^{2}) & \propto\frac{1}{\sigma_{t}^{2}},\,t=2,\ldots,s.
\end{split}\end{aligned}$$

The following lemma gives the predictive distribution of ${\mathbf{y}}({\mathbf{x}}_0):=(y_1({\mathbf{x}}_0), \ldots, y_s({\mathbf{x}}_0))^\top$ given ${\mathbf{y}}$ and ${\boldsymbol \phi}$.

\[lem: conditional predictive dists\] According to the cokriging model , we have $$\begin{aligned}
\label{eqn: joint predictive dist under noninformative priors}
\begin{split}
\pi({\mathbf{y}}({\mathbf{x}}_0)\mid {\mathbf{y}},{\boldsymbol \phi}) &=\pi(y_1({\mathbf{x}}_0)\mid {\mathbf{y}}_{1},{\boldsymbol \phi}_{1})\prod_{t=2}^{s-1}\pi(y_t({\mathbf{x}}_0)\mid {{\mathbf{y}}}_{t-1}, y_{t-1}({\mathbf{x}}_0), {\mathbf{y}}_{t},{\boldsymbol \phi}_{t})\\
&\times \pi(y_s({\mathbf{x}}_0) \mid y_{s-1}({\mathbf{x}}_0), {\mathbf{y}}_s, {\boldsymbol \phi}_s),
\end{split}\end{aligned}$$ where conditional distributions on the right-hand side are Student $t$-distributions $t_{n_t-q_t}(\mu_t({\mathbf{x}}_0), $ $ \Sigma_t({\mathbf{x}}_0))$ given by $$\begin{aligned}
\begin{split}
\mu_{t}({\mathbf{x}}_0) & :={\mathbf{X}}_t^\top({\mathbf{x}}_0)\hat{{\mathbf{b}}}_t+{\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t),\\
\Sigma_{t}({\mathbf{x}}_0) & :=\hat{\sigma}^2_t c_t^*,
\end{split}\end{aligned}$$ with $$\begin{aligned}
\begin{split}
\hat{\sigma}^2_t & :=({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t)^{\top}{\mathbf{R}}_t^{-1}({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t) / (n_t-q_t),\\
c_t^* & :=r({\mathbf{x}}_0,{\mathbf{x}}_0|{\boldsymbol \phi}_{t})-{\mathbf{r}}_t^\top({\mathbf{x}}_0){\mathbf{R}}_t^{-1}{\mathbf{r}}_t({\mathbf{x}}_0) \\
&+ [{\mathbf{X}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)]^{\top}({\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{X}}_{t})^{-1}[{\mathbf{X}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)],
\end{split}\end{aligned}$$ where $\hat{{\mathbf{b}}}_t:=(\hat{{\boldsymbol \beta}}^\top, \hat{\gamma}_{t-1})^\top = ({\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t)^{-1}{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{y}}_t$, ${\mathbf{r}}_t({\mathbf{x}}_0):=r({{{\cal X}}}_t, {\mathbf{x}}_0| {\boldsymbol \phi}_t)$, ${\mathbf{X}}_1:={\mathbf{H}}_1$, and ${\mathbf{X}}_t:=[{{\mathbf{H}}}_{t}, y_{t-1}({{\cal X}}_t)]$ for $t>1$. $q_t$ is the number of columns in ${\mathbf{X}}_t$.

See Appendix \[app: conditional predictive dist\].

This result shows that one can generate a random sample from the predictive distribution $\pi({\mathbf{y}}({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \phi})$ by sequentially sampling from a collection of conditional distributions. Notice that samples are obtained across all the code levels. The computation associated with each conditional distribution only requires $O(n_t^3)$ flops for $t=1, \ldots, s$. Before stating the next theorem that provides a convenient way to exactly compute the predictive mean and predictive variance in the predictive distribution $\pi({\mathbf{y}}({\mathbf{x}}_0) \mid {\mathbf{y}}, {\boldsymbol \phi})$, we define the cokriging predictor and cokriging variance at each code level.

Let $\hat{{\mathbf{y}}}({\mathbf{x}}_0)=(\hat{y}_1({\mathbf{x}}_0), \ldots, \hat{y}_s({\mathbf{x}}_0))^\top$ be a vector of predictive means with $\hat{y}_t({\mathbf{x}}_0):=E[y_t({\mathbf{x}}_0)\mid {\mathbf{y}}, {\boldsymbol \phi}]$ and $\hat{{\mathbf{v}}}({\mathbf{x}}_0)=(\hat{v}_1({\mathbf{x}}_0),\ldots, \hat{v}_s({\mathbf{x}}_0))^\top$ be a vector of predictive variances with $\hat{v}_t({\mathbf{x}}_0):=Var[y_t({\mathbf{x}}_0)\mid {\mathbf{y}}, {\boldsymbol \phi}]$. In what follows, $\hat{{\mathbf{y}}}({\mathbf{x}}_0)$ is called the cokriging predictor and $\hat{{\mathbf{v}}}({\mathbf{x}}_0)$ is called the cokriging variance for all levels of code at new input ${\mathbf{x}}_0$.

\[thm: predictive mean and variance\] Suppose that $n_t-q_t>2$ such that the $t$-distributions in have valid variances. Then the cokriging predictor and cokriging variance at code level $t$ are given by $$\begin{aligned}
 \label{eqn: cokriging predictor}
    \begin{split}
    \hat{y}_t({\mathbf{x}}_0) &= {\mathbf{f}}_t^\top({\mathbf{x}}_0) \hat{{\mathbf{b}}}_t + {\mathbf{r}}_t^\top({\mathbf{x}}_0) {\mathbf{R}}_t^{-1}({\mathbf{y}}_t - {\mathbf{X}}_t \hat{{\mathbf{b}}}_t), \\
    \hat{v}_t({\mathbf{x}}_0) &= \hat{\gamma}_{t-1}^2\hat{v}_{t-1}({\mathbf{x}}_0) + \frac{n_t-q_t}{n_t-q_t-2} \hat{\sigma}^2_t \{r({\mathbf{x}}_0, {\mathbf{x}}_0|{\boldsymbol \phi}_t) - {\mathbf{r}}_t^\top({\mathbf{x}}_0) {\mathbf{R}}_t^{-1} {\mathbf{r}}_t({\mathbf{x}}_0) + \kappa_t \},
    \end{split}\end{aligned}$$ where ${\mathbf{f}}_1({\mathbf{x}}_0)={\mathbf{h}}_1({\mathbf{x}}_0)$, ${\mathbf{f}}_t({\mathbf{x}}_0)=[{\mathbf{h}}_t({\mathbf{x}}_0)^\top, \hat{y}_{t-1}^\top({\mathbf{x}}_0)]^\top$ for $t>1$, $\hat{v}_0:=0$ and $$\begin{aligned}
    \kappa_t &=  [{\mathbf{f}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)]^{\top}({\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{X}}_{t})^{-1}[{\mathbf{f}}_{t}({\mathbf{x}}_0)-{\mathbf{X}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}{\mathbf{r}}_t({\mathbf{x}}_0)] \\
    & + \hat{v}_{t-1}({\mathbf{x}}_0) \{ y_{t-1}^\top({{\cal X}}_t) {\mathbf{Q}}^H_t y_{t-1}({{\cal X}}_t) \}^{-1}.\end{aligned}$$ with ${\mathbf{Q}}^H_t = {\mathbf{R}}_t^{-1} - {\mathbf{R}}_t^{-1} {\mathbf{H}}_t ({\mathbf{H}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{H}}_t)^{-1} {\mathbf{H}}_t^\top {\mathbf{R}}_t^{-1}$.

See Appendix \[app: predictive mean and variance\].

Theorem \[thm: predictive mean and variance\] shows that 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})$.

In [@Gratiet2014], a recursive predictive formula is also derived based on a recursive multifidelity model that is different from what is presented in this article. The difference between the recursive multifidelity model in [@Gratiet2014] and the one presented here is that [@Gratiet2014] represent the high-fidelity code output $y_t(\cdot)$ as a function of the Gaussian process $y_{t-1}(\cdot)$ conditional on the code output ${\mathbf{y}}_{t-1}^{\mathscr{D}} :=\{ {\mathbf{y}}_1, \ldots, {\mathbf{y}}_{t-1} \}$ at design points in the experimental design sets $\{ \mathcal{X}_{\ell}: \ell=1, \ldots, t-1\}$ and model parameters $\{{\boldsymbol \beta}_{t-1}, \gamma_{t-1}, \sigma^2_{t-1}, {\boldsymbol \phi}_{t-1} \}$. Essentially, this recursive multifidelity model models the predictive distribution of $\pi(y_t(\cdot) \mid {\mathbf{y}}_{t-1}^{\mathscr{D}}, {\boldsymbol \beta}_{t-1}, \gamma_{t-1}, \sigma^2_{t-1}, {\boldsymbol \phi}_{t-1})$. This recursive multifidelity model only needs to estimate model parameters and prediction can be made automatically based on the model, however, it is worth noting that this predictive distribution does not take into account uncertainty due to estimation of model parameters $\{ {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi}\}$ and it does not give the predictive formula to predict code output at intermediate levels, $t=1, \ldots, s-1$, where such predictions may be useful to design the experiment with multifidelity codes.

The following corollary highlights the properties of autoregressive cokriging predictors in Theorem \[thm: predictive mean and variance\].

\[cor: cokriging and UK\] Let ${\mathbf{x}}_0$ be a new input in the domain $\mathcal{X}$. Let $\hat{y}_{t}^{K}({\mathbf{x}}_0):=E\{y_t({\mathbf{x}}_0) \mid {\mathbf{y}}_t, {\boldsymbol \phi}_t\}$ and $\hat{v}_t^{K}({\mathbf{x}}_0): = Var\{ y_t({\mathbf{x}}_0) \mid {\mathbf{y}}_t, {\boldsymbol \phi}_t\}$ be the kriging predictor and kriging variance based on data $\{{\mathbf{y}}_{t}, {{\cal X}}_{t}\}$. If ${\mathbf{x}}_0 \in {{\cal X}}_{t} \setminus {{\cal X}}_{t+1}$ and with $t=1, \ldots, s-1$, we have $\hat{y}_{\ell}({\mathbf{x}}_0) = \hat{y}_{\ell}^K({\mathbf{x}}_0)=y_t({\mathbf{x}}_0)$ and $\hat{v}_{\ell}({\mathbf{x}}_0) = \hat{v}_{\ell}^K({\mathbf{x}}_0) = 0$ for $\ell = 1, \ldots, t$. If ${\mathbf{x}}_0\notin {{\cal X}}_t$, we have $\hat{v}_t({\mathbf{x}}_0)>\hat{v}_t^K({\mathbf{x}}_0)$.

[Corollary]{} \[cor: cokriging and UK\] indicates that cokriging predictors can be interpolators as a kriging predictor. If ${\mathbf{x}}_0$ belongs to the design $\mathcal{X}_s$ in the highest fidelity code, the resulting predictive variances at ${\mathbf{x}}_0$ across all levels are zeros, i.e., $\hat{v}_t=0$ for $t=1, \ldots, s$. In other words, the cokriging predictors in  are interpolators at all levels. When prediction is made at new inputs, the cokriging predictor has larger variance than simply using the one associated with the kriging predictor. The extra uncertainty in cokriging variances comes from the uncertainty from previous level and the uncertainty to estimate the scale discrepancy parameter. However, without data coming from the lower code levels, one cannot estimate the correlation parameters in the higher level very well, since in practice the higher level code is too expensive to get sufficient number of runs that can be used to obtain fairly good parameter estimates.

Objective Bayesian Analysis {#sec: objective Bayes}
===========================

In Uncertainty Quantification, Bayesian analysis of cokriging models has been focused on using conjugate priors and noninformative priors [@Kennedy2000; @Qian2008; @Gratiet2013]. In this section, we focus on objective Bayesian analysis of the cokriging model, since objective priors can be used as default priors for Bayesian analysis [@Berger2006], and have been often used in Gaussian process modeling [@Berger2001; @Paulo2005; @Gu2018].

Commonly-used Improper Priors
-----------------------------

A Gaussian process has commonly-used priors; see [@Berger2001] for detailed discussions. Following this convention, the following priors will be referred to as commonly-used priors for parameters in autoregressive cokriging models. We consider the improper prior density for ${\boldsymbol \theta}:=\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi}\} \in \Omega=\mathbb{R}^{sp} \times \mathbb{R}^s \times (0, \infty)^s\times (0, \infty)^{sd}$ of the form $$\begin{aligned}
\label{eqn: general prior}
\begin{split}
    \pi({\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2},{\boldsymbol \phi})
    & \propto \frac{\pi({\boldsymbol \phi})}{\prod_{t=1}^s(\sigma^2_t)^{a_t}}, a_t \in \mathbb{R},
            \end{split}\end{aligned}$$ for various choices of $\pi({\boldsymbol \phi})$ and $a_t$. Although Kennedy and O’Hagan [@Kennedy2000] did not assume a prior for ${\boldsymbol \phi}$, their form is the same as the form in , where inverse range priors are chosen for ${\boldsymbol \phi}$: $\pi({\boldsymbol \phi}) = \prod_{t=1}^s \prod_{\ell=1}^d \phi_{t,\ell}^{-1}$ and $a_t=1$. The parameter ${\boldsymbol \phi}$ is estimated by maximizing the distribution $p({\mathbf{y}}\mid {\boldsymbol \gamma}, {\boldsymbol \sigma}, {\boldsymbol \phi})$. In Gratiet [@Gratiet2013], no prior is assumed for ${\boldsymbol \phi}$, and ${\boldsymbol \phi}$ is estimated based on a concentrated restricted likelihood via a restricted maximum likelihood (REML) approach [@Harville1974]. In Section \[sec: integrated likelihood\], we show that constant flat prior or inverse range prior for ${\boldsymbol \phi}$ yields improper posteriors.

Commonly-used Proper Priors
---------------------------

As the commonly-used improper priors may lead to improper posteriors, one obvious way to guarantee propriety of the posterior distribution is to assume proper priors, assessed either subjectively or from previous data, however, for Gaussian processes, the correlation parameter ${\boldsymbol \phi}$ can be difficult to interpret. Another way is to assume vague proper priors, however, this can only hide the problem when the posterior concentrates its mass at zero. If the posterior impropriety is occurring because the posterior is not decreasing at infinity, then the empirical Bayes estimates of ${\boldsymbol \phi}$ can be bad. In Appendix \[app: Qian\], we discuss that proper priors when chosen to be vague lead to improper posterior in [@Qian2008]. Another choice is to use conjugate priors for parameters $\{ {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2\}$ and leave the prior for ${\boldsymbol \phi}$ unspecified and work with a concentrated restricted likelihood; see [@Gratiet2013]. We will show that this concentrated restricted likelihood is not decreasing to zero and can be maximized either at zero or infinity; see Appendix \[app: Gratiet\] for detailed discussions. This problem has been studied in great details in [@Gu2018] for Gaussian processes. The conclusions in [@Gu2018] can be potentially generalized for autoregressive cokriging models. To overcome these problems, objective priors of the form  are developed and they are shown to yield proper posteriors. Although this article is not focused on robust estimation, the independent reference prior can lead to robust estimation with parameterization given in [@Gu2018].

Integrated Likelihood {#sec: integrated likelihood}
---------------------

The integration of the marginal likelihood  and the prior  with respect to the prior over $({\boldsymbol \beta}, {\boldsymbol \sigma}^2)$ yields $$\begin{aligned}
\begin{split}
\int L({{\mathbf{y}}}\mid{\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2},{\boldsymbol \phi})\pi({\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2},{\boldsymbol \phi})\,d({\boldsymbol \beta},{\boldsymbol \gamma},{\boldsymbol \sigma}^{2})
= L^I({\boldsymbol \phi}\mid {\mathbf{y}}) \pi({\boldsymbol \phi}),
\end{split}\end{aligned}$$ with $$\begin{aligned}
 \label{eqn: integrated likelihood}
\begin{split}
    L^I({\boldsymbol \phi}\mid {\mathbf{y}})     & \propto \prod_{t=1}^s |{{\mathbf{R}}}_t|^{-1/2} |{{\mathbf{X}}}_t^\top {{\mathbf{R}}}_t^{-1} {{\mathbf{X}}}_t|^{-1/2} \{S^2({\boldsymbol \phi}_t)\}^{-(({n}_t - q_t)/2+a_t-1)}, 
\end{split}\end{aligned}$$ where ${{\mathbf{X}}}_1={{\mathbf{H}}}_1$ and ${{\mathbf{X}}}_t=[{{\mathbf{H}}}_t, {W}_{t-1}]$ with ${W}_{t-1}:=y_{t-1}({{{\cal X}}}_t)$ for $t=2, \ldots, s$. $S^2({\boldsymbol \phi}_t):={\mathbf{y}}_{t}^\top {\mathbf{Q}}_t {\mathbf{y}}_{t}$ with ${\mathbf{Q}}_t={\mathbf{R}}_t^{-1}{\mathbf{P}}$ and ${\mathbf{P}}:=\mathbf{I} - {\mathbf{X}}_t ({\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t)^{-1} {\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1}$. It immediately follows that the posterior distribution of $\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2, {\boldsymbol \phi}\}$ is proper if and only if $$\begin{aligned}
 \label{eqn: condition for properiety}
    0< \int_{(0, \infty)^{sd}} L^I({\boldsymbol \phi}\mid {\mathbf{y}}) \pi({\boldsymbol \phi}) d {\boldsymbol \phi}< \infty. \end{aligned}$$ The following lemma gives the behavior of the integrated likelihood at zero and at infinity.

\[lem: behavior of integrated likelihood\] Note that ${\boldsymbol \phi}_t=(\phi_{t,1}, \ldots, \phi_{t,d})^\top$. Let $\mathcal{C}({\mathbf{X}}_t)$ be the column space of ${\mathbf{X}}_t$. For the cokriging model with sampling distribution  and prior distribution , under mild assumptions, we have

1.  If $\exists \ell$ such that $\phi_{t,\ell} \to 0$ for at least one $t$, the integrated likelihood exists and is greater than zero.

2.  If $\phi_{t,\ell} \to \infty$ for all $\ell$ and $t$, the integrated likelihood satisfies $$\begin{aligned}
            L^I({\boldsymbol \phi}\mid {\mathbf{y}}) =             O\left( \prod_{t=1}^s \left(\sum_{\ell=1}^d \nu_{t,\ell}(\phi_{t,\ell} )\right)^{a_t-c_t} \right),
        \end{aligned}$$

    where $c_t:= 1_{\{\mathbf{1}\in \mathcal{C}({\mathbf{X}}_t)\}} + \frac{1}{2} 1_{\{\mathbf{1}\notin \mathcal{C}({\mathbf{X}}_t)\}}$ with $1_{\{\mathbf{1}\in \mathcal{C}({\mathbf{X}}_t)\}}$ being 1 if the vector $\mathbf{1}$ is in $\mathcal{C}({\mathbf{X}}_t)$ and zero otherwise. The expression for $\nu_{t,\ell}(\phi_{t, \ell})$ is given in Appendix \[app: nu expression\], which is a continuous function and $\lim_{\phi_{t, \ell} \to \infty} \nu_{t,\ell}(\phi_{t, \ell}) = 0$.

See Appendix \[app: nu expression\].

Lemma \[lem: behavior of integrated likelihood\] shows that both the flat prior $\pi({\boldsymbol \phi}_t) \propto 1$ and the noninformative prior $\pi({\boldsymbol \phi}_t)\propto \prod_{\ell=1}^d \phi_{t,\ell}^{-1}$ with $a_t=1$ in lead to improper posteriors with the condition  violated.

Objective Priors
----------------

The posterior can be improper under certain common choices of priors in . In what follows, several objective priors are derived and they are shown to yield proper posteriors.

Following [@Berger2001], the parameters of interest are chosen to be $({\boldsymbol \sigma}^2, {\boldsymbol \phi})$ and $({\boldsymbol \beta}, {\boldsymbol \gamma})$ are treated as the nuisance parameters. This specification leads to the prior factorization $\pi^R({\boldsymbol \theta})=\pi^R({\boldsymbol \beta}, {\boldsymbol \gamma}\mid {\boldsymbol \sigma}^2, {\boldsymbol \phi})\pi^R({\boldsymbol \sigma}^2, {\boldsymbol \phi})$. The Jeffreys-rule prior $\pi^R({\boldsymbol \beta}, {\boldsymbol \gamma}\mid {\boldsymbol \sigma}^2, {\boldsymbol \phi})\propto 1$ is considered for the location parameters $({\boldsymbol \beta}, {\boldsymbol \gamma})$ when other parameters are assumed known. Then the reference prior $\pi^R({\boldsymbol \sigma}^2, {\boldsymbol \phi})$ is computed based on the integrated likelihood with respect to $\pi^R({\boldsymbol \beta}, {\boldsymbol \gamma}) \propto 1$. Standard calculation yields the integrated likelihood $L^I({\boldsymbol \sigma}^2, {\boldsymbol \phi}\mid {\mathbf{y}})$: $$\begin{aligned}
    \begin{split}
    L^I({\boldsymbol \sigma}^2, {\boldsymbol \phi}\mid {\mathbf{y}}) &= \int_{\mathbb{R}^{q+1}} L({\boldsymbol \theta}\mid {\mathbf{y}}) \pi^R({\boldsymbol \beta}, {\boldsymbol \gamma}) d ({\boldsymbol \beta}, {\boldsymbol \gamma}) \\
        & \propto \prod_{t=1}^s (\sigma^2_t)^{-(n_t-q_t)/2}|{\mathbf{R}}_t|^{-1/2}|{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t|^{-1/2} \exp\left \{-\frac{S^2({\boldsymbol \phi}_t)}{2\sigma^2_t} \right\}.
    \end{split}\end{aligned}$$

\[thm: reference prior\] Consider the group of parameters ${\boldsymbol \theta}=({\boldsymbol \theta}_1, \ldots, {\boldsymbol \theta}_s)$ with ${\boldsymbol \theta}_t=({\boldsymbol \beta}_t, \gamma_{t-1}, \sigma^2_t, {\boldsymbol \phi}_t)$, where $\gamma_0:=0$. For the cokriging model with sampling distribution , the independent reference prior distribution, $\pi^R({\boldsymbol \theta})$, is of the form  with $$\begin{aligned}
    a_t = 1 & \text{ and } \pi^R({\boldsymbol \phi}) \propto \prod_{t=1}^s|I_t^R({\boldsymbol \phi}_t)|^{1/2},\end{aligned}$$ where $I_t^R({\boldsymbol \phi}_t)$ is the Fisher information matrix by fixing all parameters except ${\boldsymbol \theta}_t$: $$\begin{aligned}
I_t^R({\boldsymbol \phi}_t) = 
\begin{pmatrix}
n_t-q_t & \text{tr}({\mathbf{W}}_{t,1})   & \text{tr}({\mathbf{W}}_{t,2})      & \cdots & \text{tr}({\mathbf{W}}_{t,d}) \\
    & \text{tr}({\mathbf{W}}_{t,1}^2) & \text{tr}({\mathbf{W}}_{t,1}{\mathbf{W}}_{t,2}) & \cdots & \text{tr}({\mathbf{W}}_{t,1}{\mathbf{W}}_{t,d}) \\ 
    &                     &          & \ddots    & \vdots  \\
    &                     &         &   & \text{tr}({\mathbf{W}}_{t,d}^2)
\end{pmatrix}_{(d+1)\times (d+1)},\end{aligned}$$ with ${\mathbf{W}}_{t,k}=\dot{{\mathbf{R}}}^k_t {\mathbf{Q}}_t, k=1, \ldots, d$,

See Appendix \[app: reference prior\].

\[thm: Jeffreys prior\] Let ${\boldsymbol \theta}=({\boldsymbol \theta}_1, \ldots, {\boldsymbol \theta}_s)$ be the group of parameters with ${\boldsymbol \theta}_t=({\boldsymbol \beta}_t, \gamma_{t-1}, \sigma^2_t, {\boldsymbol \phi}_t)$, where $\gamma_0:=0$. Suppose that ${\boldsymbol \theta}_t$’s are independent. Then the independent Jeffreys prior, $\pi^{J1}$, obtained by assuming that $({\boldsymbol \beta}_t, \gamma_t)$ and $(\sigma_t^2, {\boldsymbol \phi}_t)$ are a priori independent, and the independent Jeffreys prior, $\pi^{J2}$, are of the form  with $$\begin{aligned}
    a_t = 1 & \text{ and } \pi^{J1}({\boldsymbol \phi}) \propto \prod_{t=1}^s|I_t^J({\boldsymbol \phi}_t)|^{1/2}, \\
    a_t = 1 + q_t/2 & \text{ and } \pi^{J2}({\boldsymbol \phi}) \propto \pi^{J1}({\boldsymbol \phi}) \prod_{t=1}^s |{\mathbf{X}}_t^\top {\mathbf{R}}^{-1}_t {\mathbf{X}}_t|^{1/2} ,\end{aligned}$$ where $$\begin{aligned}
I_t^J({\boldsymbol \phi}_t) = 
\begin{pmatrix}
n_t & \text{tr}({\mathbf{U}}_{t,1})   & \text{tr}({\mathbf{U}}_{t,2})      & \cdots & \text{tr}({\mathbf{U}}_{t,d}) \\
    & \text{tr}({\mathbf{U}}_{t,1}^2) & \text{tr}({\mathbf{U}}_{t,1}{\mathbf{U}}_{t,2}) & \cdots & \text{tr}({\mathbf{U}}_{t,1}{\mathbf{U}}_{t,d}) \\ 
    &                     &          & \ddots    & \vdots  \\
    &                     &         &   & \text{tr}({\mathbf{U}}_{t,d}^2)
\end{pmatrix}_{(d+1)\times (d+1)},\end{aligned}$$ with ${\mathbf{U}}_{t,k}=\dot{{\mathbf{R}}}^k {\mathbf{R}}_t^{-1}, k=1, \ldots, d$,

This result follows directly from Proposition 2.2 in [@Paulo2005].

For the cokriging model with sampling distribution , the independent reference prior $\pi^R$, independent Jeffreys priors $\pi^{J1}, \pi^{J2}$ yield proper posteriors satisfying the condition .

Following the results in [@Paulo2005], it is straightforward to show that $$\begin{aligned}
    0<\int_{\mathbf{R}^d} L^I({\boldsymbol \phi}_t\mid{\mathbf{y}}) \pi({\boldsymbol \phi}_t)\, d{\boldsymbol \phi}_t < \infty,\quad & \text{ for } \quad  \pi({\boldsymbol \phi}_t)=\pi^R({\boldsymbol \phi}_t), \pi^{J1}({\boldsymbol \phi}_t), \pi^{J2}({\boldsymbol \phi}_t),\end{aligned}$$ where $L^I({\boldsymbol \phi}_t\mid{\mathbf{y}}) \propto |{\mathbf{R}}_t|^{-1/2}|{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t|^{-1/2} \{S^2({\boldsymbol \phi}_t)\}^{-(({n}_t - q_t)/2+a_t-1)}$, $\pi^R({\boldsymbol \phi}_t) \propto I_t^R({\boldsymbol \phi}_t)$, $\pi^{J1}({\boldsymbol \phi}_t) \propto I_t^J({\boldsymbol \phi}_t)$, $\pi^{J2}({\boldsymbol \phi}_t) \propto I_t^J({\boldsymbol \phi}_t) |{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t|^{1/2}$. Then the condition  is satisfied by Fubini’s theorem.

Parameter Estimation
--------------------

With the above prior specification, the integrated posterior of ${\boldsymbol \phi}$ given ${\mathbf{y}}$ is given by $$\begin{aligned}
 \label{eqn: posterior}
    \pi({\boldsymbol \phi}\mid {\mathbf{y}}) \propto  \prod_{t=1}^s |{{\mathbf{R}}}_t|^{-1/2} |{{\mathbf{X}}}_t^\top {{\mathbf{R}}}_t^{-1} {{\mathbf{X}}}_t|^{-1/2} \{S^2({\boldsymbol \phi}_t)\}^{-(({n}_t - q_t)/2 + a_t-1)} \pi({\boldsymbol \phi}_t),\end{aligned}$$ where $\pi({\boldsymbol \phi}_t)$ refers to independent reference prior and independent Jeffreys priors. Inference based on this posterior distribution can be made in a fully Bayesian paradigm via Markov chain Monte Carlo methods. Although uncertainties in all parameters can be taken into consideration in a fully Bayesian approach, the associated computation can be too expensive in practice due to repeated evaluation of the integrated likelihood .

In what follows, we focus on empirical Bayesian inference by maximizing the posterior  to obtain the estimate of ${\boldsymbol \phi}$. In fact, the numerical optimization can be performed for each ${\boldsymbol \phi}_t$, that is, for $t=1, \ldots, s$, $$\begin{aligned}
    \hat{{\boldsymbol \phi}}_t := \underset{{\boldsymbol \phi}_t}{\text{argmax}} \left\{-\frac{1}{2}\ln |{\mathbf{R}}_t| - \frac{1}{2}\ln|{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t| - \left(\frac{n_t-q_t}{2}+a_t-1\right)\ln S^2({\boldsymbol \phi}_t) + \ln \pi({\boldsymbol \phi}_t)\right\},\end{aligned}$$ where the maximization step can be performed using standard optimization algorithms such as the Nelder-Mead algorithm [@Nocedal2006]. Once $\hat{{\boldsymbol \phi}}_t$ is obtained, the cokriging predictions and cokriging variances can be obtained based on 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 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] except for the fact that objective priors and new predictive formulas are used.

In the following numerical examples, the covariance function model is specified as product form of the Matérn covariance following previous work [@Sacks1989]: $r(h) = \sigma^2 \prod_{i=1}^d r_i(h_i)$, where $$\begin{aligned}
r_i(h) = \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \sqrt{2\nu} \frac{h}{\phi_i} \right)^{\nu} \mathcal{K}_{\nu} \left( \sqrt{2\nu} \frac{h}{\phi_i} \right),\end{aligned}$$ where $\phi_i$ is the range parameter for the $i$ input dimension. $\nu>0$ is the smoothness parameter controlling the differentiability of the Gaussian processes. In the following numerical examples, $\nu$ will be fixed at 5/2 as in [@Gratiet2013] meaning that its random process will be twice differentiable in the mean square sense. $\mathcal{K}_{\nu}$ is the modified Bessel function of the second kind. In all the numerical examples, the range parameter will be reparameterized with the log inverse range parameters $\xi_i:=\log(1/\phi_i)$. This parametrization facilitates robust estimation for Gaussian process emulation as shown in [@Gu2018; @Gu2019] for reference priors and jointly robust priors. The predictive performance is measured based on root-mean-squared-prediction error (RMSPE), coverage probability of the 95% equal-tail credible interval (CVG(95%)), and average length of the 95% equal-tail credible interval (ALCI(95%)).

Testbed with the Borehole Function
----------------------------------

The performance is investigated with the 8-dimensional borehole function that models water flow through a borehole drilled from two ground surfaces through two aquifers [@Harper1983]. Its fast evaluation makes it widely used for testing purposes in computer experiments [@Morris1993; @Xiong2013]. Let ${\mathbf{x}}= (r_w, r, T_u, H_u, T_{\ell}, H_{\ell}, L, K_w)^\top$ be a vector of input variables in the borehole function with their physical meanings given in Appendix \[app: testing function\]. The response of the model is given by $$\begin{aligned}
    y_h = \frac{2\pi T_u (H_u - H_{\ell})}{\log(r/r_{w}) \left[ 1 + \frac{2LT_u}{\log(r/r_w)r_w^2 K_w} + T_u/T_{\ell} \right] },\end{aligned}$$ and its low-fidelity output is given by $$\begin{aligned}
    y_l = \frac{5 T_u (H_u - H_{\ell})}{\log(r/r_{w}) \left[ 1.5 + \frac{2LT_u}{\log(r/r_w)r_w^2 K_w} + T_u/T_{\ell} \right] }.\end{aligned}$$

To setup the experiment, 100 inputs are selected via Latin hypercube design with the `DiceDesign` package [@DiceDesign]. Then 20 inputs are randomly held out to evaluate predictive performance. The remaining 80 inputs are used to run the low-fidelity code $y_l(\cdot)$, and 30 inputs are randomly selected from these 80 inputs to run the high-fidelity code $y_h(\cdot)$. The predictive performance of the autoregressive cokriging model based on the proposed new formulas is compared with the approach in [@Gratiet2013], where the approach in [@Gratiet2013] is implemented in the `MuFiCokriging` package [@MuFiCokriging]. For all methods, the mean function is chosen to be constant and the covariance function is chosen to be the Matérn covariance with smoothness parameter fixed at 2.5 at each level.

Figure \[fig: borehole example\] compares predictive means against the high-fidelity code output at 20 held-out inputs. It indicates that the proposed approach gives better prediction than the approach in [@Gratiet2013], since the predictive values are concentrated along the 45 degree line. Table \[table: borehole example\] shows that the new cokriging predictors and cokriging variances give better predictive performance than the approach in [@Gratiet2013] in terms of RMSPE and ALCI(95%), although the empirical coverage probability is below 0.95. The result based on the jointly robust prior is obtained by fixing the hyperparameter $a_0$ at 0.2 after trying several different values. It is worth noting that the proposed cokriging variances take into account uncertainties in estimating all model parameters except the range parameters, but they still give much shorter predictive intervals. This also reveals that maximizing the posterior with the uniform improper priors for correlation parameters or maximizing the concentrated restricted likelihood are less preferred than maximizing the posterior with the proposed objective priors or the jointly robust prior. It is worth noting that the large predictive intervals are obtained based on the approach in [@Gratiet2013] even though uncertainties are not accounted for due to estimation of ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$. This is occurring because the concentrated restricted likelihood in [@Gratiet2013] can have nonrobust parameter estimates as discussed in Appendix \[app: nonrobust estimation\], and we found in this example that several correlation parameters are estimated to be very large, resulting in a nearly singular correlation matrix. The independent reference prior and independent Jeffreys prior lead to very similar predictive performance. Both of them have empirical coverage probability 0.85. In contrast, the jointly robust prior seems to give better predictive performance than these two objective priors.

\[table: borehole example\]

Application to Fluidized-Bed Processes
--------------------------------------

This section studies the predictive performance of the proposed cokriging formulas with objective priors for the fluidized-bed process experiment analyzed in [@Qian2008; @Gratiet2013]. The computer model named “Topsim" simulates the temperature of the steady-state thermodynamic operation point for a fluidized-bed process based on eight physical parameters: fluid velocity of the fluidization air, temperature of the air from the pump, flow rate of the coating solution, temperature of the coating solution, coating solution dry matter content, pressure of atomized air, room temperature, and humidity. [@Dewettinck1999] consider 28 different process conditions with coating solution used for distilled water (i.e., coating solution dry matter content is 0) and the room temperature at 20$^\circ$C. For each input configuration, one physical experiment $T_{exp}$ and three computer model runs ($T_1, T_2, T_3)$ were conducted, where $T_{exp}$ is the experimental response, $T_3$ is the most accurate code modeling the experiment, $T_2$ is a simplified version of $T_3$, and $T_1$ is the lowest accurate code modeling the experiment. The six inputs and corresponding outputs $T_1, T_2, T_3$, and $T_{exp}$ for these 28 runs are given in [@Qian2008]. The following numerical study mainly follows the procedure in [@Gratiet2013], and compares the proposed approach with that in [@Gratiet2013], since Gratiet [@Gratiet2013] demonstrates that the approach in [@Gratiet2013] performs better than the approach in [@Qian2008]. Following [@Gratiet2013], each input parameter is also scaled to the unit interval $[0, 1]$.

The predictive performance of a 2-level cokriging model is investigated based on 20 randomly selected $T_{exp}$ runs and all 28 $T_2$ runs. The remaining eight $T_{exp}$ runs are used for model validation. The mean function is chosen to be constant and the covariance function is chosen to be Matérn with smoothness parameter 2.5 at each level. In the jointly robust prior, the hyperparameter $a_0$ is chosen to be the default setting. The results in Table \[table: two-level example\] show that the 2-level cokriging model under the objective priors and jointly robust priors yields much smaller RMSPE and ALCI than the approach in [@Gratiet2013], while the 2-level cokriging models with all different priors have the empirical coverage probability smaller than 0.95.

Then the predictive performance of a 3-level cokriging model is investigated. Following [@Gratiet2013], 10 $T_{exp}$ runs (with the row number given by 1, 3, 8, 10, 12, 14, 18, 19, 20, 27 in Table 4 in [@Qian2008]), 20 $T_3$ runs (with the row number given by 1, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 19, 20, 22, 24, 27), and all the 28 $T_2$ runs are used. The remaining 18 $T_{exp}$ runs are used for mode validation. Table \[table: three-level example\] shows that the 3-level cokriging model gives much smaller RMSPE and ALCI with independent reference prior, independent Jeffreys prior, and jointly robust prior than the 3-level cokriging model with the approach in [@Gratiet2013]. In the jointly robust prior, the hyperparameter $a_0$ is chosen to be the default setting. Meanwhile, the independent reference prior and jointly robust prior result in much shorter predictive intervals than the approach in [@Gratiet2013]. The Jeffreys prior seems to perform very poor, since it only has empirical coverage probability at 0.5, which is much smaller than the nominal coverage probability 0.95. Notice that the jointly robust prior leads to better predictive performance than the independent reference prior with default settings for its hyperparameter, since it leads to better empirical coverage probability that the independent reference prior.

To briefly summarize, the objective priors and the independent jointly robust prior yield very similar predictive performance, and they yield much better predictive performance than the approach in [@Gratiet2013] in the 2-level cokriging model and the 3-level cokriging model. It is again worth noting that the approach in [@Gratiet2013] always gives larger predictive uncertainties than the proposed approach in these examples. This is due to the fact that the estimated parameters in [@Gratiet2013] are not robust, since several correlation parameters are estimated to be near zero. Such nonrobust estimates can mess up with the predictions. This suggests that inference with the objective priors should always be preferred over inference with the uniform improper priors or the approach in [@Gratiet2013] in autoregressive cokriging models when reliable experts’ opinions do not exist or are challenging to obtain and that the proposed cokriging predictors and cokriging variances are recommended for predictive inference.

\[table: two-level example\]

\[table: three-level example\]

Discussion {#sec: discussion}
==========

This article presents a unifying view in making prediction and parameter estimation in a computationally efficient way in the sense that the computational cost of both prediction and parameter estimation in an $s$-level cokriging model is the same as that in $s$ independent kriging models. The formulas in the predictive distributions account for uncertainties in all model parameters except the correlation parameters. The objective Bayesian analysis performed in the autoregressive cokriging model can be used as a default choice when prior information is challenging to obtain. In addition, the objective priors can also encourage robust estimation of correlation parameters.

The independent reference prior and independent Jeffreys priors are shown to yield proper posterior distributions. The predictive performance under these objective priors and independent jointly robust priors are also compared based on frequentist properties under various numerical studies. The numerical examples show that the objective priors and the jointly robust prior yield very similar predictive performance in 2-level and 3-level cokriging models. We also found that the jointly robust prior can provides better predictive performance than the independent reference prior sometimes, but this requires tuning its hyperparameters. The determination of optimal values for its hyperparameters is still an open question.

When designs are not hierarchically nested, there is no closed-form expression for the marginal likelihood function, and hence objective Bayesian analysis could be very challenging. However, the independent jointly robust prior could be a promising choice for this situation, since it has comparable predictive performance and it is a proper prior that allows fast computation. This has been used in [@Ma2019PPCokriging] for parameter estimation in a cokriging model that emulates high-dimensional output from multiple computer models.

Proofs {#app: Fisher Info}
======

Proof of Lemma \[lem: conditional predictive dists\] {#app: conditional predictive dist}
----------------------------------------------------

As $\{{{\cal X}}_t\cup\{{\mathbf{x}}_0\}: t=1, \ldots, s\}$ forms a collection of nested design, it follows form the cokriging model  that $$\begin{aligned}
    \pi({\mathbf{y}}({\mathbf{x}}_0), {\mathbf{y}}| {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \phi}) &= \pi(y_1({\mathbf{x}}_0), {\mathbf{y}}_1 | {\boldsymbol \beta}_1, \sigma^2_1, {\boldsymbol \phi}_1) \prod_{t=2}^s \pi(y_t({\mathbf{x}}_0), {\mathbf{y}}_t | y_{t-1}({\mathbf{x}}_0), {\boldsymbol \beta}_t, \gamma_{t-1}, \sigma^2_t, {\boldsymbol \phi}_t), \end{aligned}$$ where each joint distribution on the right hand side follows a multivariate normal distribution. It is straightforward to show that $$\begin{aligned}
    \pi({\mathbf{y}}({\mathbf{x}}_0) |{\mathbf{y}}, {\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \phi}) &= \pi(y_1({\mathbf{x}}_0) | {\mathbf{y}}_1, {\boldsymbol \beta}_1, \sigma^2_1, {\boldsymbol \phi}_1) \prod_{t=2}^s \pi(y_t({\mathbf{x}}_0) | {\mathbf{y}}_t, y_{t-1}({\mathbf{x}}_0), {\boldsymbol \beta}_t, \gamma_{t-1}, \sigma^2_t, {\boldsymbol \phi}_t), \end{aligned}$$ where each predictive distribution on the right hand side also follows a multivariate normal distribution. With location-scale priors for ${\boldsymbol \beta}_t, \gamma_{t-1}, \sigma^2_t$, standard calculations yield the formulas given in Lemma \[lem: conditional predictive dists\].

Proof of Theorem \[thm: predictive mean and variance\] {#app: predictive mean and variance}
------------------------------------------------------

The formula for the cokriging predictor at level $t$ follows from the law of total expectation as follows: $$\begin{aligned}
          E_{y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}}\{y_t({\mathbf{x}}_0) \} &=  E_{[y_{t-1}({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}]}\{E_{[y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}, y_{t-1}({\mathbf{x}}_0)]}[y_t({\mathbf{x}}_0)] \} \\
    &=E_{[y_{t-1}({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}]}\{{\mathbf{X}}_t^\top({\mathbf{x}}_0)\hat{{\mathbf{b}}}_t+r({\mathbf{x}}_0,{{{\cal X}}}_{t}|{\boldsymbol \phi}_{t}){\mathbf{R}}_t^{-1}({\mathbf{y}}_{t}-{\mathbf{X}}_t\hat{{\mathbf{b}}}_t) \} = \hat{y}_t({\mathbf{x}}_0),\end{aligned}$$ where standard calculations yield the last equality.

The formula for the cokriging variance at level $t$ follows from the law of total variance as follows: $$\begin{aligned}
   Var\{y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}\} &= Var\{ E[y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}, y_{t-1}({\mathbf{x}}_0)] \mid {\mathbf{y}}, {\boldsymbol \phi}\} \\
   &+ E\{ Var[y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}, y_{t-1}({\mathbf{x}}_0)] \mid {\mathbf{y}}, {\boldsymbol \phi}\},
        $$ with $$\begin{aligned}
    Var\{ E[y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}, y_{t-1}({\mathbf{x}}_0)] \mid {\mathbf{y}}, {\boldsymbol \phi}\} &= Var\{\hat{\gamma}_{t-1} y_{t-1}({\mathbf{x}}_0) | {\mathbf{y}}, {\boldsymbol \phi}\} = \hat{\gamma}_{t-1}^2 \hat{v}_{t-1}({\mathbf{x}}_0), \\
    E\{ Var[y_t({\mathbf{x}}_0)|{\mathbf{y}}, {\boldsymbol \phi}, y_{t-1}({\mathbf{x}}_0)] \mid {\mathbf{y}}, {\boldsymbol \phi}\} &= \frac{n_t-q_t}{n_t-q_t-2} E\{\hat{\sigma}_t^2 c_t^* | {\mathbf{y}}, {\boldsymbol \phi}\} ,\end{aligned}$$ where standard calculations yield that $E(c_t^*|{\mathbf{y}}, {\boldsymbol \phi}) = r({\mathbf{x}}_0, {\mathbf{x}}_0|{\boldsymbol \phi}_t) - {\mathbf{r}}_t^\top({\mathbf{x}}_0) {\mathbf{R}}_t^{-1} {\mathbf{r}}_t({\mathbf{x}}_0) + \kappa_t$.

Proof of Lemma \[lem: behavior of integrated likelihood\] {#app: nu expression}
---------------------------------------------------------

Given $1\leq t\leq s$ and $1\leq \ell \leq d$, suppose that $K(u) = r(u/\phi_{t,\ell})$ is a continuous function of $\phi_{t, \ell}>0$ for any $u>0$ such that:

-   $K(u) = r(u/\phi_{t,\ell})$, where $r(\cdot)$ is a correlation function satisfies $\lim_{u\to \infty} r(u) = 0$.

-   As $\phi_{t,\ell} \to \infty$, the correlation matrix satisfies $R_{t, \ell} = \mathbf{1} \mathbf{1}^\top + \nu_{t,\ell}(\phi_{t, \ell}) D_{t, \ell} + \nu_{t,\ell}(\phi_{t, \ell}) $ $\omega_{t,\ell}(\phi_{t, \ell}) (D^*_{t, \ell} + B_{t,\ell}(\phi_{t,\ell})) $, where $\nu_{t, \ell}(\phi_{t, \ell})>0$ is a continuous function of $\phi_{t,\ell}$, $D_{t,\ell}$ is a fixed nonsingular matrix with $\mathbf{1}^\top D_{t, \ell}^{-1} \mathbf{1}\neq 0$, $D^*_{t, \ell} $ is fixed matrix, and $B_{t,\ell}(\phi_{t,\ell})$ is a differential matrix satisfying $$\begin{aligned}
     \nu_{t,\ell}(\phi_{t, \ell}) \to 0, \quad \omega_{t, \ell}(\phi_{t,\ell}) \to 0, \quad \frac{\omega'_{t, \ell}(\phi_{t, \ell})}{ \frac{\partial \log \nu_{t, \ell}(\phi_{t, \ell}) }{\partial \phi_{t, \ell}} } \to 0, \\
     \| B_{t,\ell}(\phi_{t,\ell}) \|_{\infty} \to 0, \quad \frac{\| \frac{\partial B_{t,\ell}(\phi_{t,\ell})}{\partial \phi_{t,\ell}}  \|_{\infty}}{\frac{\partial}{\partial \phi_{t, \ell}} \log(\omega_{t, \ell}(\phi_{t, \ell}))} \to 0, 
     \end{aligned}$$ where $\| B\|_{\infty} = max_{i,j} |a_{i,j}|$ with $a_{i,j}$ being the $(i,j)$ entry of the matrix $B$.

These assumptions hold for all the correlation functions including power exponential, spherical, rational quadratic, and Matérn; see Table 1 in [@Gu2018]. Let the integrated likelihood at level $t$ be $$\begin{aligned}
\begin{split}
    L^I({\boldsymbol \phi}_t\mid {\mathbf{y}}) 
    & \propto  |{{\mathbf{R}}}_t|^{-1/2} |{{\mathbf{X}}}_t^\top {{\mathbf{R}}}_t^{-1} {{\mathbf{X}}}_t|^{-1/2} \{S^2({\boldsymbol \phi}_t)\}^{-(({n}_t - q_t)/2+a_t-1)}. 
\end{split}\end{aligned}$$ It follows from [@Gu2018] that

-   If $\exists \ell$ such that $\phi_{t,\ell} \to 0$ for at least one $t$, the integrated likelihood $L^I({\boldsymbol \phi}_t \mid {\mathbf{y}})$ at level $t$ exists and is greater than zero.

-   If $\phi_{t,\ell} \to \infty$ for all $\ell$ and $t$, the integrated likelihood at level $t$ satisfies $$\begin{aligned}
            L^I({\boldsymbol \phi}_{t} \mid {\mathbf{y}}) =  
            \begin{cases}
                O\left(  \left(\sum_{\ell=1}^d \nu_{t,\ell}(\phi_{t,\ell})\right)^{a_t-1/2} \right), & \mathbf{1} \notin  \mathcal{C}({\mathbf{X}}_t), \\
                O\left(  \left(\sum_{\ell=1}^d \nu_{t,\ell}(\phi_{t,\ell} )\right)^{a_t-1} \right), & \mathbf{1} \in \mathcal{C}({\mathbf{X}}_t).
            \end{cases}
        \end{aligned}$$

As $L^I({\boldsymbol \phi}\mid {\mathbf{y}}) = \prod_{t=1}^s L^I({\boldsymbol \phi}_t\mid {\mathbf{y}})$, the results in Lemma \[lem: behavior of integrated likelihood\] follow straightforwardly.

Proof of Theorem \[thm: reference prior\] {#app: reference prior}
-----------------------------------------

Arranging the parameters in the order ${\boldsymbol \vartheta}:=({\boldsymbol \vartheta}_1, \ldots, {\boldsymbol \vartheta}_s)$ with ${\boldsymbol \vartheta}_t:=(\sigma^2_t, {\boldsymbol \phi}_t^\top)^\top$, the Fisher information matrix $I^I({\boldsymbol \vartheta}_1, \ldots, {\boldsymbol \vartheta}_s)$ is computed from $\ell^I({\boldsymbol \vartheta}\mid {\mathbf{y}})$, whose $(i,j)$ entry is $$\begin{aligned}
    [I^I({\boldsymbol \vartheta}\mid {\mathbf{y}})]_{t,k} = E\left\{ \frac{\partial}{\partial \vartheta_{i}} \ell^I({\boldsymbol \vartheta}\mid {\mathbf{y}}) \times \frac{\partial}{\partial \vartheta_{j}} \ell^I({\boldsymbol \vartheta}\mid {\mathbf{y}})  \right\}.\end{aligned}$$ Differentiation with respect to $\sigma^2_t, \phi_{t,\ell}$ yields that $$\begin{aligned}
    \begin{split}
        \frac{\partial}{\partial \sigma^2_t} \ell^I({\boldsymbol \vartheta}\mid {\mathbf{y}}) = \frac{S^2_t-E(S^2_t)}{2\sigma^4},\quad & 
        \frac{\partial}{\partial\phi_{t,\ell}} \ell^I({\boldsymbol \vartheta}\mid {\mathbf{y}}) = 
        \frac{\Sigma^{\ell}_t - E(\Sigma^{\ell}_t)}{2\sigma^2_t}
    \end{split}\end{aligned}$$ where $S^2_t:={\mathbf{y}}_t^\top {\mathbf{Q}}_t{\mathbf{y}}_t$ with $S^2_t/\sigma^2_t \sim \chi_{n_t-q_t}^2$. $\Sigma^{\ell}_t$ is quadratic form on ${\mathbf{P}}_t{\mathbf{y}}_t\sim N(\mathbf{0}, \sigma^2_t{\mathbf{P}}_t{\mathbf{R}}_t)$ associated with the matrix ${\mathbf{R}}_t^{-1} \dot{{\mathbf{R}}}_t^{\ell} {\mathbf{R}}^{-1}$, where $\dot{{\mathbf{R}}}_t^{\ell}=\frac{\partial}{\partial \phi_{t,\ell}} {\mathbf{R}}_t$ is element-wise differentiation. Using results in [@Berger2001], the $(t,t)$ block diagonal matrix in the Fisher information matrix $I^I({\boldsymbol \vartheta})$ is $I^R_t({\boldsymbol \phi}_t)$.

Nonrobust Estimation {#app: nonrobust estimation}
====================

Posterior in Qian [@Qian2008] {#app: Qian}
-----------------------------

This section gives an example to show that posterior impropriety is occurring when vague priors in [@Qian2008] are chosen. As an illustrating example, we only discuss the posterior for correlation parameters at the first level. According to [@Qian2008], the following priors are assumed: $$\begin{aligned}
\pi({\boldsymbol \beta}\mid \sigma^2_1) &\sim \mathcal{N}({\mathbf{u}}_1, v_1\mathbf{I} \sigma^2_1) \\
\pi(\sigma^2_1) &\sim \mathcal{IG}(\alpha_1, \gamma_1^0) \\
\pi(\phi_{1,\ell}) &\sim \text{Gamma}(a_1^0, b_1^0), \ell=1, \ldots, d,\end{aligned}$$ where ${\mathbf{u}}_1, v_1, \alpha_1, \gamma_1^0, a_1^0, b_1^0$ are hyperparameters. With similar notations in [@Qian2008], the posterior distribution of ${\boldsymbol \phi}_1$ is, $$\pi({\boldsymbol \phi}_1 \mid {\mathbf{y}}_1) \propto \pi({\boldsymbol \phi}_1) |{\mathbf{R}}_1|^{-1/2} |{\mathbf{A}}_1|^{-1/2} |^{-1/2} \left\{ \gamma_1^0 + \frac{4c_1 - {\mathbf{B}}^\top_1{\mathbf{A}}_1^{-1} {\mathbf{B}}_1}{8}  \right \}^{-(\alpha_1+n/2)} ,$$ where ${\mathbf{A}}_1 = v_1^{-1} \mathbf{I} + {\mathbf{H}}_1^\top {\mathbf{R}}_1^{-1} {\mathbf{H}}_1$, ${\mathbf{B}}_1=-2v_1^{-1} {\mathbf{u}}_1 - 2{\mathbf{H}}_1^\top {\mathbf{R}}_1^{-1} {\mathbf{y}}_1$, $c_1= v_1^{-1} ({\mathbf{u}}^\top_1 {\mathbf{u}}_1) + {\mathbf{y}}_1^\top {\mathbf{R}}_1^{-1} {\mathbf{y}}_1$.

Thus, when $v_1 \to \infty, \alpha_1\to 0, \gamma_1^0\to 0$ such that priors for ${\boldsymbol \beta}_1$ and $\sigma^2_1$ become vague, the marginal posterior of ${\boldsymbol \phi}_1$ will be proportional to a product of the prior $\prod_{\ell=1}^d \mathcal{IG}({\boldsymbol \phi}_{1, \ell} \mid a_1^0, b_1^0)$ and the marginal likelihood $L({\boldsymbol \phi}_1\mid {\mathbf{y}})=|{\mathbf{R}}_1|^{-1/2} |{\mathbf{H}}_1^\top {\mathbf{R}}_1^{-1} {\mathbf{H}}_1|^{-1/2} (S^2)^{-n/2}$ with $S^2={\mathbf{y}}_1{\mathbf{Q}}_1 {\mathbf{y}}_1$, where ${\mathbf{Q}}_1={\mathbf{R}}_1^{-1} - {\mathbf{R}}_1^{-1}{\mathbf{H}}_1({\mathbf{H}}_1^\top {\mathbf{R}}_1^{-1}$ ${\mathbf{H}}_1)^{-1}{\mathbf{H}}_1^\top {\mathbf{R}}_1^{-1}$. This posterior will concentrate all its mass near 0 as $a_1^0 \to 0$ and $b_1^0\to 0$, resulting in nonrobust estimation according to Lemma 3.3 in [@Gu2018].

Marginal likelihood in Gratiet [@Gratiet2013] {#app: Gratiet}
---------------------------------------------

This section shows that the concentrated restricted likelihood can be maximized either at zero or infinity when noninformative priors or informative priors (when they are chosen to be vague) in [@Gratiet2013] are used. Gratiet [@Gratiet2013] considers two different types of priors for ${\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2$: noninformative priors and informative priors. The noninformative priors are chosen to be $$\begin{aligned}
\pi({\boldsymbol \beta}_1 \mid \sigma^2_1, {\boldsymbol \phi}_1) \propto 1, & & \pi(\sigma^2_1) \propto 1/\sigma^2_1, \\
\pi({\boldsymbol \beta}_t, {\boldsymbol \gamma}_{t-1} \mid \sigma^2_t, {\boldsymbol \phi}_t) \propto 1, & &\pi(\sigma^2_t) \propto 1/\sigma^2_t, t=2, \ldots, s. \end{aligned}$$ and the informative priors in [@Gratiet2013] are chosen to be $$\begin{aligned}
\pi({\boldsymbol \beta}_1 \mid \sigma^2_1, {\boldsymbol \phi}_1) \sim \mathcal{N}({\mathbf{b}}_1^0, \sigma^2_1 {\mathbf{V}}_1^0), & & \pi(\sigma^2_1 \mid {\boldsymbol \phi}_1) \sim \mathcal{IG}(\alpha_1^0, \gamma_1^0) \\
\pi(({\boldsymbol \beta}_t, \gamma_{t-1}) \mid \sigma_t^2, {\boldsymbol \phi}_t) \sim \mathcal{N}({\mathbf{b}}_t^0, \sigma^2_t {\mathbf{V}}_t^0), & & \pi(\sigma^2_t \mid {\boldsymbol \phi}_t) \sim \mathcal{IG}(\alpha_t^0, \gamma_t^0)\end{aligned}$$ Without further assuming a prior for $\pi({\boldsymbol \phi})$ as in [@Gratiet2013], Gratiet [@Gratiet2013] proposes to maximize the following concentrated restricted likelihood: $$\begin{aligned}
L_1({\boldsymbol \phi}_1 \mid {\mathbf{y}}, \hat{\sigma}^2_1) &\propto |{\mathbf{R}}_1|^{-1/2} (\hat{\sigma}^2_1)^{-(n_1-p_1)/2}, \\
L_t({\boldsymbol \phi}_t \mid {\mathbf{y}}, \hat{\sigma}^2_t) &\propto |{\mathbf{R}}_t|^{-1/2} (\hat{\sigma}^2_t)^{-(n_1-p_1-1)/2}, t=1, \ldots, s.\end{aligned}$$ For noninformative priors, $\hat{\sigma}^2_t \propto S^2({\boldsymbol \phi}_t)$. According to Lemma 3.3 in [@Gratiet2013], these marginal likelihood functions can have modes at ${\mathbf{R}}= \mathbf{I}_n$ and ${\mathbf{R}}= \mathbf{I}_n \mathbf{I}_n^\top$, resulting in nonrobust estimates for ${\boldsymbol \phi}$. For informative priors, the expression for $\hat{\sigma}^2_t$ is of the following form: $$\begin{aligned}
\hat{\sigma}^2_t \propto \gamma_t^0 + ({\mathbf{b}}_t - \bar{{\mathbf{b}}}_t)^\top \{ {\mathbf{V}}_t + ({\mathbf{X}}_t^\top {\mathbf{R}}^{-1} {\mathbf{X}}_t)^{-1} \}^{-1} ({\mathbf{b}}_t - \bar{{\mathbf{b}}}_t) + S^2({\boldsymbol \phi}_t),\end{aligned}$$ where ${\mathbf{b}}_t :=({\boldsymbol \beta}_t, \gamma_{t-1})^\top$ with $\gamma_0:=0$. $\bar{{\mathbf{b}}}_t$ is the generalized least square estimates for ${\mathbf{b}}_t$. When $\gamma_t^0 \to 0$ and ${\mathbf{V}}_t^{-1} \to \mathbf{0}$, $\hat{\sigma}^2_t \propto S^2({\boldsymbol \phi}_t)$. This reduces to the case when noninformative priors are used. Thus, estimates of the parameters ${\boldsymbol \phi}_t$ can be nonrobust. It is worth noting that [@Gratiet2013] chooses this proper prior to be informative instead of vague. It is crucial to perform sensitivity analysis whenever this prior is chosen to be informative.

Parameter Estimation for $\{{\boldsymbol \beta}, {\boldsymbol \gamma}, {\boldsymbol \sigma}^2 \}$ {#app: parameter estimation}
=================================================================================================

Let ${\mathbf{b}}_1={\boldsymbol \beta}_1$, ${\mathbf{b}}_t=({\boldsymbol \beta}^\top_t, \gamma_{t-1})^\top$ for $t>1$, and ${\mathbf{b}}=({\mathbf{b}}_1^\top, \ldots, {\mathbf{b}}_s)^\top$. The posterior distribution of ${\mathbf{b}}$ given ${\mathbf{y}}$ and $\hat{{\boldsymbol \phi}}$ with objective priors $\pi({\mathbf{b}}, {\boldsymbol \sigma}^2) \propto \prod_{t=1}^s \sigma^{-2}_{t}$ is $$\begin{aligned}
    \pi({\mathbf{b}}\mid {\mathbf{y}}, \hat{{\boldsymbol \phi}}) & \propto \int \sigma^{-2}_1 \pi({\mathbf{y}}_1\mid {\mathbf{b}}_1, \sigma^2_1, \hat{{\boldsymbol \phi}}_1) \prod_{t=2}^s \pi({\mathbf{y}}_t\mid {\mathbf{y}}_{t-1}, {\mathbf{b}}_t, \sigma^2_t, \hat{{\boldsymbol \phi}}_t) \sigma^{-2}_t \, d(\prod_{t=1}^s\sigma^2_t) \\
    & \propto \prod_{t}^s \{({\mathbf{y}}_t - {\mathbf{X}}_t{\mathbf{b}}_t)^\top {\mathbf{R}}_t^{-1}({\mathbf{y}}_t - {\mathbf{X}}_t{\mathbf{b}}_t)\}^{-n_t/2} |{\mathbf{R}}_t|^{-1/2}. \end{aligned}$$ Maximization with respect to this posterior distribution yields that $$\hat{{\mathbf{b}}}_t = ({\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t)^{-1} {\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{y}}_t.$$

Similarly, the posterior distribution of ${\boldsymbol \sigma}^2$ given ${\mathbf{y}}$ and $\hat{{\boldsymbol \phi}}$ with objective priors $\pi({\mathbf{b}}, {\boldsymbol \sigma}^2) \propto \prod_{t=1}^s \sigma^{-2}_{t}$ is $$\begin{aligned}
    \pi({\boldsymbol \sigma}^2 \mid {\mathbf{y}}, \hat{{\boldsymbol \phi}}) & \propto \int \sigma^{-2}_1 \pi({\mathbf{y}}_1\mid {\mathbf{b}}_1, \sigma^2_1, \hat{{\boldsymbol \phi}}_1) \prod_{t=2}^s \pi({\mathbf{y}}_t\mid {\mathbf{y}}_{t-1}, {\mathbf{b}}_t, \sigma^2_t, \hat{{\boldsymbol \phi}}_t) \sigma^{-2}_t \, d(\prod_{t=1}^s{\mathbf{b}}_t) \\
    & \propto \prod_{t=1}^s (\sigma_t^2)^{-(n_t-q_t)/2-1} |{\mathbf{R}}_t|^{-1/2} |{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t|^{-1/2} \exp\{-S^2(\hat{{\boldsymbol \phi}}_t)\}. \end{aligned}$$ It is easy to recognize that $\pi(\sigma_t^2\mid {\mathbf{y}}_{t-1}, {\mathbf{y}}_t, {\boldsymbol \phi}_t) = \mathcal{IG}((n_t-q_t)/2, S^2(\hat{{\boldsymbol \phi}}_t)/2)$. Hence, maximizing this posterior distribution with respect to $\sigma^2_t$ yields that $$\begin{aligned}
     \hat{\sigma}^2_t = S^2(\hat{{\boldsymbol \phi}}_t) / (n_t-q_t+2).
 \end{aligned}$$

Testing Function {#app: testing function}
================

\[table: input in borehole\]

Acknowledgement {#acknowledgement .unnumbered}
===============

This material is supported by the U.S. NSF under Grant DMS-1638521 to the Statistical and Applied Mathematical Sciences Institute. The author would like to thank Professor James O. Berger for his valuable comments and suggestions.

[^1]: Corresponds to: Pulong Ma, The Statistical and Applied Mathematical Sciences Institute and Duke University, 4051 Research Commons, Suite 300, 79 T.W. Alexander Drive, P.O. Box 110207, Durham, NC 27709, USA. Email: pulong.ma@duke.edu
