{#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 