Impact of existence and nonexistence of pivot on the coverage of empirical best linear prediction intervals for small areas

The following two-level model, commonly referred to as the area-level model, has been extensively used in small area applications.

Level 1 (Sampling model): $y_{i}\,|\,\theta_{i}\stackrel{{\scriptstyle ind}}{{\sim}}\mathcal{N}(\theta_{i},D_{i})$;

Level 2 (Linking model): $\theta_{i}\stackrel{{\scriptstyle ind}}{{\sim}}\mathcal{G}(\bm{x_{i}^{\prime}\beta},A,\bm{\phi})$.

In the above model, level 1 is used to account for the sampling distribution of the direct estimates $y_{i}$, which are weighted averages of observations from small area $i$. Level 2 links the true small area means $\theta_{i}$ to a vector of $p$ known auxiliary variables $\bm{x_{i}}=(x_{i1},\cdots,x_{ip})^{\prime}$, often obtained from various administrative records. The parameters $\bm{\beta}$ and $A$ of the linking model are generally unknown and are estimated from the available data. As in other papers on the area-level model, the sampling variances $D_{i}$ are assumed to be known. In practice, $D_{i}$’s are estimated using a smoothing technique such as the ones given in Fay and Herriot (1979), Otto and Bell (1995), Ha et al. (2014) and Hawala and Lahiri (2018).

The two-level model can be rewritten as the following simple linear mixed model

where the random effect $u_{i}\overset{\mathrm{iid}}{\sim}\mathcal{G}(0,A,\bm{\phi})$ and the sampling error $e_{i}\overset{\mathrm{ind}}{\sim}\mathcal{N}(0,D_{i})$ are independent. When $\mathcal{G}$ is different from normal distribution, the best predictor (BP) of $\theta_{i}$, $\tilde{\theta_{i}}^{\rm{BP}}=\text{E}(\theta_{i}|y_{i})$ may not have a closed form. Instead, the best linear predictor (BLP) of $\theta_{i}$ always has the explicit form as below:

where $B_{i}=D_{i}/(A+D_{i})$, and the BLP minimizes the mean squared prediction error (MSPE) among all linear predictors of $\theta_{i}$. The variance of the prediction error $(\theta_{i}-\tilde{\theta_{i}}^{\rm{BLP}})$ is denoted by $g_{1i}(A)\coloneq\text{Var}(\theta_{i}-\tilde{\theta_{i}}^{\rm{BLP}})=AD_{i}/(A+D_{i})$. If $A$ is known, one can obtain the standard weighted least squares estimator of $\bm{\beta}$, denoting as $\tilde{\bm{\beta}}(A)$. Replacing $\bm{\beta}$ in (1.2) by $\tilde{\bm{\beta}}=\tilde{\bm{\beta}}(A)$, a best linear unbiased prediction (BLUP) estimator of $\theta_{i}$ is given by

which does not require normality assumptions typically used in small area estimation models. In practice, it is common that both $\bm{\beta}$ and $A$ are unknown and need to be estimated based on data. After plugging in their estimators, an empirical best linear unbiased predictor (EBLUP) of $\theta_{i}$ is given by:

where $\hat{B}_{i}=D_{i}/(\hat{A}+D_{i})$ and $\hat{\bm{\beta}}=\tilde{\bm{\beta}}(\hat{A})$, and $\hat{A}$ is a consistent estimator of $A$ for large $m$. After plugging in the estimator $\hat{A}$, we have $g_{1i}(\hat{A})=\hat{A}D_{i}/(\hat{A}+D_{i})$. To simplify the notation, throughout the remainder of this paper, we will use $g_{1i}$ and $\hat{g}_{1i}$ to denote $g_{1i}(A)$ and $g_{1i}(\hat{A})$, respectively. Point prediction using EBLUP and the associated mean squared prediction error (MSPE) estimation have been studied extensively. See Rao and Molina (2015) and Jiang and Nguyen (2007) for a detailed discussion on this subject.

In this paper, we consider prediction interval approximates of small area means $\theta_{i}$. A prediction interval of $\theta_{i}$, denoted by $I_{i}$, is called a $100(1-\alpha)\%$ interval for $\theta_{i}$ if $\mathrm{P}(\theta_{i}\in I_{i}|\bm{\beta},A,\bm{\phi})=1-\alpha$, for any fixed $\bm{\beta}\in\mathcal{R}^{p}$, $A\in\mathcal{R}^{+}$ and $\bm{\phi}\in\Theta_{\bm{\phi}}$, the parameter space of $\bm{\phi}$. The probability $\mathrm{P}$ is with respect to the area-level model. There are several options for constructing interval estimates for $\theta_{i}$. The prediction interval based on only the Level 1 model for the observed data is given by $y_{i}\pm z_{\alpha/2}\sqrt{D_{i}}$, where $z_{\alpha/2}$ is the $100(1-\alpha/2)$th standard normal percentile. While the coverage probability for this direct interval is $1-\alpha$, it is not efficient when $D_{i}$ is large as in the case of small area estimation. Hall and Maiti (2006) considered the interval based the regression synthetic estimator: $(\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{l}^{s}}\sqrt{\hat{A}},\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{u}^{s}}\sqrt{\hat{A}})$, where $b_{\hat{\alpha}_{l}^{s}}$ and $b_{\hat{\alpha}_{u}^{s}}$ are obtained using a parametric bootstrap method (described in detail in a later section). However, this interval is synthetic in the sense that it is constructed using synthetic regression estimator of $\theta_{i}$ and its associated uncertainty measure, which are not area specific in the outcome variable and hence may cause larger length of the confidence interval.

It is of importance to combine both levels of the model in the interval estimation. A effective approach is to use empirical best methodology. We call an interval empirical best (EB) prediction interval if it is based on empirical best predictor of $\theta_{i}$. For a special case of the area-level model that $\theta_{i}$ follows a normal distribution, Cox (1975) initiated the exact empirical Bayes interval: $\hat{\theta}_{i}^{\rm{EB}}\pm z_{\alpha/2}\sqrt{\hat{g}_{1i}}$, where $\hat{\theta}_{i}^{\rm{EB}}$ is the empirical Bayes estimator of $\theta_{i}$. Although the Cox interval always has smaller length than that of the direct interval, its coverage error is of the order $O(m^{-1})$, not accurate enough in most small area applications. Yoshimori and Lahiri (2014) improved the cox-type empirical Bayes interval by using a carefully devised adjusted residual maximum likelihood (ARML) estimator of $A$. Their interval has a coverage error of order $O(m^{-3/2})$. Additionally, they analytically showed that their interval always produces shorter length than the corresponding direct interval. However, the properties of both the ARML estimator of $A$ and the associated interval have not yet been explored for cases involving non-normally distributed random effects.

A function of the data and parameters is called to be a pivot if its distribution does not depend on any unknown quantity (Shao, 2008; Hall, 2013). When we have a linear mixed normal model on (1.1) , $(\theta_{i}-\tilde{\theta}_{i}^{\rm{BLP}})/\sqrt{g_{1i}}$ is a standard normal pivot. The traditional method of interval estimation for $\theta_{i}$ is of the form $\rm{EBLUP}\pm z_{\alpha/2}\sqrt{\text{mspe}}$, where mspe is an estimate of the true mean squared prediction error (MSPE) of the EBLUP. Unfortunately, $(\theta_{i}-\hat{\theta}_{i}^{\text{EBLUP}})/\sqrt{\hat{g}_{1i}}$ is not a pivot and this traditional approach produces too short or too long intervals. The coverage error of such interval is of the order $O(m^{-1})$, not accurate enough in most small area applications. Recognizing that $(\theta_{i}-\hat{\theta}_{i}^{\text{EBLUP}})/\sqrt{\hat{g}_{1i}}$ does not follow a standard normal distribution, Chatterjee et al. (2008) and Li and Lahiri (2010) developed a parametric bootstrap method to approximate its distribution and obtained a EB prediction interval for $\theta_{i}$ in linear mixed normal models. They showed such interval has the coverage error of the order $O(m^{-3/2})$. However, the property remains unknown for non-normal distributed random effects.

In this paper, one main aim is to bring out the virtues of pivoting or rescaling, which can decrease the dependence of our statistics on unknown parameters and yield improved prediction interval approximation for small areas under the general model (1.1). Analogous to Chatterjee et al. (2008), we propose parametric bootstrap methods to approximate the distribution of a suitably centered and scaled EBLUP under the general model (1.1), and apply it to construct a $100(1-\alpha)\%$ prediction interval estimation of $\theta_{i}$. Here, we define an interval based on the EBLUP of $\theta_{i}$ as an Empirical Best Linear (EBL) prediction interval. Specifically, we introduce these two key quantities: one is the centered and scaled EBLUP following the $F_{1i}$ distribution, which can be expressed as below

and the other is based on BLP which has the $F_{2i}$ distribution and can be expressed as

If $F_{2i}$ does not depend on any unknown parameters, $H_{i}(\bm{\beta},A)$ can be referred as a pivot; otherwise, $H_{i}(\bm{\beta},A)$ is not a pivot and we rewrite $F_{2i}$ as $F_{2i}(\bm{\nu})$, where $\bm{\nu}=(\nu_{1},\cdots,\nu_{s})^{\prime}$ is the unknown parameter vector determining the distribution of $H_{i}(\bm{\beta},A)$. We define $q_{\alpha}$ as the $\alpha$-level quantile of the distribution $F_{2i}(\bm{\nu})$, if $\mathrm{P}(H_{i}(\bm{\beta},A)\leq q_{\alpha}|\;\bm{\nu})=\alpha$ for any fixed $\bm{\nu}$.

In Section 2, we introduce a list of notation and regularity conditions used throughout the paper. In Section 3, we propose a single parametric bootstrap EBL prediction interval for small area means in the context of an area-level model, where the random effects follow a general (non-normal) known distribution, except possibly for unknown hyperparameters. We analytically demonstrate that the coverage error order of the EBL prediction interval remains the same as in Chatterjee et al. (2008), even when relaxing the normality of random effects through the existence of a pivot for suitably standardized random effects when hyperparameters are known. However, when a pivot does not exist, the EBL prediction interval fails to achieve the desired coverage error order, i.e. $O(m^{-3/2})$. Surprisingly, we find that the order $O(m^{-1})$ term is always positive under certain conditions, indicating potential overcoverage in the current single parametric bootstrap EBL prediction interval. Recognizing the challenge of showing existence of a pivot, we develop a simple moment-based method to claim non-existence of pivot. In Section 4, we propose a double parametric bootstrap method under a general area level model and, for the first time, analytically show that this approach can correct the coverage problem. In Section 5, we compare our proposed EBL prediction interval methods with the direct method, other traditional approaches, and the parametric bootstrap prediction interval proposed by Hall and Maiti (2006).

We use the following notations throughout the paper:

$\bm{Y}=(y_{1},\cdots,y_{m})^{\prime}$, a $m\times 1$ column vector of direct estimates;

$\bm{X^{\prime}}=(\bm{x_{1}},\cdots,\bm{x_{m}})$, a $p\times m$ known matrix of rank $p$;

$\bm{\Sigma}={\rm diag}(A+D_{1},\cdots,A+D_{m})$, a $m\times m$ diagonal matrix;

$\bm{\tilde{\beta}}=(\bm{X^{\prime}\bm{\Sigma}^{-1}X)^{-1}X^{\prime}\bm{\Sigma}^{-1}Y}$, weighted least square estimator of $\bm{\beta}$ with known $A$;

$\bm{P=\bm{\Sigma}^{-1}-\bm{\Sigma}^{-1}X(X^{\prime}\bm{\Sigma}^{-1}X)^{-1}X^{\prime}\bm{\Sigma}^{-1}}$;

$\partial F_{2i}(x;\tilde{\bm{\nu}})/\partial\bm{\nu}=\partial F_{2i}(x;\bm{\nu})/\partial\bm{\nu}|_{\bm{\nu}=\tilde{\bm{\nu}}}$, derivative with respect to $\bm{\nu}$ evaluated at $\tilde{\bm{\nu}}$.

We assume following regularity conditions for proving various results presented in this paper:

• R1:

  The rank of $\bm{X}$, $\text{rank}(\bm{X})=p$, is bounded for a large $m$;

• R2:

  $0<\inf_{i\geq 1}D_{i}\leq\sup_{i\geq 1}D_{i}<\infty$, $A\in(0,\infty)$ and the true $\bm{\phi}\in\Theta_{\bm{\phi}}^{0}$, the interior of $\Theta_{\bm{\phi}}$;

• R3:

  $\text{sup}_{i\geq 1}h_{ii}\equiv\bm{x_{i}^{\prime}(X^{\prime}X)^{-1}x_{i}}=O(m^{-1})$;

• R4:

  $\text{E}|u_{i}|^{8+\delta}<\infty$ for $0<\delta<1$;

• R5:

  The distribution function $F_{2i}(\cdot)$ is three times continuously differentiable with respect to ($\cdot$), and third derivative is uniformly bounded. When $H_{i}(\bm{\beta},A)$ is not a pivot, $F_{2i}(\cdot;\bm{\nu})$ is three times continuously differentiable with respect to the parameter vector ${\bm{\nu}}$, and its third derivative is uniformly bounded;

• R6:

  When having a non-pivot $H_{i}(\bm{\beta},A)$, we assume the estimator $\bm{\hat{\nu}}$ satisfies:
  

  $$\begin{split}\operatorname{E}[(\bm{\hat{\nu}}-\bm{\nu})^{\prime}\partial F_{2i}(x;\bm{\nu})/\partial\bm{\hat{\nu}}]=&\;O(m^{-1}),\\ \operatorname{E}\left[(\bm{\hat{\nu}}-\bm{\nu})^{\prime}\Bigl{\{}\partial^{2}F_{2i}(x;\bm{\nu})/\partial\bm{\hat{\nu}}\partial\bm{\hat{\nu}^{\prime}}\Bigr{\}}(\bm{\hat{\nu}}-\bm{\nu})\right]=&\;O(m^{-1}),\\ \text{E}\lVert(\bm{\hat{\nu}}-\bm{\nu})\rVert^{3}=\sum_{i,j,k}\text{E}[(\hat{\nu}_{i}-\nu_{i})(\hat{\nu}_{j}-\nu_{j})(\hat{\nu}_{k}-\nu_{k})]=&\;o(m^{-1}),\end{split}$$

  where the expectation is taken at the true $\bm{\nu}$. For a special case that $\bm{\nu}=A$, using the method of moment from Prasad and Rao (1990) to estimate $A$, Lahiri and Rao (1995) showed that $\text{E}(\hat{A}_{\rm{PR}}-A)=o(m^{-1})$. Moreover, they provided the Lemma C.1 that under the regularity conditions R1 - R5, $\text{E}|\hat{A}_{\rm{PR}}-A|^{2q}=O(m^{-q})$ for any $q$ satisfying $0<q\leq 2+\delta^{\prime}$ with $0<\delta^{\prime}<\frac{1}{4}\delta$.

For the remainder of this paper, without further explicit mention, we will use $\tilde{\theta}_{i}$ to represent the BLP of $\theta_{i}$, and $\hat{\theta}_{i}$ to denote the EBLUP of $\theta_{i}$. The single parametric bootstrap method has been widely studied for its simplicity in obtaining prediction intervals directly from the bootstrap histogram. Ideally, a prediction interval of $\theta_{i}$ can be constructed based on the distribution of $(\theta_{i}-\hat{\theta}_{i})/\sqrt{\hat{g}_{1i}}$, although this distribution function $F_{1i}$ is complex and difficult to approximate analytically. In this paper, we firstly follow the procedure introduced by Chatterjee et al. (2008) and Li and Lahiri (2010) to provide a bootstrap approximation of $F_{1i}$ by using a parametric bootstrap method. The implementation is straightforward, following these steps:

• 1.

  Conditionally on the data $(\bm{X,Y})$, draw $\theta_{i}^{*}$ for $i=1,\cdots,m$, independently from the distribution $\mathcal{G}(\bm{x_{i}}^{\prime}\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})$;

• 2.

  Given $\theta_{i}^{*}$, draw $y_{i}^{*}$ from the distribution $\mathcal{N}(\theta_{i}^{*},D_{i})$;

• 3.

  Construct the bootstrap estimators $\hat{\bm{\beta}}^{*}$, $\hat{A}^{*}$ using the data $(\bm{X,Y^{*}})$ and then obtain $\hat{\theta}_{i}^{*}=(1-\hat{B}^{*}_{i})y_{i}+\hat{B}^{*}_{i}\bm{x_{i}^{\prime}\hat{\beta}^{*}}$ and $\hat{g}_{1i}^{*}=\hat{A}^{*}D_{i}/(\hat{A}^{*}+D_{i})$;

• 4.

  Calibrate on $\alpha$ using the bootstrap method. Let $(q_{\alpha_{l}},q_{\alpha_{u}})=(q_{\hat{\alpha}_{l}^{s}},q_{\hat{\alpha}_{u}^{s}})$ be the solution of the following equations:

  $$\begin{split}\mathrm{P}^{*}(H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\leq q_{\alpha_{u}}\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})&=1-\alpha/2\\ \mathrm{P}^{*}(H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\leq q_{\alpha_{l}}\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})&=\alpha/2,\end{split}$$

  (3.1)

  where $H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})=(\theta_{i}^{*}-\hat{\theta}^{*}_{i})/\sqrt{\hat{g}^{*}_{1i}}$;

• 5.

  The single bootstrap calibrated prediction interval is constructed by:

  $$\hat{I}_{i}^{s}=\left(\hat{\theta}_{i}+q_{\hat{\alpha}_{l}^{s}}\sqrt{\hat{g}_{1i}}\,,\,\hat{\theta}_{i}+q_{\hat{\alpha}_{u}^{s}}\sqrt{\hat{g}_{1i}}\right).$$

One of our main results from the algorithm above is that we relax the normality of the random effects by the existence of a pivot $H_{i}(\bm{\beta},A)$, the prediction interval $\hat{I}_{i}^{s}$ obtained above still has a high degree of coverage accuracy. That is, it brings the coverage error down to $O(m^{-3/2})$.

The proof of Theorem 3.1 is given in the Appendix A. An example of Theorem 3.1 is that when $\mathcal{G}(\bm{x_{i}^{\prime}\beta},A,\bm{\phi})$ is a normal distribution, $\mathcal{N}(\bm{x_{i}^{\prime}\beta},A)$, using Theorem 3.2 of Chatterjee et al. (2008), we have

Hall and Maiti (2006) considered parametric bootstrap methods to approximate the distribution of $(\theta_{i}-\bm{x_{i}^{\prime}\hat{\beta}})/\sqrt{\hat{A}}$. Then, the prediction interval can be constructed as $(\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{l}^{s}}\sqrt{\hat{A}},\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{u}^{s}}\sqrt{\hat{A}})$, where $b_{\hat{\alpha}_{l}^{s}}$ and $b_{\hat{\alpha}_{u}^{s}}$ are obtained from single bootstrap approximation based on their algorithm. Their prediction interval is based on a synthetic model or the regression model, which does not permit approximation of the conditional distribution of $\theta_{i}$ given the data $\bm{Y}$. As a consequence, it is likely to underweight the area specific data. When the level 2 distribution $\mathcal{G}$ determined only by $\bm{x_{i}^{\prime}\beta}$ and $A$, it is easy to know that the quantity, $(\theta_{i}-\bm{x_{i}^{\prime}\beta})/\sqrt{A}$, is a pivot. As Hall and Maiti (2006) stated, their prediction interval is effective as $\mathrm{P}(\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{l}^{s}}\sqrt{\hat{A}}\leq\theta_{i}\leq\bm{x_{i}^{\prime}}\hat{\bm{\beta}}+b_{\hat{\alpha}_{u}^{s}}\sqrt{\hat{A}})=1-\alpha+O(m^{-2})$. However, when additional parameters involved into the distribution of random effects, $(\theta_{i}-\bm{x_{i}^{\prime}\beta})/\sqrt{A}$ might not a pivot and we may lose the effectiveness.

Hall and Maiti (2006) proposed a double-bootstrap method to calibrate $(\hat{\alpha}_{l}^{s},\hat{\alpha}_{u}^{s})$ from single parametric bootstraps, which can achieve a high degree of coverage accuracy $O(m^{-3})$. However, their calibration approach can overcorrect and produce a calibrated $\hat{\alpha}$ greater than 1, which makes the implementation infeasible.

While the order $O(m^{-1})$ term $T_{1}$ in (3.5) is theoretically positive under certain conditions, it remains unclear whether this positiveness holds when $u_{i}$ is asymmetrically distributed. Unlike Hall and Maiti (2006), we introduce a new double bootstrap method, which does not require a pivot or symmetric $u_{i}$. This method reduces the coverage error to $o(m^{-1})$, even when $u_{i}$ is asymmetrically distributed. Our double bootstrap approach is based on the algorithm from Shi (1992), where the double bootstrap is proposed to obtain accurate and efficient confidence interval for parameters of interest in both nonparametric and univariate parametric distribution settings. Later on, McCullough and Vinod (1998) discussed the theory of the double bootstrap both with and without pivoting, and provided the implementations in some nonlinear production functions. In this paper, we develop the double bootstrap in the context of our mixed effect model and apply it to obtain the EBL prediction intervals of $\theta_{i}$. The framework of our double parametric bootstrap is as below:

• 1.

  First-stage bootstrap

  • 1.1

    Conditionally on the data $(\bm{X,Y})$, draw $\theta_{i}^{*},\,i=1,\cdots,m$, from the distribution $\mathcal{G}(\bm{x_{i}}^{\prime}\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})$;

  • 1.2

    Given $\theta_{i}^{*}$, draw $y_{i}^{*}$ from the distribution $\mathcal{N}(\theta_{i}^{*},D_{i})$;

  • 1.3

    Compute $\hat{\bm{\beta}}^{*}$, $\hat{A}^{*}$, $\hat{\bm{\phi}}^{*}$ from the data $(\bm{X,Y^{*}})$, and obtain $\hat{\theta}_{i}^{*}$ and $\hat{g}_{1i}^{*}$;

• 2.

  Second-stage bootstrap

  • 2.1

    Given $\hat{\bm{\beta}}^{*}$, $\hat{A}^{*}$, $\hat{\bm{\phi}}^{*}$, draw $\theta_{i}^{**}$ from the distribution $\mathcal{G}(\bm{x_{i}}^{\prime}\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\bm{\phi}}^{*})$;

  • 2.2

    Given $\theta_{i}^{**}$, draw $y_{i}^{**}$ from the distribution $\mathcal{N}(\theta_{i}^{**},D_{i})$;

  • 2.3

    Compute $\hat{\bm{\beta}}^{**}$, $A^{**}$ from the data $(\bm{X,Y^{**}})$, also obtain $\hat{\theta}_{i}^{**}$ and $\hat{g}_{1i}^{**}$;

  • 2.4

    Consider to obtain a $(1-\alpha)$-level, two-sided, equal-tailed prediction interval. Define

    $$\hat{G}^{*}(z)\equiv\mathrm{P}^{**}(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq z\,|\,\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\phi}^{*})$$

    (4.2)

    For seeking the upper limit, we firstly solve the following system of equations in order to obtain $\hat{\alpha}_{u}$ such that

    $$\left\{\begin{array}[]{@{}l@{}}\mathrm{P}^{*}(H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\leq q_{\hat{\alpha}_{u}}\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})=1-\alpha/2\\ \mathrm{P}^{**}(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq q_{\hat{\alpha}_{u}}\,|\,\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\bm{\phi}}^{*})={\hat{\alpha}_{u}}.\end{array}\right.\,,$$

    (4.3)

    Using the definition (4.2), we have $q_{\hat{\alpha}_{u}}=\hat{G}^{*-1}(\hat{\alpha}_{u})$. Then, the above system of equations is equivalent to

    $$\mathrm{P}^{*}\left(H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\leq\hat{G}^{*-1}(\hat{\alpha}_{u})\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}}\right)=1-\alpha/2.$$

    (4.4)

    Rewriting (4.4) gives

    $$\mathrm{P}^{*}\left[\mathrm{P}^{**}\left(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\,|\,\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\bm{\phi}}^{*}\right)\leq\hat{\alpha}_{u}\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}}\right]=1-\alpha/2.$$

    (4.5)

    Note that the inner probability, $\mathrm{P}^{**}(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\,|\,\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\bm{\phi}}^{*})$, is a function of the first-stage bootstrap resample. More specifically, on the $j$th first-stage bootstrap sample, after all $K$ second-stage bootstrap operations are completed, let $Z_{j}$ be the proportion of times that $H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})$, i.e.,

    $$Z_{j}=\mathrm{P}^{**}\left(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\,|\,\hat{\bm{\beta}}^{*},\hat{A}^{*},\hat{\bm{\phi}}^{*}\right)\approx\#(H_{i}(\hat{\bm{\beta}}^{**},\hat{A}^{**})\leq H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*}))/K$$

    We will use $Z_{j}$ to adjust the first-stage intervals. After all bootstrapping operations are complete, we have estimates $Z_{j}$, $j=1,2,\cdots,J$, where $J$ is the number of first-stage bootstrap operations. Sort $Z_{j}$ and choose the $1-\alpha/2$ quantile as the percentile point $\hat{\alpha}_{u}$ for defining the double-bootstrap upper limit.

  • 2.5

    After completing the step 2.4, we have all estimates $H_{i}^{(j)}(\hat{\bm{\beta}}^{*},\hat{A}^{*})$, $j=1,\cdots,J$ and $\hat{\alpha}_{u}$. Choose $q_{\alpha_{u}}=q_{\hat{\alpha}^{d}_{u}}$ such that

    $$\mathrm{P}^{*}(H_{i}(\hat{\bm{\beta}}^{*},\hat{A}^{*})\leq q_{\alpha_{u}}\,|\,\hat{\bm{\beta}},\hat{A},\hat{\bm{\phi}})=\hat{\alpha}_{u}.$$

    (4.6)

  • 2.6

    Take $q_{\hat{\alpha}^{d}_{u}}$ as the upper limit of the double bootstrap prediction interval. Similar operations determine the lower limit, $q_{\hat{\alpha}^{d}_{l}}$. Finally, construct the prediction interval of $\theta_{i}$ as

    $$\hat{I}_{i}^{d}=\left(\hat{\theta}_{i}+q_{\hat{\alpha}^{d}_{l}}\sqrt{\hat{g}_{1i}}\,,\,\hat{\theta}_{i}+q_{\hat{\alpha}^{d}_{u}}\sqrt{\hat{g}_{1i}}\right).$$

The algorithm above shows that the single parametric bootstrap $\hat{I}^{s}_{i}$ is calibrated by the second-stage bootstrap. In general, such a calibration improves the coverage accuracy to $o(m^{-1})$, even when the pivot does not exist. One more advantage of our double bootstrap algorithm is that it avoids the problem of over correction and so make it more practical than that of Hall and Maiti (2006).

See proof of the theorem in appendix E.

In this section, we compare the performance of the proposed parametric bootstrap with their competitors where available, using Monte Carlo simulation studies. To maintain comparability with existing studies, we adopt part of the simulation framework of Chatterjee et al. (2008). We consider an area-level model with $\bm{x_{i}^{\prime}\beta}=0$, and five groups of small areas. Within each group, the sampling variances $D_{i}$’s remain the same. Two patterns for the $D_{i}$’s are considered: (i) $(4.0,0.6,0.5,0.4,0.2)$; (ii) $(8.0,1.2,1.0,0.8,0.4)$. For pattern (i), we take $A=1$; for pattern (ii), we take $A=2$ doubling the variances of pattern (i) while preserving the $B_{i}=D_{i}/(A+D_{i})$ ratios. To examine the effect of $m$, we consider $m=15$ and 50. With the increase of $m$, all methods improve and get closer to one another, supporting our asymptotic theory. Since we obtained virtually identical results under these two patterns, for the full study we confined attention to the pattern (i) and the results for pattern (ii) are provided in the Appendix G.

The Prasad-Rao method-of-moments (Prasad and Rao, 1990), and the Fay-Herriot method of estimating the variance component $A$ are considered. For the Fay-Herriot estimator of $A$, we employ the method of scoring to solve the estimating equation, which has showed to be more stable (Datta et al., 2005) than the Newton-Raphson method originally used in Fay and Herriot (1979). The estimation equation of the Fay-Herriot estimator is

Here, $f(A)$ is a decreasing function with respect to $A$, and the expectation of its first derivative $\text{E}(f^{\prime}(A))=-\rm{tr}(\bm{P})<0$. To improve computational efficiency, we implemented a slight modification to the algorithm. Specifically, the revised Fay-Herriot algorithm begins by calculating $f(A)$ at the initial point $A_{0}=0$, as in Fay and Herriot (1979). If $f(A_{0})<0$, we truncate the estimator to $\hat{A}_{\rm{FH}}=0.01$. Otherwise, the iterative process continues to search for a positive solution, with $\hat{A}_{\rm{FH}}=0.01$ also applied if no positive solution is found. This revised Fay-Herriot algorithm further enhances computational efficiency compared to the original method.