Phase Transition Unbiased Estimation in High Dimensional Settings

TL;DR

This paper introduces a new property called phase transition unbiasedness for estimators in high-dimensional statistics, and develops an estimator that achieves this property using simulation techniques, improving bias control in complex models.

Contribution

The paper proposes a novel phase transition unbiasedness property and constructs an estimator that satisfies it, applicable in high-dimensional settings without requiring initial estimator consistency.

Findings

The estimator converges exponentially fast via iterative bootstrap.

It is unbiased for all sample sizes above a finite threshold.

The approach enhances robustness and avoids issues like separability in logistic regression.

Abstract

An important challenge in statistical analysis concerns the control of the finite sample bias of estimators. For example, the maximum likelihood estimator has a bias that can result in a significant inferential loss. This problem is typically magnified in high-dimensional settings where the number of variables $p$ is allowed to diverge with the sample size $n$. However, it is generally difficult to establish whether an estimator is unbiased and therefore its asymptotic order is a common approach used (in low-dimensional settings) to quantify the magnitude of the bias. As an alternative, we introduce a new and stronger property, possibly for high-dimensional settings, called phase transition unbiasedness. An estimator satisfying this property is unbiased for all $n$ greater than a finite sample size $n^\ast$. Moreover, we propose a phase transition unbiased estimator built upon the idea of matching an initial estimator computed on the sample and on simulated data. It is not required for this initial estimator to be consistent and thus it can be chosen for its computational efficiency and/or for other desirable properties such as robustness. This estimator can be computed using a suitable simulation based algorithm, namely the iterative bootstrap, which is shown to converge exponentially fast. In addition, we demonstrate the consistency and the limiting distribution of this estimator in high-dimensional settings. Finally, as an illustration, we use our approach to develop new estimators for the logistic regression model, with and without random effects, that also enjoy other properties such as robustness to data contamination and are also not affected by the problem of separability. In a simulation exercise, the theoretical results are confirmed in settings where the sample size is relatively small compared to the model dimension.

Tables

Table 1. Table 1: Simulation settings for the LRM.

Parameters	Setting I	Setting II
$p=$	$200$	$200$
$n=$	$2000$	$3000$
${\sum }_{i=1}^n{y}_i\approx$	$1000$	$750$
EPV $\approx$	5	3.75
$H=$	$500$	$500$
${\beta }_1={\beta }_2=$	$5$	$5$
${\beta }_3={\beta }_4=$	$-7$	$-7$
${\beta }_5=\dots ={\beta }_{200}=$	$0$	$0$
$\delta =$	$0.01$	$0.01$
Simulations	$500$	$500$

Table 2. Table 2: Simulation settings for the LRM with a random intercept.

Parameters	Setting I	Setting II
$p=q+2=$	$31$	$31$
$m=$	$5$	$50$
$\forall i,n_i=n_{\circ }=$	$50$	$5$
$n=$	$250$	$250$
${\sum }_{i=1}^m{\sum }_{j=1}^{n_i}{y}_{i,j}\approx$	$125$	$125$
EPV $\approx$	$4$	$4$
$H=$	$200$	$200$
${\beta }_0=$	$0$	$0$
${\beta }_1={\beta }_2=$	$5$	$5$
${\beta }_3={\beta }_4=$	$-7$	$-7$
${\beta }_5=\dots ={\beta }_{30}=$	$0$	$0$
${\sigma }^2=$	$1.5$	$1.5$
$\delta =$	$0.01$	$0.01$
Simulations	$1,000$	$1,000$

Table 3. Table 3: Dictionary between notations used in the main text and the appendices.

Main		Appendices
$𝐗({𝜽}_0,n)$	$\equiv$	$𝐗({𝜽}_0,n,{𝝎}_0)$
${𝐗}_i^{\ast }(𝜽,n)$	$\equiv$	$𝐗(𝜽,n,{𝝎}_i)$
$\widehat{𝝅}({𝜽}_0,n)$	$\equiv$	$\widehat{𝝅}({𝜽}_0,n,{𝝎}_0)$
${\widehat{𝝅}}_h^{\ast }(𝜽,n)$	$\equiv$	$\widehat{𝝅}({𝜽}_0,n,{𝝎}_h)$
${\widehat{𝜽}}^{\left(k\right)}$	$\equiv$	${\widehat{𝜽}}_{(n,H)}^{\left(k\right)}$
$\widehat{𝜽}$	$\equiv$	${\widehat{𝜽}}_{(n,H)}$
$\widehat{𝚯}$	$\equiv$	${\widehat{𝚯}}_{(n,H)}$
$𝐯(𝜽,n)$	$\equiv$	$𝐯(𝜽,n,𝝎)$

Equations

$${\hat{\mathbf{\theta }}}^{(k)}\equiv {\hat{\mathbf{\theta }}}^{(k-1)}+\left[\hat{\mathbf{\pi }}({\mathbf{\theta }}_0,n)-\frac{1}{H}\sum _{h=1}^H{\hat{\mathbf{\pi }}}_h^{\ast }\left({\hat{\mathbf{\theta }}}^{(k-1)},n\right)\right],$$

$$\hat{\mathbf{\Theta }}\equiv \underset{\mathbf{\theta }\in \mathbf{\Theta }}{\mathrm{argzero}}\hat{\mathbf{\pi }}({\mathbf{\theta }}_0,n)-\frac{1}{H}\sum _{h=1}^H{\hat{\mathbf{\pi }}}_h^{\ast }(\mathbf{\theta },n),$$

$$\hat{\bm{\theta}}^{\ast}\in\widehat{\bm{\Theta}}^{\ast}\equiv\operatorname*{argmin}_{\bm{\theta}\in\bm{\Theta}}\;\Big{\|}\hat{\bm{\pi}}(\bm{\theta}_{0},n)-\frac{1}{H}\sum_{h=1}^{H}\hat{\bm{\pi}}^{\ast}_{h}(\bm{\theta},n)\Big{\|}^{2}_{\bm{\Phi}}\;,$$

$$\underset{\mathbf{\theta }\in {\mathrm{I}\mathrm{R}}^p\setminus \mathbf{\Theta }}{\mathrm{argzero}}\hat{\mathbf{\pi }}({\mathbf{\theta }}_0,n)-\frac{1}{H}\sum _{h=1}^H{\hat{\mathbf{\pi }}}_h^{\ast }(\mathbf{\theta },n)=\varnothing .$$

$$\hat{\mathbf{\Theta }}=\underset{\mathbf{\theta }\in {\mathrm{I}\mathrm{R}}^p}{\mathrm{argzero}}\hat{\mathbf{\pi }}({\mathbf{\theta }}_0,n)-\frac{1}{H}\sum _{h=1}^H{\hat{\mathbf{\pi }}}_h^{\ast }(\mathbf{\theta },n).$$

$$\hat{\mathbf{\pi }}(\mathbf{\theta },n)=\mathbf{\pi }\left(\mathbf{\theta },n\right)+\mathbf{v}\left(\mathbf{\theta },n\right),$$

$${\mathbf{v}}_j(\mathbf{\theta },n)={\mathcal{O}}_{\mathrm{p}}(n^{-\alpha })\text{and}\underset{n\to \infty }{\lim }\frac{p^{\text{nicefrac}12}}{n^{\alpha }}=0.$$

$$\hat{\mathbf{\pi }}(\mathbf{\theta },n)=\mathbf{\theta }+\mathbf{d}\left(\mathbf{\theta },n\right)+\mathbf{v}\left(\mathbf{\theta },n\right).$$

$$\mathbf{d}\left(\mathbf{\theta },n\right)=\mathbf{a}(\mathbf{\theta })+\mathbf{c}(n)+\mathbf{b}\left(\mathbf{\theta },n\right),$$

$$\mathbf{b}\left(\mathbf{\theta },n\right)=\mathbf{L}(n)\mathbf{\theta }+\mathbf{r}\left(\mathbf{\theta },n\right),$$

$$\big{\lVert}\mathbf{a}(\bm{\theta}_{2})-\mathbf{a}(\bm{\theta}_{1})\big{\rVert}_{2}<\big{\rVert}\bm{\theta}_{2}-\bm{\theta}_{1}\big{\lVert}_{2}\,.$$

$${\mathbf{L}}_{j,l}(n)=\mathcal{O}(n^{-\beta }),{\mathbf{r}}_j\left(\mathbf{\theta },n\right)=\mathcal{O}(n^{-\gamma }),\underset{n\to \infty }{\lim }\frac{p^{\text{nicefrac}32}}{n^{\beta }}=0\text{and}\underset{n\to \infty }{\lim }\frac{p^{\text{nicefrac}12}}{n^{\gamma }}=0.$$

$$\underset{n\to \infty }{\lim }{p}^{\text{nicefrac}12}{c}_n=0.$$

$$\Big{\lVert}\mathbf{a}(\bm{\theta}_{2})-\mathbf{a}(\bm{\theta}_{1})\Big{\rVert}_{2}=\Big{\lVert}\mathbf{M}(\bm{\theta}_{2}-\bm{\theta}_{1})\Big{\rVert}_{2}\leq\big{\lVert}\mathbf{M}\big{\rVert}_{F}\big{\lVert}\bm{\theta}_{2}-\bm{\theta}_{1}\big{\rVert}_{2}<\big{\lVert}\bm{\theta}_{2}-\bm{\theta}_{1}\big{\rVert}_{2}.$$

$${\hat{\mathbf{\theta }}}^{(k)}\equiv {\hat{\mathbf{\theta }}}^{(k-1)}+{\epsilon }_k\left[\hat{\mathbf{\pi }}({\mathbf{\theta }}_0,n)-\frac{1}{H}\sum _{h=1}^H{\hat{\mathbf{\pi }}}_h^{\ast }({\hat{\mathbf{\theta }}}^{(k-1)},n)\right],$$

$$\mathbf{d}(\mathbf{\theta },n)=\sum _{j=1}^m\frac{{\mathbf{h}}^{(j)}(\mathbf{\theta })}{n^j}+\mathbf{g}(\mathbf{\theta },n),$$

$$\underset{n\to \infty }{\lim }\max \left(\frac{p^{\text{nicefrac}12}}{n^{\min (\alpha ,\gamma )}},\frac{p^{\text{nicefrac}32}}{n^{\beta }}\right)=0,$$

$$\underset{n\to \infty }{\lim }\frac{p^{\text{nicefrac}32}}{n}=0.$$

$$\underset{k\to \infty }{\lim }{\hat{\mathbf{\theta }}}^{(k)}=\hat{\mathbf{\theta }}.$$

$${\Vert {\hat{\mathbf{\theta }}}^{(k)}-\hat{\mathbf{\theta }}\Vert }_2={\mathcal{O}}_{\mathrm{p}}(p^{\text{nicefrac}12}{ϵ}^k).$$

$${\Vert {\hat{\mathbf{\theta }}}^{(k_n^{\ast })}-{\mathbf{\theta }}_0\Vert }_2\le {\Vert {\hat{\mathbf{\theta }}}^{(k_n^{\ast })}-\hat{\mathbf{\theta }}\Vert }_2+{\Vert \hat{\mathbf{\theta }}-{\mathbf{\theta }}_0\Vert }_2=o_{\mathrm{p}}(1),$$

$$\mathbf{a}(\mathbf{\theta })=\mathbf{M}\mathbf{\theta }+\mathbf{s},$$

$${\mathbf{L}}_{j,l}(n)=\mathcal{O}(n^{-\beta }),{\mathbf{r}}_j\left(\mathbf{\theta },n\right)=\mathcal{O}(n^{-\gamma }),\underset{n\to \infty }{\lim }\frac{p^{\text{nicefrac}32}}{n^{\beta }}=0\text{and}\underset{n\to \infty }{\lim }\frac{p^2}{n^{\gamma }}=0.$$

$$\underset{n\to \infty }{\lim }{p}^{\text{nicefrac}12}{c}_n=0.$$

$$\exists n^{\ast }>0,\forall \delta \in \mathcal{A},\exists M_{\delta }>0\text{such that}\forall n\ge n^{\ast },|f(n)|\le M_{\delta }g(n{)}^{\delta }.$$

$${\mathcal{O}}_{\delta \in \mathcal{A}}\left(g(n{)}^{\delta }\right)⟺\mathcal{O}\left(g(n{)}^{\delta }\right)\text{for all}\delta \in \mathcal{A},$$

$$f(n)=\mathcal{O}\left(g(n{)}^{\delta }\right)\text{for all}\delta \in \mathcal{A},$$

$$\forall \delta \in \mathcal{A},\exists n_{\delta }>0,\exists M_{\delta }>0\text{such that}\forall n\ge n_{\delta }|f(n)|\le M_{\delta }g(n{)}^{\delta }.$$

$${\mathcal{O}}_{\delta \in \mathcal{A}}\left(g(n{)}^{\delta }\right)⟺{\mathcal{O}}_{\delta \in \mathbb{N}}\left(g(n{)}^{\delta }\right).$$

$$\bar{M}_{\delta}\equiv\inf\left\{M>0\;\;\Bigg{|}\;\;\begin{array}[]{l}\exists n^{*}>0,\,\text{such that}\;\;|f(n)|\leq M\,g(n)^{\delta},\;\;\forall n\geq n^{\ast}\\ \text{and}\;\;\forall\delta^{\prime}\in\mathbb{N}\backslash\left\{\delta\right\},\,\exists M_{\delta^{\prime}}>0,\,|f(n)|\leq M_{\delta^{\prime}}\,g(n)^{\delta^{\prime}},\;\;\forall n\geq n^{\ast}\end{array}\right\}.$$

Full text

Keywords: Finite sample bias, Iterative bootstrap, Two-step estimators, Indirect inference, Robust estimation, Logistic regression.

1.   Introduction

An important challenge in statistical analysis concerns the control of the (finite sample) bias of estimators. For example, the Maximum Likelihood Estimator (MLE) has a bias that can result in a significant inferential loss. Indeed, under usual regularity conditions, the MLE is consistent and has a bias of asymptotic order $\mathcal{O}(n^{-1})$, where $n$ denotes the sample size. More generally, we qualify an estimator as asymptotically unbiased of order $\alpha$ if it has a bias of asymptotic order $\mathcal{O}(n^{-\alpha})$, elementwise, where $\alpha>0$. Thus, the MLE is typically asymptotically unbiased of order 1 and its bias vanishes as $n$ diverges. However, when $n$ is finite, the bias of the MLE may be large, and this bias is transferred to quantities such as test statistics where the MLE is used as a plug-in value (see e.g. Kosmidis, 2014a for a review). This problem is typically magnified in situations where the number of variables $p$ is large and possibly allowed to increase with $n$. For example, Sur and Candès (2019) show that the MLE for the logistic regression model can be severely biased in the case where $n$ and $p$ become increasingly large (given a fixed ratio).

The bias of an estimator can often be reduced by using bias correction methods that have recently received substantive attention. Bias correction is often achieved by simulation methods such as the jackknife (Efron, 1982) or the bootstrap (Efron, 1979), or can be approximated using asymptotic expansions (see e.g. Cordeiro and Vasconcellos, 1997, Cordeiro and Toyama Udo, 2008 and the references therein ). For example, the bootstrap bias correction (Efron and Tibshirani, 1994) allows to obtain an asymptotically unbiased estimator of order 2, provided that the initial estimator is consistent and asymptotically unbiased of order 1 (see Hall and Martin, 1988). An alternative approach is to correct the bias of an estimator by modifying its associated estimating equations. For example, Firth (1993) provides an adjustment of the MLE score function and this approach has been successively adapted and extended to Generalized Linear Models (GLM), among others, by Mehrabi and Matthews (1995); Bull et al. (2002); Kosmidis and Firth (2009, 2011); Kosmidis (2014b) . Several other bias correction methods have been proposed and a review on existing methods can, for example, be found in Kosmidis (2014a).

It is in general difficult to establish whether an estimator (possibly resulting from a bias correction technique) is unbiased for all $n$ (assuming that $n$ is large enough so that the estimator can be computed), and therefore the asymptotic order of the bias is used to “quantify” its magnitude. However, the asymptotic order of the bias has mainly been studied in low-dimensional settings and its extension in asymptotic regimes where both $n$ and $p$ can tend to infinity is often unclear. To address this, we introduce a new property to quantify the bias of an estimator in possibly high dimensional settings called Phase Transition (PT) unbiasedness. This property, presented in Definition 1 (below), can serve as a middle-ground between unbiasedness and the asymptotic order of the bias.

Definition 1 (Phase Transition unbiasedness):

An estimator $\hat{\bm{\theta}}$ of $\bm{\theta}_{0}$ is said to be PT-unbiased if there exists a $n^{\ast}\in\mathbb{N}^{*}\equiv\mathbb{N}\backslash\left\{0\right\}$ such that for all $n\in\mathbb{N}^{*}$ with $n\geq n^{*}$, we have \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0.

In short, this property implies that if an estimator is PT-unbiased, we have \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0 for all $n$ greater than a finite sample size $n^{\ast}$. Among others, this interpretation also explains the name of this property since, starting from a certain $n^{\ast}$, the estimator transitions from a biased phase to an unbiased phase. Therefore, PT-unbiasedness is weaker than the classical notion of unbiasedness but stronger than the asymptotic unbiasedness of order $\alpha$. Indeed, suppose that $\tilde{\bm{\theta}}$ is an asymptotically unbiased estimator of order $\alpha$, then we can write $\mathbb{E}[\tilde{\bm{\theta}}]=\bm{\theta}_{0}+\mathcal{O}\left(n^{-\alpha}\right)$, elementwise, for a specific $\alpha>0$. If $\hat{\bm{\theta}}$ is a PT-unbiased estimator, we have $\mathbb{E}[\hat{\bm{\theta}}]=\bm{\theta}_{0}+\mathcal{O}\left(n^{-\beta}\right)$, elementwise, for all $\beta>0$ and thus for all $\beta\geq\alpha$.

In this article, we propose a class of PT-unbiased estimators, called IB-estimators, building on the ideas of the Iterative Bootstrap (IB) (Kuk, 1995) and of indirect inference (Smith, 1993; Gourieroux et al., 1993). In short, these simulation-based methods allow to obtain consistent estimators from inconsistent initial ones. Moreover, Guerrier et al. (2019) study the order of the asymptotic bias of these estimators when the initial estimator is consistent. In particular, in low-dimensional settings, the latter work shows that these techniques present some strong similarities and can lead to asymptotically unbiased estimators of order $\delta$ where $\delta\in(2,3]$. In Section 3, we further study these simulation-based techniques in high dimensional settings and establish a connection between IB-estimators and PT-unbiasedness. In particular, we demonstrate, under suitable (and reasonable) conditions given in Section 2 that when the initial estimator is consistent or “slightly” asymptotically biased, the IB-estimator is a PT-unbiased estimator. In addition, we propose a two-step approach to attain PT-unbiasedness from inconsistent initial estimators. These new results have noticeable advantages in practical applications. As an illustration, in Section 4 we present a simulation study in a logistic regression setting (with and without random effects). We consider several slightly inconsistent initial estimators with different properties (e.g. “robustness” to separation and/or data contamination) and illustrate the effectiveness, in terms of finite sample bias and mean squared error, of the resulting IB-estimators. These simulation results are in line with the PT-unbiasedness property of IB-estimators. Moreover, they illustrate that IB-estimators can have additional properties (such as robustness) that are transferred from their initial estimator.

2.   Mathematical Setup

Let $\mathbf{X}(\bm{\theta}_{0},n)\in{\rm I\!R}^{n}$ denote a random sample generated under model $F_{\bm{\theta}_{0}}$ (possibly conditional on a set of fixed covariates), where $\bm{\theta}_{0}\in\bm{\Theta}\subset{\rm I\!R}^{p}$ is the parameter vector we wish to estimate. Moreover, we consider simulated samples that we denote as $\mathbf{X}_{h}^{*}(\bm{\theta},n)\in{\rm I\!R}^{n}$ for $\bm{\theta}\in\bm{\Theta}$ and where the subscript $h=1,\ldots,H$ identifies the distinct samples. Next, we define $\hat{\bm{\pi}}\left(\bm{\theta}_{0},n\right)$ and $\hat{\bm{\pi}}_{h}^{*}\left(\bm{\theta},n\right)$ as estimators of $\bm{\theta}_{0}$ and $\bm{\theta}$ based on the samples $\mathbf{X}(\bm{\theta}_{0},n)$ and $\mathbf{X}^{*}_{h}(\bm{\theta},n)$, respectively. In general, we assume that $\hat{\bm{\pi}}\left(\bm{\theta}_{0},n\right)$ is a biased estimator of $\bm{\theta}_{0}$ where this bias may be asymptotic and/or finite sample in nature.

In order to correct the bias of $\hat{\bm{\pi}}(\bm{\theta}_{0},n)$, similarly to Guerrier et al. (2019), we propose to use the IB to produce the sequence \Big{\{}\hat{\bm{\theta}}^{(k)}\Big{\}}_{k\in\mathbb{N}} defined as

[TABLE]

where $\hat{\bm{\theta}}^{(0)}\in\bm{\Theta}$ and $H$ is an integer chosen a priori. For example, when $\hat{{\bm{\pi}}}(\bm{\theta}_{0},n)\in\bm{\Theta}$, then this estimator can be used as the initial value of the above sequence. When this iterative procedure converges, we define the IB-estimator $\hat{\bm{\theta}}\in{\rm I\!R}^{p}$ as the limit in $k$ of $\hat{\bm{\theta}}^{(k)}$. It is possible that the sequence is stationary in that there exists a $k^{*}\in\mathbb{N}$ such that for all $k\geq k^{*}$ we have $\hat{\bm{\theta}}=\hat{\bm{\theta}}^{(k)}$. Moreover, when $\hat{\bm{\pi}}(\bm{\theta}_{0},n)$ in (1) is a consistent estimator of ${\bm{\theta}}_{0}$, the first step of the IB sequence (i.e. $k=1$) is equivalent to the standard bootstrap bias correction proposed by Efron and Tibshirani (1994) which can be used to reduce the bias of a consistent estimator (under some appropriate conditions).

Remark A:

Before proceeding to the presentation of the main results of this article, the authors would like to make the reader aware of the fact that the notation used in the article is slightly different from the one used in its appendices where most of the proofs are presented. Indeed, a simplified notation is used in the article to enhance readability while a more detailed and precise notation is used in the appendices to prove the results and avoid misleading conclusions. The notation used in the appendices is presented and explained in Appendix A.

In order to present the approach, it is useful to define the set $\widehat{\bm{\Theta}}$ as

[TABLE]

which will allow us to study the properties of the sequence $\hat{\bm{\theta}}^{(k)}$. As one may expect, under some suitable conditions, we have that $\hat{\bm{\theta}}\in\widehat{\bm{\Theta}}$ thereby explaining why some of our assumptions concern $\bm{\Theta}$.

Remark B:

The elements of $\widehat{\bm{\Theta}}$ are in fact special cases of an indirect inference estimator. Indeed, using the previous notation, an indirect inference estimator can be defined as

[TABLE]

where $\bm{\Phi}$ is a positive-definite matrix. In the case where $\widehat{\bm{\Theta}}\neq\emptyset$, we have that $\widehat{\bm{\Theta}}=\widehat{\bm{\Theta}}^{\ast}$, which clearly shows the equivalence between (2) and (3), for any positive-definite matrix $\bm{\Phi}$. Note that $\hat{\bm{\pi}}(\bm{\theta},n)$ is referred to as the auxiliary estimator in the indirect inference framework.

Before presenting the properties of the IB sequence, we first describe the assumptions we consider.

Assumption A:

Let $\bm{\Theta}$ be a convex and compact subset of ${\rm I\!R}^{p}$ such that ${\bm{\theta}}_{0}\in\operatorname{Int}(\bm{\Theta})$ and

[TABLE]

Assumption A is quite mild although stronger than necessary. Indeed, the convexity condition and the fact that $\bm{\theta}_{0}$ is required to be in the interior of $\bm{\Theta}$ are convenient to ensure that expansions can be made between $\bm{\theta}_{0}$ and an arbitrary point in $\bm{\Theta}$. Similarly, the compactness is particularly convenient as it allows to bound certain quantities. Finally, the last part of the assumption essentially ensures that no solution outside of $\bm{\Theta}$ exists for the estimator defined in (2). Therefore, the solution set $\widehat{\bm{\Theta}}$ may equivalently be written as

[TABLE]

Next, we impose some conditions on the initial estimator $\hat{\bm{\pi}}(\bm{\theta},n)$. Denoting $\mathbf{a}_{j}$ as being the $j$-th entry of a generic vector $\mathbf{a}\in{\rm I\!R}^{p}$, the assumption is as follows.

Assumption B:

For all $\left(\bm{\theta},\,n\right)\in\bm{\Theta}\times\mathbb{N}^{\ast}$, the expectation $\bm{\pi}\left(\bm{\theta},n\right)\equiv\mathbb{E}\left[\hat{\bm{\pi}}(\bm{\theta},n)\right]$ exists and is finite, i.e. $\lvert\bm{\pi}_{j}\left(\bm{\theta},n\right)\rvert<\infty$ for all $j=1,\ldots,p$. Moreover, for all $j=1,\ldots,p$ $\bm{\pi}_{j}\left(\bm{\theta}\right)\equiv\displaystyle{\lim_{n\to\infty}}\bm{\pi}_{j}\left(\bm{\theta},n\right)$ exists.

Assumption B is likely to be satisfied in the majority of practical situations and is particularly useful as it allows to decompose the estimator $\hat{\bm{\pi}}(\bm{\theta},n)$ into a non-stochastic component $\bm{\pi}\left(\bm{\theta},n\right)$ and a random term $\mathbf{v}\left(\bm{\theta},n\right)$. Indeed, using Assumption B, we can write:

[TABLE]

where $\mathbf{v}\left(\bm{\theta},n\right)\equiv\hat{\bm{\pi}}(\bm{\theta},n)-\bm{\pi}\left(\bm{\theta},n\right)$ is a zero-mean random vector. It must be noticed that this assumption does not imply that the functions $\bm{\pi}\left(\bm{\theta},n\right)$ and $\hat{\bm{\pi}}(\bm{\theta},n)$ are continuous in ${\bm{\theta}}$. Consequently, when using the IB (see Theorem 1 below), it is possible to make use of initial estimators that are not continuous in ${\bm{\theta}}$. In Guerrier et al. (2019), the continuity of these functions was (implicitly) assumed while the weaker assumptions presented in this article extend the applicability of this framework also to models for discrete data such as logistic regression, which is discussed in Section 4. We now move onto Assumption C (below) that imposes some restrictions on the random vector $\mathbf{v}\left(\bm{\theta},n\right)$.

Assumption C:

The second moment of $\mathbf{v}(\bm{\theta},n)$ exists and there is a real $\alpha>0$ such that for every $\bm{\theta}\in\bm{\Theta}$ and for all $j=1,\,\ldots,\,p$, we have

[TABLE]

Assumption C is frequently employed and typically very mild as it simply requires that the variance of $\mathbf{v}_{j}\left(\bm{\theta},n\right)$ exists and goes to zero as $n$ increases. For example, if $\hat{\bm{\pi}}(\bm{\theta},n)$ is $\sqrt{n}$-consistent (towards $\bm{\pi}\left(\bm{\theta}\right)$) then we would have $\alpha=\nicefrac{{1}}{{2}}$ and Assumption C would simply require that $p=o(n)$. Next, we consider the bias of $\hat{\bm{\pi}}(\bm{\theta},n)$ for $\bm{\theta}$ and we let $\mathbf{d}\left(\bm{\theta},n\right)\equiv\bm{\pi}(\bm{\theta},n)-\bm{\theta}$. Using this definition we can rewrite (4) as follows

[TABLE]

Moreover, the bias function $\mathbf{d}\left(\bm{\theta},n\right)$ can always be expressed as follows

[TABLE]

where $\mathbf{a}({\bm{\theta}})$ is defined as the asymptotic bias in the sense that $\mathbf{a}_{j}(\bm{\theta})\equiv\displaystyle{\lim_{n\to\infty}}\;\mathbf{d}_{j}\left(\bm{\theta},n\right)$ for $j=1,\ldots,p$, while $\mathbf{c}(n)$ and $\mathbf{b}\left(\bm{\theta},n\right)$ are used to represent the finite sample bias. More precisely, $\mathbf{a}({\bm{\theta}})$ contains all the terms that are strictly functions of ${\bm{\theta}}$, $\mathbf{c}(n)$ the terms that are strictly functions of $n$ and $\mathbf{b}\left(\bm{\theta},n\right)$ the rest. This definition implies that if $\mathbf{d}\left(\bm{\theta},n\right)$ contains a constant term, it is included in $\mathbf{a}({\bm{\theta}})$. Moreover, the function $\mathbf{b}\left(\bm{\theta},n\right)$ can always be decomposed into a linear and a non-linear term in ${\bm{\theta}}$, i.e.

[TABLE]

where $\mathbf{L}(n)\in{\rm I\!R}^{p\times p}$ and $\mathbf{r}\left(\bm{\theta},n\right)$ does not contain any linear term in $\bm{\theta}$. Denoting $\mathbf{A}_{j,l}$ as the entry in the $j$-th row and $l$-th column of a generic matrix $\mathbf{A}\in{\rm I\!R}^{p\times p}$, we consider Assumption D which imposes some restrictions on the different terms of the bias function.

Assumption D:

The bias function $\mathbf{d}\left(\bm{\theta},n\right)$ is as follows:

The function $\mathbf{a}(\bm{\theta})$ is a contraction map in that for any $\bm{\theta}_{1},\bm{\theta}_{2}\in\bm{\Theta}$ such that $\bm{\theta}_{1}\neq\bm{\theta}_{2}$ we have

[TABLE] 2. 2.

There exist real $\beta,\gamma>0$ such that for all $\bm{\theta}\in\bm{\Theta}$ and any $j,l=1,\dots,p$, we have

[TABLE] 3. 3.

Defining $c_{n}\equiv\displaystyle{\max_{j=1,\dots,p}}\mathbf{c}_{j}(n)$ for all $n\in\mathbb{N}^{*}$, we require that the sequence $\left\{c_{n}\right\}_{n\in\mathbb{N}^{*}}$ is such that

[TABLE]

The first part of Assumption D is reasonable provided that the asymptotic bias is relatively “small” compared to ${\bm{\theta}}$ (up to a constant term). For example, if $\mathbf{a}({\bm{\theta}})$ is sublinear in ${\bm{\theta}}$, i.e. $\mathbf{a}({\bm{\theta}})=\mathbf{M}{\bm{\theta}}+\mathbf{s}$, then the first part of Assumption D would be satisfied if the Frobenius norm is such that $||\mathbf{M}||_{F}<1$ since

[TABLE]

Moreover, the function ${\bm{\pi}}(\bm{\theta})$ (as defined in Assumption B) is often called the asymptotic binding function in the indirect inference literature (see e.g. Gourieroux et al., 1993). To ensure the consistency of such estimators, it is typically required for ${\bm{\pi}}(\bm{\theta})$ to be continuous and injective. In our setting, this function is given by ${\bm{\pi}}(\bm{\theta})={\bm{\theta}}+\mathbf{a}({\bm{\theta}})$ and its continuity and injectivity are directly implied by the first part of Assumption D. Indeed, since $\mathbf{a}(\bm{\theta})$ is a contraction map, it is continuous. Moreover, taking $\bm{\theta}_{1},\bm{\theta}_{2}\in\bm{\Theta}$ with ${\bm{\pi}}(\bm{\theta}_{1})={\bm{\pi}}(\bm{\theta}_{2})$, then \big{\lVert}{\mathbf{a}}(\bm{\theta}_{1})-{\mathbf{a}}(\bm{\theta}_{2})\big{\rVert}=\big{\lVert}\bm{\theta}_{1}-\bm{\theta}_{2}\big{\rVert}, which is only possible if ${\bm{\theta}}_{1}={\bm{\theta}}_{2}$. Thus, $\bm{\pi}(\bm{\theta})$ is injective.

Remark C:

In situations where $\mathbf{a}(\bm{\theta})$ is not a contraction map, a possible solution is to modify the sequence considered in the IB as follows:

[TABLE]

with $\varepsilon_{k}\in(0,1]$ for all $k\in\mathbb{N}$. If $\varepsilon_{k}=\varepsilon$ (i.e. a constant), $\mathbf{a}(\bm{\theta})$ does not need to be a contraction map. Indeed, if $\bm{\Theta}$ is bounded and $\mathbf{a}(\bm{\theta})$ is differentiable, it is always possible to find an $\varepsilon$ such that $\varepsilon\,\mathbf{a}(\bm{\theta})$ is a contraction map. A formal study on the influence of $\varepsilon_{k}$ on the IB algorithm is, however, left for further research.

While more general, the second part of Assumption D would be satisfied, for example, if $\mathbf{b}(\bm{\theta},n)$ is a sufficiently smooth function in $\bm{\theta}$ and/or $n$, thereby allowing a Taylor expansion, as considered, for example, in Guerrier et al. (2019). Moreover, this assumption is typically less restrictive than the approximations that are commonly used to describe the bias of estimators. Indeed, a common assumption is that the bias of a consistent estimator (including the MLE), can be expanded in a power series in $n^{-1}$ (see e.g. Kosmidis, 2014a and Hall and Martin, 1988 in the context of the iterated bootstrap), i.e.

[TABLE]

where $\mathbf{h}^{(j)}(\bm{\theta})$ is $\mathcal{O}(1)$ elementwise, for $j=1,\ldots,m$, and $\mathbf{g}({\bm{\theta}},n)$ is $\mathcal{O}\left(n^{-(m+1)}\right)$ elementwise, for some $m\geq 1$. The bias function $\mathbf{d}(\bm{\theta},n)$ given in (8) clearly satisfies the requirements of Assumption D. Moreover, under the form of the bias postulated in (8), we have that $\beta,\gamma\geq 1$. If the initial estimator is $\sqrt{n}$-consistent, we have $\alpha=\nicefrac{{1}}{{2}}$, and therefore the requirements of Assumptions C and D, i.e.

[TABLE]

are satisfied if

[TABLE]

The last part of Assumption D is particularly mild. Indeed, it simply requires that $\mathbf{c}(n)$, i.e. the part of bias that only depends on the sample size $n$, goes to $\mathbf{0}$ faster than $\sqrt{p}$. For the vast majority of estimators the function $\mathbf{c}(n)$ is either very small or equal to $\mathbf{0}$ as the bias generally depends on $\bm{\theta}$.

The assumptions presented in this section are mild and likely satisfied in most practical situations. Moreover, compared to conditions considered in many instances (see e.g. Guerrier et al., 2019 and the references therein), our assumption framework allows to relax various requirements: (i) the initial estimator may have an asymptotic bias that depends on $\bm{\theta}$; (ii) the initial estimator may also be discontinuous in $\bm{\theta}$ (as is commonly the case when considering discrete data models); (iii) the finite sample bias of the initial estimator is allowed to have a more general expression with respect to previously defined bias functions; (iv) some of the technical requirements (for example on the topology of $\bm{\Theta}$) are relaxed (v) $p$ is allowed to increase with the sample size $n$. Nonetheless, our assumption framework is not necessarily the weakest possible in theory and may be further relaxed (as discussed in Remark E of Appendix H). However, we do not attempt to pursue the weakest possible conditions to avoid overly technical treatments in establishing the theoretical results of the following section.

3.   Main results

In this section, we study the convergence of the IB sequence along with the properties of the IB-estimator. In particular, in Section 3.1 we discuss the convergence of the IB sequence together with the consistency of the IB-estimator. In Section 3.2, under the conditions set in Section 2, we show that the IB-estimator is PT-unbiased while the asymptotic distribution of the estimator is discussed in Section 3.3. The asymptotic results presented in this section are somewhat unusual as we always consider arbitrarily large but finite $n$ and $p$ and our results may not be valid when taking the limit in $n$ (and therefore in $p$). Indeed, the (usual) notion of limit we are considering in this article comes from the usual topology (and its induced metric) of ${\rm I\!R}^{p}$ for finite $p$. In order to ensure that these results hold when $p$ is $\infty$, a more detailed topological discussion is needed. However, the difference between the considered framework and others where limits are studied is rather subtle. A more detailed discussion on our asymptotic framework is provided in Appendix C.

3.1.   Consistency

Theorem 1 below shows that the IB sequence converges (exponentially fast) to the IB-estimator and that the latter is identifiable and consistent.

Theorem 1:

Under Assumptions A, B, C and D, for all $H\in\mathbb{N}^{\ast}$,

There exists a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{\ast}$ with $n\geq n^{*}$, $\left\{\hat{\bm{\theta}}\right\}=\widehat{\bm{\Theta}}$, i.e. the set $\widehat{\bm{\Theta}}$ is a singleton. 2. 2.

There exists a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{\ast}$ with $n\geq n^{*}$, the sequence $\left\{\hat{\bm{\theta}}^{(k)}\right\}_{k\in\mathbb{N}}$ has the following limit

[TABLE]

Moreover, there exists a real $\epsilon\in(0,\,1)$ such that for any $k\in\mathbb{N}^{\ast}$

[TABLE] 3. 3.

$\hat{\bm{\theta}}$* is a consistent estimator of $\bm{\theta}_{0}$, i.e., \big{\lVert}\hat{\bm{\theta}}-\bm{\theta}_{0}\big{\rVert}_{2}=o_{\rm p}(1).*

The proof of Theorem 1 is given in Appendix H. The three results of this theorem are derived separately. Indeed, the first and the second result of the theorem, which correspond to Lemma 2 and Proposition 1 in Appendix D, are proved under a certain set of assumptions that is different from the set of assumptions needed to prove consistency (third result of the theorem). In addition, consistency is proved in two different ways, using two different sets of assumptions (see Proposition 2 and Corollary 1 in Appendices E and F, respectively). As a result, Assumptions A to D condense all the requirements from the previously mentioned assumptions that are found in the appendices, removing possible redundancies between them. A precise account of the dependence structure between the assumptions is given in Figure 4 of Appendix H.

An important consequence of Theorem 1 is that the IB provides a computationally efficient algorithm to solve the optimization problem in (2). Moreover, in practical settings, the IB can often be applied to the estimation of complex models where standard optimization procedures used to solve (2) may fail to converge numerically (see e.g. Guerrier et al., 2019). In practice, the IB procedure is computationally efficient, which is in line with Theorem 1 showing that $\hat{\bm{\theta}}^{(k)}$ converges to $\hat{\bm{\theta}}$ (in norm) at an exponential rate. Even though the convergence of the algorithm may be slower when $p$ is large, in practical situations, the number of iterations necessary to reach a suitable neighbourhood of the solution appears to be relatively small. For example, if we define $k_{n}\equiv\mathcal{O}(\log(pn))$, then we have

[TABLE]

and therefore $\hat{\bm{\theta}}^{(k_{n}^{*})}$ is also a consistent estimator for $\bm{\theta}_{0}$, where $k_{n}^{*}\equiv\lfloor k_{n}\rfloor$.

3.2.   Phase Transition Unbiasedness

In the previous section, we studied the consistency of the IB-estimator $\hat{\bm{\theta}}$ and provided a computationally efficient strategy to obtain it. We now investigate the bias of $\hat{\bm{\theta}}$ and we show that this estimator achieves PT-unbiasedness. This result, combined with the guarantee of consistency, is of particular interest in many practical settings. To obtain this property, additional requirements are needed on the function $\mathbf{d}(\bm{\theta},n)$. For this reason, we introduce Assumption D∗ which combines Assumption D and these additional requirements.

Assumption D∗:

The bias function $\mathbf{d}\left(\bm{\theta},n\right)$ is such that:

The asymptotic bias function $\mathbf{a}(\bm{\theta})$ can be written as

[TABLE]

where $\mathbf{M}\in{\rm I\!R}^{p\times p}$ with $||\mathbf{M}||_{F}<1$ and $\mathbf{s}\in{\rm I\!R}^{p}$. 2. 2.

There exists a $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$, the matrix $(\mathbf{M}+\mathbf{L}(n)+\mathbf{I})^{-1}$ exists. 3. 3.

There exist real $\beta,\gamma>0$ such that for all $\bm{\theta}\in\bm{\Theta}$ and any $j,l=1,\dots,p$, we have

[TABLE] 4. 4.

The Jacobian matrix $\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}(\bm{\theta},n)\in{\rm I\!R}^{p\times p}$ exists and is continuous in ${\bm{\theta}}\in\bm{\Theta}$ for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$. 5. 5.

Defining $c_{n}\equiv\displaystyle{\max_{j=1,\dots,p}}\mathbf{c}_{j}(n)$ for all $n\in\mathbb{N}^{*}$, the sequence $\left\{c_{n}\right\}_{n\in\mathbb{N}^{*}}$ is such that

[TABLE]

To avoid redundancy, we only discuss the additional requirement of Assumption D∗. The first one requires $\mathbf{a}({\bm{\theta}})$ to be a sublinear function of ${\bm{\theta}}$ with its linear term $\mathbf{M}\in{\rm I\!R}^{p\times p}$ satisfying $||\mathbf{M}||_{F}<1$. As discussed after Assumption D, this condition implies that $\mathbf{a}({\bm{\theta}})$ is a contraction map. Moreover, the PT-unbiasedness of $\hat{\bm{\theta}}$ is still guaranteed even if the condition $||\mathbf{M}||_{F}<1$ is not satisfied but the IB sequence may not converge. More details are given in Appendices E and F.

The second additional assumption requires that the matrix $(\mathbf{M}+\mathbf{L}(n)+\mathbf{I})^{-1}$ exists when $n$ is sufficiently large. This requirement is quite general and is, for example, satisfied if $\mathbf{a}({\bm{\theta}})$ is relatively “small” compared to $\bm{\theta}$ or if $\hat{\bm{\pi}}(\bm{\theta}_{0},n)$ is a consistent estimator of $\bm{\theta}_{0}$. Interestingly, this part of the assumption can be interpreted as requiring that the matrix $(\mathbf{M}+\mathbf{I})^{-1}$ exists, which directly implies that the binding function $\bm{\pi}({\bm{\theta}})$ is injective.

The third one requires that $\displaystyle\lim_{n\to\infty}p^{2}n^{-\gamma}=0$. In many practical settings it is reasonable111For example using the bias function proposed in (8) and assuming that the first term of the expansion $\mathbf{h}^{(1)}({\bm{\theta}})$ is linear. to assume that $\gamma=2$ which implies that this condition is satisfied if $\displaystyle\lim_{n\to\infty}p/n=0$. This latter condition is less strict than (9). It is also easy to see that $\displaystyle\lim_{n\to\infty}p^{2}n^{-\gamma}=0$ implies the similar requirement used in Assumption D, namely $\displaystyle\lim_{n\to\infty}\;p^{\nicefrac{{1}}{{2}}}n^{-\gamma}=0$.

Finally, the last additional requirement implies that $\mathbf{R}(\bm{\theta},n)$ is continuous in ${\bm{\theta}}\in\bm{\Theta}$ when $n$ is large enough. A useful consequence of this requirement combined with the compactness of $\bm{\Theta}$ (Assumption A) is that $\mathbf{r}_{j}(\hat{\bm{\theta}},n)$ and $\mathbf{R}_{j,l}(\hat{\bm{\theta}},n)$ are bounded random variables for all $j,l=1,\dots,p$.

Before stating Theorem 2 which ensures that $\hat{\bm{\theta}}$ is a PT-unbiased estimator, we introduce a particular asymptotic notation and state a lemma that provides a strategy for proving PT-unbiasedness. Let $\mathcal{A}\subset\mathbb{N}$ and let $f(n)$ and $g(n)$ be real-valued functions with $g(n)$ being strictly positive for all $n\in\mathbb{N}^{\ast}$. We write $f(n)~{}=~{}\mathcal{O}_{\delta\in\mathcal{A}}\left(g(n)^{\delta}\right)$ if and only if :

[TABLE]

When $\mathcal{A}$ is a finite set we have

[TABLE]

but this equivalence is not preserved when $\mathcal{A}$ is an infinite set. Indeed, the notation

[TABLE]

means that

[TABLE]

Therefore, this definition does not guarantee the existence of a $n^{\ast}\in\mathbb{N}^{*}$ that bounds the sequence $\left\{n_{\delta}\right\}_{\delta\in\mathcal{A}}$ as it may diverge. However, if $\mathcal{A}\subset\mathbb{N}^{\ast}$ is an infinite set, the following equivalence holds true

[TABLE]

Lemma 1 below shows how this particular notation is relevant to show PT-unbiasedness. The proof of this lemma is given below the statement as it is short and provides insight on this property.

Lemma 1:

*Let $g(n)$ be a strictly positive real-valued function such that $\displaystyle\lim_{n\to\infty}g(n)=0$. If

$f(n)~{}=~{}\mathcal{O}_{\delta\in\mathbb{N}}\left(g(n)^{\delta}\right)$ then there exists a $n^{\ast}\in\mathbb{N}^{\ast}$ such that $f(n)=0$ for all $n\geq n^{\ast}$.*

Proof: Let $\delta\in\mathbb{N}$. Since $f(n)=\mathcal{O}_{\delta\in\mathbb{N}}\left(g(n)^{\delta}\right)$ there exists $n^{\ast}\geq 0$ and $M_{\delta}>0$ such that $|f(n)|\leq M_{\delta}\,g(n)^{\delta}$ for all $n\geq n^{\ast}$. We can therefore consider the following quantity

[TABLE]

Defining $\bar{M}_{\delta+1}$ in a similar fashion, we have that for all $n\geq n^{\ast}$ (where $n^{\ast}$ is the max of the ones in the definitions of $\bar{M}_{\delta}$ and $\bar{M}_{\delta+1}$),

[TABLE]

By definition of $\bar{M}_{\delta}$, we have

[TABLE]

for all $n\geq n^{\ast}$, which implies that $\bar{M}_{\delta}=0$. Indeed, supposing the contrary of (12), there exists $n_{0}\geq n^{\ast}$ such that $\bar{M}_{\delta}>\bar{M}_{\delta+1}\,g(n_{0})$. Without loss of generality we can suppose that the function $g(n)$ is decreasing222Indeed, we can define a decreasing step function $\tilde{g}(n)$ such that $\displaystyle\lim_{n\to 0}\tilde{g}(n)=0$ and $g(n)\leq\tilde{g}(n)$, $\forall n\in\mathbb{N}^{\ast}$.. Therefore, setting $M^{\prime}_{\delta}\equiv\bar{M}_{\delta+1}\,g(n_{0})>0$ and redefining $n^{\ast}$ to be $n_{0}$, we have that for all $n\geq n^{\ast}$

[TABLE]

which contradicts the minimality of $\bar{M}_{\delta}$. Since $\delta$ was chosen arbitrarily in $\mathbb{N}$, we have that $\bar{M}_{\delta}=0$ for all $\delta\in\mathbb{N}$, which concludes the proof. ∎

The conditions of Lemma 1 may be considered as being too strong per se, but as the following two examples show, these conditions are actually necessary.

Suppose $\displaystyle\lim_{n\to\infty}g(n)=0$ and $f(n)~{}=~{}\mathcal{O}\left(g(n)^{\delta}\right)$ for all $\delta\in\mathbb{N}$.

In this case, we can consider $g(n)\equiv n^{-1}$ and $f(n)\equiv\exp(-n)$ for all $n\in\mathbb{N}^{\ast}$. Clearly, $f(n)~{}=~{}\mathcal{O}\left(g(n)^{\delta}\right)$ for all $\delta\in\mathbb{N}$ but $f(n)\neq 0$ for all $n\in\mathbb{N}^{\ast}$. 2. 2.

Suppose $\displaystyle\lim_{n\to\infty}g(n)\neq 0$ and $f(n)~{}=~{}\mathcal{O}_{\delta\in\mathbb{N}}\left(g(n)^{\delta}\right)$.

In this case, we consider a decreasing function $g(n)$ such that $0<\displaystyle\lim_{n\to\infty}g(n)=c<1$ and we define $f(n)\equiv d>0$ for all $n\in\mathbb{N}^{\ast}$. Since $g(n)$ is decreasing and $\displaystyle\lim_{n\to\infty}g(n)<1$, there exists $n^{\ast}\in\mathbb{N}$ such that $0<g(n)<1$ for all $n\geq n^{\ast}$. Setting $M_{\delta}\equiv dc^{-\delta}$ for all $\delta\in\mathbb{N}$, we have for all $n\geq n^{\ast}$

[TABLE]

Therefore, $f(n)~{}=~{}\mathcal{O}_{\delta\in\mathbb{N}}\left(g(n)^{\delta}\right)$ but $f(n)\neq 0$ for all $n\in\mathbb{N}^{\ast}$.

Theorem 2 (below) delivers the conditions that guarantee $\hat{\bm{\theta}}$ is a consistent and a PT-unbiased estimator of ${\bm{\theta}}_{0}$.

Theorem 2:

Under Assumptions A, B, C and D∗, and for all $H\in\mathbb{N}^{*}$, the IB-estimator $\hat{\bm{\theta}}$ is consistent and PT-unbiased. Therefore:

$\hat{\bm{\theta}}$* is a consistent estimator of $\bm{\theta}_{0}$, i.e., \big{\lVert}\hat{\bm{\theta}}-\bm{\theta}_{0}\big{\rVert}_{2}=o_{\rm p}(1).* 2. 2.

There exists $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ with $n\geq n^{*}$, we have \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0.

The proof of this result is given in Appendix H where we show

[TABLE]

By Assumption D∗ we have that $p^{2}n^{-\gamma}\to 0$ as $n\to\infty$ and therefore Lemma 1 can be applied to (13) to conclude that \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0 for all $n\geq n^{\ast}$. In practice, as illustrated within the simulation settings in Section 4, the value $n^{*}$ appears to be quite small as the IB-estimator does indeed appear to be unbiased.

In Appendix F we demonstrate the PT-unbiasedness result presented in Theorem 2 under weaker assumptions. With this set of weaker assumptions the consistency of the estimator is not guaranteed. For example, the rate of convergence of the initial estimator (discussed in Assumption C) has to be considered to establish consistency results (in particular when $p$ diverges) but is not needed to study the bias of $\hat{\bm{\theta}}$.

Finally, Theorem 2 relies on the additional assumption that the asymptotic bias $\mathbf{a}({\bm{\theta}})$ is sublinear, which can be quite restrictive. However, if $\hat{{\bm{\pi}}}(\bm{\theta}_{0},n)$ is a consistent estimator of $\bm{\theta}_{0}$, then this restriction is automatically satisfied. Nevertheless, considering the third part of Theorem 1, a possible approach that would guarantee a PT-unbiased estimator is to obtain the IB-estimator $\hat{\bm{\theta}}$ from an inconsistent initial estimator $\hat{\bm{\pi}}(\bm{\theta},n)$ and then to compute a new IB-estimator with the initial estimator $\hat{\bm{\pi}}^{\text{NEW}}(\bm{\theta}_{0},n)\equiv\hat{\bm{\theta}}$. In practice, however, this computationally intensive approach is probably unnecessary since the (one step) IB-estimator appears to eliminate the bias almost completely when considering inconsistent initial estimators $\hat{\bm{\pi}}(\bm{\theta},n)$ (see e.g. Section 4.1). In fact, we conjecture that for a larger class of asymptotic bias functions, Theorem 2 remains true, although the verification of this conjecture is left for further research.

3.3.   Asymptotic Normality

In this section we study the asymptotic distribution of the IB-estimator. In the case where $p\to c<\infty$ as $n\to\infty$, the results of Gourieroux et al. (1993) can be directly applied to obtain

[TABLE]

where

[TABLE]

and where

[TABLE]

Apart from technical requirements (such as the continuity of $\mathbf{B}(\bm{\theta})$ in $\bm{\theta}_{0}$), the main condition ensuring the validity of (14) is the asymptotic normality of $\mathbf{v}(\bm{\theta}_{0},n)$ (or $\hat{\bm{\pi}}(\bm{\theta}_{0},n)$), i.e.

[TABLE]

However, in the case where $p\to\infty$ as $n\to\infty$, such a condition requires to be considerably modified by assuming instead that $\mathbf{v}(\bm{\theta}_{0},n)$ satisfies a specific Gaussian approximation. Moreover, several other additional technical requirements are needed to adapt (14) to the case where $p$ diverges. To avoid an overly technical discussion in the main text, these conditions are stated and discussed in Appendix G together with Theorem 3 in Appendix H which proposes a possible high dimensional counterpart of (14). Unfortunately, many of the conditions needed to derive these results can be strong and difficult to verify for a specific model and initial estimator. In the case where $p\to c<\infty$ as $n\to\infty$, Theorem 3 becomes equivalent to (14) and relies on very similar requirements as those needed for the asymptotic results in Gourieroux et al. (1993). Informally, Theorem 3 states that for all $\mathbf{u}\in{\rm I\!R}^{p}$ such that $||\mathbf{u}||_{2}=1$, we have

[TABLE]

where

[TABLE]

One interesting difference in the conditions needed to derive (14) and (15), is that the former is valid for all $H\in\mathbb{N}^{*}$ while the latter requires $H=\mathcal{O}(p^{4})$. Therefore, it suggests that the “quality” (and validity) of the approximation depends on $H$ when $p$ diverges. However, the conditions of Theorem 3 are sufficient but may not be necessary thereby implying that $H=\mathcal{O}(p^{4})$ may not always be needed as a condition for this approximation.

In practice, the estimation of the variance of the IB-estimator can be obtained through different methods. A simple approach takes advantage of the results of the last iteration of the IB sequence to construct a parametric bootstrap estimator of the variance of the initial estimator $\hat{\bm{\pi}}(\bm{\theta}_{0},n)$. Indeed, in the last iteration of the IB sequence, $H$ samples are simulated under $\hat{\bm{\theta}}$ (which is a consistent estimator), which allows to compute the following quantity (assuming $H$ to be sufficiently large)

[TABLE]

where $\bar{\bm{\pi}}(\hat{\bm{\theta}},n)\equiv\frac{1}{H}\sum_{h=1}^{H}\hat{\bm{\pi}}^{*}_{h}(\hat{\bm{\theta}},n)$. Then, an estimator of the covariance matrix of $\hat{\bm{\theta}}$ can be obtained as follows

[TABLE]

where $\widehat{\mathbf{B}}$ can be obtained by numerical derivation of $\hat{\bm{\pi}}(\bm{\theta},n)$ evaluated at $\hat{\bm{\theta}}$.

4.   Application: Logistic Regression Model

In this section, we apply the methodology developed in Sections 2 and 3 to investigate the performance of IB-estimators in three different settings. As previously mentioned, our conditions allow the initial estimator to be discontinuous in $\bm{\theta}$ and it is therefore of interest to consider the logistic regression model, which may be the most commonly used model for binary (response) data. To illustrate the flexibility of the IB estimation approach, as initial estimators we select slightly modified versions of the MLE and of a robust estimator. As explained further on, these modifications were introduced to allow the estimators to be “robust” to the problem of data separation. To compare these IB-estimators, we consider the MLE, the robust estimator as well as the bias corrected estimator proposed by Kosmidis and Firth (2009). Their respective performance is studied when the data is generated at the model but also under slight model misspecification where “outliers” are randomly created. Then, we extend the logistic regression model to include a random intercept, a special case of Generalized Linear Mixed Models (GLMM) for which there is no closed form expression for the likelihood function. In this case, the initial estimator is selected to be a numerically simple approximation to the MLE, and its finite sample behaviour is compared to the ones of the MLE computed using several (more precise) approximation methods.

4.1.   Bias Corrected Estimators for the Logistic Regression Model

4.1.1.   Introduction

The Logistic Regression Model (LRM) (Nelder and Wedderburn, 1972; McCullagh and Nelder, 1989) is one of the most frequently used models for binary response variables conditioned on a set of covariates. However, it is well known that in some quite frequent practical situations, the MLE is biased or its computation can become very unstable, especially when performing some type of resampling scheme for inference. The underlying reasons are diverse, but the main ones are the possibly large $p/n$ ratio, separability (leading to regression slope estimates of infinite value) and data contamination (robustness).

The first two sources are often confounded and practical solutions are continuously sought to overcome the difficulty in performing “reasonable” inference. For example, in medical studies, the bias of the MLE together with the problem of separability has led to a rule of thumb called the number of Events Per Variable (EPV), that is the number of occurrences of the least frequent event over the number of predictors, which is used in practice to choose the maximal number of predictors one is “allowed” to use in a LRM (see e.g. Austin and Steyerberg, 2017 and the references therein ).

The problem of separation or near separation in logistic regression is linked to the existence of the MLE which is not always guaranteed (see Silvapulle, 1981; Albert and Anderson, 1984; Santner and Duffy, 1986; Zeng, 2017). Alternative conditions for the existence of the MLE have been recently developed in Candès and Sur (2019) (see also Sur and Candès, 2019). In order to detect separation, several approaches have been proposed (see for instance, Lesaffre and Albert, 1989; Kolassa, 1997; Christmann and Rousseeuw, 2001). The adjustment to the score function proposed by Firth (1993) has been implemented in Kosmidis and Firth (2009) (see also Kosmidis and Firth, 2010) to propose a bias corrected MLE for GLM which has the additional natural property that it is not subject to the problem of separability . Moreover, the MLE is known to be sensitive to slight model deviations that take the form of outliers in the data leading to the proposal of several robust estimators for the LRM and more generally for GLM (see e.g. Cantoni and Ronchetti, 2001; Cíźek, 2008; Heritier et al., 2009 and the references therein .

Despite all the proposals for finite sample bias correction, separation and data contamination problems, no estimator has so far been able to handle the three potential sources of bias jointly. In this section, we make use of the IB-estimator which is built through a simple adaptation of available estimators. Although the latter estimators might possibly not be the best ones for this problem, they at least guarantee, at a reasonable computational cost, a reduced finite sample bias which is comparable, for example, to that of the bias reduced MLE of Kosmidis and Firth (2009) in uncontaminated data settings as well as reducing the bias in contaminated data settings. Moreover, in both cases, we adapt the initial estimator so that it is not affected by the problem of separability.

4.1.2.   Bias Corrected Estimators

Consider the LRM with response $\mathbf{Y}(\bm{\beta}_{0},n)$ and linear predictor $\mathbf{X}\bm{\beta}$, where $\mathbf{X}$ is an $n\times p$ matrix of fixed covariates with rows $\mathbf{x}_{i},i=1,\ldots,n$, and with logit link $\mu_{i}(\bm{\beta})\equiv\mathbb{E}[\mathbf{Y}_{i}(\bm{\beta},n)]=\exp(\mathbf{x}_{i}\bm{\beta})/(1+\exp(\mathbf{x}_{i}\bm{\beta}))$. The MLE for ${\bm{\beta}}$ is given by

[TABLE]

and can be used as an initial estimator in order to obtain the IB-estimator in (1).

To avoid the (potential) problem of separation, we follow the suggestion of Rousseeuw and Christmann (2003) to transform the observed responses $\mathbf{Y}_{i}(\bm{\beta}_{0},n)$ to get pseudo-values

[TABLE]

for all $i=1,\dots,n$, where $\delta\in\left[0,0.5\right)$ is a (fixed) “small” (i.e. close to 0) scalar. For a discussion on the choice of $\delta$ and also possible asymmetric transformations, see Rousseeuw and Christmann (2003). Jointly using these two approaches (i.e. MLE in (16) computed on the pseudo-values) as an initial estimator, we denote the estimator resulting from the IB procedure as the IB-MLE.

As a robust initial estimator, we consider the robust $M$-estimator proposed by Cantoni and Ronchetti (2001), with general estimating function (for GLMs) given by

[TABLE]

with $r\left({\bm{\beta}},\mathbf{Y}_{i}(\bm{\beta}_{0},n)\right)\equiv\left(\mathbf{Y}_{i}(\bm{\beta}_{0},n)-\mu_{i}({\bm{\beta}})\right)V^{-1/2}\left(\mu_{i}({\bm{\beta}})\right)$ being the Pearson residuals and with consistency correction factor

[TABLE]

where the expectation is taken over the (conditional) distribution of the responses $\mathbf{Y}_{i}(\bm{\beta}_{0},n)$ (given $\mathbf{x}_{i}$). For the LRM, we have $V\left(\mu_{i}({\bm{\beta}})\right)\equiv\mu_{i}({\bm{\beta}})(1-\mu_{i}({\bm{\beta}}))$. We compute the robust initial estimator $\hat{{\bm{\pi}}}(\bm{\beta}_{0},n)$ on the pseudo-values (17), using the implementation in the glmrob function of the robustbase package in R (Maechler et al., 2019), with $\psi_{c}$ in (18) being the Huber loss function (with default parameter $c$) (see Huber, 1964) and $w\left(\mathbf{x}_{i}\right)=\sqrt{1-h_{ii}}$, $h_{ii}$ being the $i$th diagonal element of $\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-1}\mathbf{X}^{T}$. The resulting robust estimator is taken as the initial estimator in (1) to get a robust IB-estimator that we call IB-ROB. The resulting IB-estimator is in fact robust in the sense that it has a bounded influence function (see Hampel et al., 1986; Hampel, 1974).

Both initial estimators are not consistent estimators, even if we expect the asymptotic bias to be very small (for small values of $\delta$ in (17)), but both IB-estimators should have a reduced finite sample bias where, in addition, the IB-ROB is also robust to data contamination.

4.1.3.   Simulation Study

We perform a simulation study to validate the properties of the IB-MLE and IB-ROB and compare their finite sample performance to other well established estimators. In particular, as a benchmark, we also compute the MLE, the bias reduced MLE (MLE-BR) using the brglm function (with default parameters) of the brglm package in R (Kosmidis, 2019), as well as the robust estimator (18) using the glmrob function in R without data transformation (ROB). We consider four situations that can occur with real data which result from the combinations of balanced outcome classes (Setting I) and unbalanced outcome classes (Setting II) with and without data contamination. We also consider a large model with $p=200$ and chose $n$ as to provide EPV of respectively $5$ and $3.75$, which are below the usually recommended value of 10. The parameter values for the simulations are provided in Table 1.

The covariates were simulated independently from distributions $\mathcal{N}(0,4/\sqrt{n})$ for Setting I and $\mathcal{N}(0.6,4/\sqrt{n})$ for Setting II, in order to ensure that the size of the log-odds ratio $\mathbf{x}_{i}{\bm{\beta}}$ does not increase with $n$, so that $\mu_{i}({\bm{\beta}})$ is not trivially equal to either [math] or $1$. To contaminate the data, we chose a rather extreme misclassification error to show a noticeable effect on the different estimators, which consists in permuting 2% of the responses with corresponding larger (smaller) fitted probabilities (expectations). The simulation results are presented in Figure 1 as boxplots of the finite sample distribution, and in Figure 2 as the bias and Root Mean Squared Error (RMSE) of the different estimators.

The finite sample distributions presented in Figure 1, as well as the summary statistics given by the bias and RMSE presented in Figure 2, allow us to draw the following conclusions that support the theoretical results. In the uncontaminated case, as expected the MLE is biased (except when the slope parameters are zero) and the MLE-BR, IB-MLE and IB-ROB are all unbiased which is not the case for the robust estimator ROB. Moreover, the variability of all estimators is comparable, except for ROB which makes it rather inefficient in these settings. With 2% of contaminated data (missclassification error), the only unbiased estimator is IB-ROB and its behaviour remains stable compared to the uncontaminated data setting. This is in line with a desirable property of robust estimators, that is stability with or without (slight) data contamination. The behaviour of all estimators remains the same in both settings, that is, whether or not the responses are balanced. Finally, as argued above, a better proposal for a robust, bias reduced and consistent estimator, as an alternative to IB-ROB, could in principle be proposed, but this is left for further research.

4.2.   Bias Corrected Estimator for the Random Intercept Logistic Regression Model

An interesting way of extending the LRM to account for the dependence structure between the observed responses is to use the GLMM family (see e.g. Nelder and Wedderburn, 1972; Breslow and Clayton, 1993; Lee and Nelder, 2001; McCulloch and Searle, 2001; Jiang, 2007). We consider here the case of the random intercept model, which is a suitable model in many applied settings. The binary response is denoted by $\mathbf{Y}_{ij}(\bm{\beta}_{0},n)$, where $i=1,\ldots,m$ and $j=1,\ldots,n_{i}$. The expected value of the response is expressed as

[TABLE]

where $\beta_{0}$ is the intercept, $\mathbf{x}_{ij}$ is a $q$-vector of covariates, $\bm{\beta}$ is a $q$-vector of regression coefficients and the random effect $U_{i},\;i=1,...,m$ is a normal random variable with zero mean and (unknown) variance $\sigma^{2}$.

Because the random effects are not observed, the MLE is derived on the marginal likelihood function where the random effects are integrated out. These integrals have no known closed-form solutions, so approximations to the (marginal) likelihood function have been proposed, including the pseudo- and Penalized Quasi- Likelihood (PQL) (see Schall, 1991; Wolfinger and O’connell, 1993; Breslow and Clayton, 1993), LAplace approximations (LA) (Raudenbush et al., 2000; Huber et al., 2004) and adaptive Gauss-Hermite Quadrature (GHQ) (Pinheiro and Chao, 2006). It is well known that PQL methods lead to biased estimators while LA and GHQ are more accurate (for extensive accounts of methods across software and packages, see e.g. Bolker et al., 2009; Kim et al., 2013). The lme4 R package (Bates et al., 2015) uses both the LA and GHQ to compute the likelihood function, while the glmmPQL R function of the MASS library (that supports Venables and Ripley, 2002) uses the PQL.

In this section, we propose an IB-estimator that has a reduced finite sample bias compared to the (approximated) MLE using either the LA, GHQ or PQL approaches. More precisely, as an initial estimator, for computational efficiency, we consider an estimator defined through penalized iteratively reweighted least squares as described, for example, in Bates et al. (2015) and implemented in the glmer function (with argument nAGQ set to 0) of the lme4 R package. As done in Section 4.1.2 for the LRM, we transform the observed responses to get pseudo-values as in (17) with $\delta=0.01$. Hence, the initial estimator is a less accurate approximation to the MLE (compared e.g. to the LA or the GHQ approximations) although it is not affected by data separation which makes it asymptotically biased. However, we expect the IB to correct the bias induced through the transformed responses.

To study the behaviour of the IB-estimator and compare its performance in terms of bias and variance in finite samples to different approximations of the MLE, we perform a simulation study using the two settings described in Table 2. Both settings can be considered as high dimensional in the sense that $p$ is relatively large with respect to $m$. Moreover, while Setting II ($m=50$, $n_{i}=n_{\circ}=5,\forall i$) reflects a possibly more frequent situation, Setting I concerns the case where $m$ is small (and much smaller than $n_{\circ}$), a situation frequently encountered in cluster randomised trials (see e.g. Huang et al., 2016; Leyrat et al., 2017, and the references therein). As for the simulation study in Section 4.1.3 studying the LRM, the covariates are simulated independently from the distribution $\mathcal{N}(0,4/\sqrt{n})$.

Figure 3 presents the finite sample bias and RMSE of the approximated MLE estimators and the IB-estimator. It would appear evident that the proposed IB-estimator drastically reduces (removes) the finite sample estimation bias, especially for the random effect variance estimator. Moreover, the IB-estimator also achieves the lowest RMSE. These simulation results confirm the theoretical properties of the IB-estimator, namely finite sample (near) unbiasedness, which has important advantages when performing inference in practice, and/or when the parameter estimates, such as the random intercept variance estimate, are used, for example, to evaluate the sample size needed in subsequent randomized trials. The IB-estimator can eventually be based on another initial estimators in order to possibly improve efficiency for example, however this study is left for future research.

Appendix A Notation and Organisation of the Appendices

As mentioned in Remark A, in order to avoid potential confusing situations, we will use more precise notations in the sequel. Namely, we make our notations more precise and we use, as presented in Table 3, a dictionary between the notations used in the main text and in the appendices. The justification for these precision in notation is given in Appendix B.

Another notable difference between the main text and the appendices is that the results presented in the latter are proved independently in the following appendices. Indeed, in Appendix D we prove the converenge of the IB sequence, in Appendix F the unbiasedness of the IB-estimator, in Appendix E its consistency and finally in Appendix G we state and prove its asymptotic normality. In each Appendix, the assumptions are tailored to the result we aim to prove resulting in a consequent number of assumptions. For this reason, an assumption framework is defined in Appendix C to classify assumptions of the same type. In the last one, Appendix H, these assumptions are combined to prove the main results of the article, namely Theorem 1 and Theorem 2 that are stated in the main text, as well as Theorem 3 which formally states the asymptotic normality of the IB-estimator.

Although the results presented in Appendices D to F are closely related to the ones of Theorem 1 and 2, we nevertheless provide a detailed discussion that is, in parts, redundant, but makes it complete within the material provided in these appendices.

Appendix B Mathematical Setup

In this appendix, we redefine the mathematical setup presented in the main text in order to make it more precise, by using the notations presented in Table 3. Let us define $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}\right)\in{\rm I\!R}^{n}$ as being a random sample generated under model $F_{\bm{\theta}}$ (possibly conditional on a set of fixed covariates), where $\bm{\theta}\in\bm{\Theta}\subset{\rm I\!R}^{p}$ is the parameter vector of interest and $\bm{\omega}\in{\rm I\!R}^{m},\,m\geq n>0$, represents a random variable explaining the source of randomness of the sample. More specifically, $\bm{\omega}$ can be considered as a random seed that produces, conditionally on its value and given a value of $\bm{\theta}$, the sample $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}\right)$ of size $n$ in a deterministic manner. Indeed, $\bm{\omega}$ can be conceived as a random variable issued from a model $G$, thereby justifying the definition of the random sample as $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}\right)$. With this in mind, denoting “$\overset{d}{=}$” as “equality in distribution”, there is no requirement for $\bm{\omega}$ to be unique since it is possible to have $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}\right)\overset{d}{=}\mathbf{X}\left(\bm{\theta},n,\bm{\omega}^{\ast}\right)$ even though $\bm{\omega}\overset{d}{\neq}\bm{\omega}^{\ast}$. For simplicity, throughout this article we will consider $\bm{\omega}$ as only belonging to a fixed set ${\bm{\Omega}}\equiv\left\{\bm{\omega}_{h}\right\}_{h=1,\dots,H}\subset{\rm I\!R}^{m}$. In this setting, defining $\bm{\theta}_{0}\in\bm{\Theta}$, we let $\bm{\omega}_{0}$ and $\bm{\theta}_{0}$ denote respectively the unknown random vector and fixed parameter vector used to generate the random sample $\mathbf{X}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)$ that will be used to estimate $\bm{\theta}_{0}$. Knowing that the sample $\mathbf{X}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)$ is generated from a random seed $\bm{\omega}_{0}$, we can also consider other samples that we denote as $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}\right)$, where $(\bm{\theta},n,\bm{\omega})\in\bm{\Theta}\times\mathbb{N}^{\ast}\times\bm{\Omega}$, which can therefore be simulated based on different values of $\bm{\omega}$ (and $\bm{\theta})$. Notice that, based on this premise, we have that $\bm{\omega}_{0}\not\in{\bm{\Omega}}$ which therefore implies that, in the context of this article, it is not possible to generate samples using $\bm{\omega}_{0}$. Hence, the use of $\bm{\omega}$ allows to explicitly display the randomness of a sample. In addition, with this notation it is not only possible to clearly distinguish the observed sample and the simulated ones but also to define the difference (or equivalence) between two simulated samples, say $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}_{j}\right)$ and $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}_{l}\right)$.

Now, we focus on the estimation procedure aimed at targeting the value of $\bm{\theta}_{0}$ from $\mathbf{X}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)$. For this purpose, let us define $\hat{\bm{\pi}}\left(\mathbf{X}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)\right)$ as a biased estimator of $\bm{\theta}_{0}$ which, despite being possibly inconsistent, is either readily available or can easily be computed. For simplicity of notation, from this point onward we will refer to the estimator $\hat{\bm{\pi}}\left(\mathbf{X}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)\right)$ as $\hat{\bm{\pi}}\left(\bm{\theta}_{0},n,\bm{\omega}_{0}\right)$. Similarly, to the definitions of Section 2 we define, for all $(n,H)\in\mathbb{N}^{\ast}\times\mathbb{N}^{\ast}$, the following set

[TABLE]

where $\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega}_{h})$ is obtained on $\mathbf{X}\left(\bm{\theta},n,\bm{\omega}_{h}\right)$. This definition assumes that the solution set $\widehat{\bm{\Theta}}_{(n,H)}$ is not empty for all $(n,H)\in\mathbb{N}^{\ast}\times\mathbb{N}^{\ast}$, which is reasonable when $\dim(\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega}))=\dim(\bm{\theta})$. In similar fashion, the IB sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ can defined using the notation of this appendix as

[TABLE]

where $\hat{\bm{\theta}}_{(n,H)}^{(0)}\in\bm{\Theta}$. When this iterative procedure converges, we define the IB-estimator $\hat{\bm{\theta}}_{(n,H)}\in{\rm I\!R}^{p}$ as the limit in $k$ of $\hat{\bm{\theta}}_{(n,H)}^{(k)}$.

Appendix C Assumptions and Asymptotic Framework

Various assumptions are considered and, in order to provide a general reference for the type of assumption, we employ the following conventions to name them:

• •

Assumptions indexed by “A” refer to the topology of $\bm{\Theta}$;

• •

Assumptions indexed by “B” refer to the existence of the expected value of the initial estimator;

• •

Assumptions indexed by “C” refer to the random component of the initial estimator;

• •

Assumptions indexed by “D” refer to the bias of the initial estimator.

When needed, sub-indices or upper-indices will be used to distinguish assumptions of the same type (for example Assumptions A1 and A are different but both related to the topology of $\bm{\Theta}$).

As it was mentioned in the main text, the asymptotic framework we use here is somewhat unusual as we always consider arbitrarily large but finite $n$ and $p$. Indeed, if $p$ is such that $p\to\infty$ as $n\to\infty$, our asymptotic results may not be valid when taking the limit in $n$ (and therefore in $p$) but are valid for all finite $n$ and $p$. The difference between the considered framework and others where limits are studied is rather subtle as, in practice, asymptotic results are typically used to derive approximations in finite samples for which the infinite (dimensional) case is not necessarily informative. Indeed, the topology of ${\rm I\!R}^{p}$ for any finite $p$ and ${\rm I\!R}^{\infty}$ are profoundly different. For example, the consistency of an estimator is generally dependent on the assumption that $\bm{\Theta}$ is compact. In the infinite dimensional setting closed and bounded sets are not necessarily compact. Therefore, this assumption becomes rather implausible for many statistical models and would imply, among other things, a detailed topological discussion of requirements imposed on $\bm{\Theta}$. Similarly, many of the mathematical arguments presented in this article may not apply in the infinite dimensional case. Naturally, when $p\to c<\infty$ as $n\to\infty$, limits in $n$ are allowed. Although not necessary, we abstain from using statements or definitions such as “$\,\xrightarrow{p}$” (convergence in probability) or “$\,\xrightarrow{d}\,$” (convergence in distribution) which may lead one to believe that the limit in $n$ exists. This choice is simply made to avoid confusion.

Appendix D Converenge of the IB Sequence

In this appendix, we define new assumptions that are weaker than the ones presented in the main text. These assumptions are used to study the properties of the IB sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ defined in (B.2). Our first assumption concerns the topology of the parameter space $\bm{\Theta}$.

Assumption A1:

Let $\bm{\Theta}$ be such that

[TABLE]

Assumption A1 is clearly weaker than Assumption A. It is very mild and essentially ensures that no solution outside of $\bm{\Theta}$ exists for the IB-estimator defined in (B.1). In addition, Assumption A1 doesn’t require $\bm{\Theta}$ to be a compact set since, for example, it can represent an open set. Next, we simply restate Assumption B with the notation used in the appendices.

Assumption B:

For all $\left(\bm{\theta},\,n,\,\bm{\omega}\right)\in\bm{\Theta}\times\mathbb{N}^{\ast}\times{\bm{\Omega}}$, the expectation $\bm{\pi}\left(\bm{\theta},n\right)\equiv\mathbb{E}\left[\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega})\right]$ exists and is finite, i.e. $\lvert\bm{\pi}_{j}\left(\bm{\theta},n\right)\rvert<\infty$ for all $j=1,\ldots,p$. Moreover, for all $j=1,\ldots,p$ $\bm{\pi}_{j}\left(\bm{\theta}\right)\equiv\displaystyle{\lim_{n\to\infty}}\bm{\pi}_{j}\left(\bm{\theta},n\right)$ exists.

Using Assumption B, we can write:

[TABLE]

where $\mathbf{v}\left(\bm{\theta},n,\bm{\omega}\right)\equiv\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega})-\bm{\pi}\left(\bm{\theta},n\right)$ is a zero-mean random vector.

Similarly, we restate Assumption C using the notation of the appendices.

Assumption C:

The random vector $\mathbf{v}(\bm{\theta},n,\bm{\omega})$ is such that its second moment exists and such that there exists a real $\alpha>0$ such that for all $\left(\bm{\theta},\,\bm{\omega}\right)\in\bm{\Theta}\times{\bm{\Omega}}$ and all $j=1,\,\ldots,\,p$, we have

[TABLE]

Next, we adapt the bias terms to the notation used in the appendices. Indeed, we consider the bias of $\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega})$ for $\bm{\theta}$ and we let $\mathbf{d}\left(\bm{\theta},n\right)\equiv\mathbb{E}\left[\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega})\right]-\bm{\theta}$, allowing to write

[TABLE]

The functions $\mathbf{a}\left({\bm{\theta}}\right),\mathbf{b}(\bm{\theta},n),\mathbf{c}(n),\mathbf{L}(n)$ and $\mathbf{r}\left({\bm{\theta}},n\right)$ are defined in the same way as in (5) and (6). Assumption D1 (below) imposes restrictions on the bias function that are weaker than the ones proposed in Assumption D.

Assumption D1:

The bias function $\mathbf{d}\left(\bm{\theta},n\right)$ is such that:

The function $\mathbf{a}(\bm{\theta})$ is a contraction map in that for all $\bm{\theta}_{1},\bm{\theta}_{2}\in\bm{\Theta}$ such that $\bm{\theta}_{1}\neq\bm{\theta}_{2}$ we have

[TABLE] 2. 2.

There exist real $\beta,\gamma>0$ such that for all $\bm{\theta}\in\bm{\Theta}$ and all $j,l=1,\,\ldots,\,p$, we have

[TABLE]

There are two differences between Assumptions D and D1. The first one is that Assumption D requires

[TABLE]

whereas in Assumption D1 this requirement is simply

[TABLE]

This modification of Assumption D implies that under the form of the bias postulated in (8), where we have that $\beta,\gamma\geq 1$ and for $\sqrt{n}$-consistent initial estimator, where we have $\alpha=\nicefrac{{1}}{{2}}$, the requirements of Assumptions C and D1 are satisfied if

[TABLE]

This is a weaker condition than (9), obtained under Assumption D. The second difference between Assumptions D and D1 is that the requirement on the term $\mathbf{c}(n)$ are not needed.

We now study the convergence properties of the IB sequence defined in (B.2) when used to obtained the IB-estimator presented in (B.1). In Lemma 2, we show that when $n$ is sufficiently large, the solution space $\widehat{\bm{\Theta}}_{(n,H)}$ contains only one element. In other words, this result ensures that the function $\frac{1}{H}\sum_{h=1}^{H}\hat{\bm{\pi}}(\bm{\theta},n,\bm{\omega}_{h})$ in (B.1) is injective for fixed (but possibly large) sample size $n$ and fixed $H$. Lemma 2 formally states this result.

Lemma 2:

Under Assumptions A1, B, C and D1, for all $H\in\mathbb{N}^{\ast}$ and $n\in\mathbb{N}^{\ast}$ sufficiently large, the solution space $\widehat{\bm{\Theta}}_{(n,H)}$ is a singleton.

Proof: We define the function

[TABLE]

where

[TABLE]

We define the set of fixed points of $T_{(n,H)}$ as follows

[TABLE]

Clearly, we have that

[TABLE]

as defined in (B.1). In the remainder of the proof we will show that an extension of the function $T_{(n,H)}$ admits a unique fixed point in ${\rm I\!R}^{p}$ by applying Kirszbraun theorem (see Federer, 2014) and Banach fixed-point theorem. Then, using (D.3) and Assumption A1 we will be able to conclude that the set $\widehat{\bm{\Theta}}_{(n,H)}$ only contains this fixed point. Let us start by using (D.1) and Assumptions B, C and D1, we have that

[TABLE]

where $\bm{\delta}_{j}^{(1)}=\mathcal{O}\left(n^{-\gamma}\right)+\mathcal{O}_{\rm p}\left(n^{-\alpha}\right)$ for $j=1,\,\ldots,\,p$ and

[TABLE]

where $\bm{\delta}^{(2)}_{j}=\mathcal{O}\left(n^{-\gamma}\right)+\mathcal{O}_{\rm p}\left(H^{-\nicefrac{{1}}{{2}}}n^{-\alpha}\right)$ for $j=1,\,\ldots,\,p$. Therefore, we have

[TABLE]

Next, we consider the following quantity for $\bm{\theta}_{1},\bm{\theta}_{2}\in\bm{\Theta}$,

[TABLE]

where $\bm{\delta}_{j}^{(3)}=\mathcal{O}\left(n^{-\gamma}\right)+\mathcal{O}_{\rm p}\left(n^{-\alpha}\right)$ for $j=1,\,\ldots,\,p$. Then, using the fact that

[TABLE]

where $\lVert\cdot\rVert_{F}$ is the Frobenius norm, we obtain

[TABLE]

We now consider each term of the above equation separately. First, since $\bm{a}({\bm{\theta}})$ is a contraction map, there exists $\varepsilon\in(0,1)$ such that

[TABLE]

Secondly, we have

[TABLE]

where $\Delta=\mathcal{O}\left(p^{2}n^{-2\beta}\right)$. Indeed, by writing $\mathbf{B}=\mathbf{L}(n)^{T}\mathbf{L}(n)$ we obtain

[TABLE]

Finally, we have

[TABLE]

By combining these results, we have

[TABLE]

Since $\alpha,\beta,\gamma>0$ by Assumptions C and D1 , for sufficiently large $n$ we have that for all $\bm{\theta}_{1},\,\bm{\theta}_{2}\in\bm{\Theta}$

[TABLE]

Using Kirszbraun theorem, we can extend $T_{(n,H)}$ to a contraction map from ${\rm I\!R}^{p}$ to itself. Therefore, applying Banach fixed-point theorem, there exists a unique fixed point $\bm{\theta}^{\ast}_{(n,H)}\in{\rm I\!R}^{p}$ such that

[TABLE]

By Assumption A1 we have that $\bm{\theta}^{\ast}_{(j,n,H)}\in\bm{\Theta}$, which implies by (D.3) that

[TABLE]

which concludes the proof.

∎

Building on the identification result given in the above lemma, the following proposition states that the IB-estimator is the limit of the IB sequence.

Proposition 1:

Under Assumptions A1, B, C and D1, for all $H\in\mathbb{N}^{\ast}$ and $n\in\mathbb{N}^{\ast}$ sufficiently large, the IB sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ has a limit which is the IB-estimator, i.e.

[TABLE]

Moreover, there exists a real $\epsilon\in(0,\,1)$ such that for any $k\in\mathbb{N}^{\ast}$

[TABLE]

Proof: We start by recalling that $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ is defined as

[TABLE]

where $\hat{\bm{\theta}}_{(n,H)}^{(0)}\in\bm{\Theta}$ and the function $T_{(n,H)}$ is defined in (D.2). Using the same arguments used in the proof of Lemma 2 we have that $T_{(n,H)}$ allows to apply Banach fixed-point theorem implying that

[TABLE]

which concludes the first part of the proof.

For the second part, consider a $k\in\mathbb{N}^{\ast}$. By the inequality (D.4) of the proof of Lemma 2, there exists an $\epsilon\in(0,\,1)$ such that

[TABLE]

By construction of the sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ and since $\hat{\bm{\theta}}_{(n,H)}$ is a fixed point of $T_{(n,H)}$, we can derive that

[TABLE]

which concludes the second part of the proof. ∎

Appendix E Consistency of the IB-Estimator

In this appendix, we consider the consistency property of the IB-estimator and set weaker versions of Assumptions A and D, used in Theorem 1.

Assumption A2:

Preserving the same requirement given in Assumption A1, we add the condition that $\bm{\Theta}$ is compact.

Assumption A2 ensures that $\bm{\Theta}$ is compact but remains weaker than Assumption A. We now introduce Assumption D2.

Assumption D2:

Preserving the same definitions given in Assumption D1 and defining $c_{n}\equiv\displaystyle{\max_{j=1,\dots,p}}\mathbf{c}_{j}(n)$, we require that:

$\mathbf{a}({\bm{\theta}})$* is continuous and such that the function ${\bm{\theta}}+\mathbf{a}({\bm{\theta}})$ is injective,* 2. 2.

the constant $\beta$ and the sequence $\left\{c_{n}\right\}_{n\in\mathbb{N}^{*}}$ are such that

[TABLE]

The other requirements of Assumption D1 remain unchanged.

As stated in Assumption D2, some requirements are the same as those in Assumption D1, namely that there exist real $\beta,\gamma>0$ such that for all $\bm{\theta}\in\bm{\Theta}$, we have, for all $j,l=1,\,\ldots,\,p$,

[TABLE]

Compared to Assumptions D and D1, Assumption D2 relaxes the condition on $\mathbf{a}({\bm{\theta}})$. Apart from this difference, Assumptions D2 and D impose the same requirements. Compared to Assumption D1, Assumption D2 imposes stronger requirements on $\beta$ and finally imposes a condition on the vector $\mathbf{c}(n)$. Clearly, Assumption A2 implies Assumption A1 while, on the contrary, Assumptions D1 and D2 don’t mutually imply each other. In the situation where $\alpha=\nicefrac{{1}}{{2}}$, $\beta=1$, $\gamma=2$ and $\mathbf{c}(n)=\mathbf{0}$, Assumption D2 (on which Proposition 2 is based) is satisfied if

[TABLE]

Proposition 2:

Under Assumptions A2, B, C and D2, $\hat{\bm{\theta}}_{(n,H)}$ is a consistent estimator of $\bm{\theta}_{0}$ for all $H\in\mathbb{N}^{\ast}$ in that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have:

[TABLE]

Proof: This proof is directly obtained by verifying the conditions of Theorem 2.1 of Newey and McFadden (1994) on the functions $\hat{Q}(\bm{\theta},n)$ and $Q(\bm{\theta})$ defined as follow:

[TABLE]

and

[TABLE]

Reformulating the requirements of this theorem to our setting, we have to show that (i) $\bm{\Theta}$ is compact, (ii) ${Q}(\bm{\theta})$ is continuous, (iii) ${Q}(\bm{\theta})$ is uniquely minimized at ${\bm{\theta}}_{0}$, (iv) $\hat{Q}(\bm{\theta},n)$ converges uniformly in probability to $Q(\bm{\theta})$.

On one hand, Assumptions A2 ensures that $\bm{\Theta}$ is compact. On the other hand, Assumption D2 guarantees that ${Q}(\bm{\theta})$ is continuous and uniquely minimized at ${\bm{\theta}}_{0}$ since $\pi({\bm{\theta}})={\bm{\theta}}+\mathbf{a}({\bm{\theta}})$ is required to be continuous and injective. What remains to be shown is that $\hat{Q}(\bm{\theta},n)$ converges uniformly in probability to $Q(\bm{\theta})$, which is equivalent to show that: $\forall\varepsilon>0$ and $\forall\delta>0$, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\geq n^{\ast}$

[TABLE]

Fix $\varepsilon>0$ and $\delta>0$. Using the above definitions, we have that

[TABLE]

where

[TABLE]

We now consider each term of (E.2) separately. For simplicity, we define

[TABLE]

and, considering the first term of (E.2), we have

[TABLE]

by Assumption C. Similarly, we have

[TABLE]

by Assumption D2. Therefore, we obtain

[TABLE]

By Assumptions C and D2, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{\ast}$ we have

[TABLE]

Therefore, the four condition of Theorem 2.1 of Newey and McFadden (1994) are verified implying the result. ∎

Proposition 2 implies that when $p\to c<\infty$ as $n\to\infty$ we can simply write:

[TABLE]

However, in the case where $p\to\infty$ as $n\to\infty$, the statement is weaker in the sense that we cannot conclude that the limit exist but only that $\Pr(||\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}||_{2}>\varepsilon)$ is arbitrarily small.

Appendix F Phase Transition Unbiasedness of the IB-Estimator

In this appendix, we study the finite sample bias of the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ and show that it is a PT-unbiased estimator. For this purpose, a slightly modified assumption framework for Theorem 2 is considered, where a modified version of Assumption D∗ is introduced, while Assumption A and B remains unchanged. Assumption C is not needed here since we are considering the bias of an estimator. For completness, we first restate Assumption A using our notation.

Assumption A:

Let $\bm{\Theta}$ be a convex and compact subset of ${\rm I\!R}^{p}$ such that ${\bm{\theta}}_{0}\in\operatorname{Int}(\bm{\Theta})$ and

[TABLE]

Next, we introduce Assumption D3.

Assumption D3:

Consider the functions $\mathbf{d}(\bm{\theta},n)$, $\mathbf{a}({\bm{\theta}})$, $\mathbf{b}(\bm{\theta},n)$ and $\mathbf{r}(\bm{\theta},n)$ defined in (5) and (6) and let $\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}(\bm{\theta},n)\in{\rm I\!R}^{p\times p}$ be the Jacobian matrix of $\mathbf{r}(\bm{\theta},n)$ in ${\bm{\theta}}$. These functions are such that for all ${\bm{\theta}}\in\bm{\Theta}$:

There exists a matrix $\mathbf{M}\in{\rm I\!R}^{p\times p}$ and a vector $\mathbf{s}\in{\rm I\!R}^{p}$ such that $\mathbf{a}({\bm{\theta}})=\mathbf{M}{\bm{\theta}}+\mathbf{s}$. 2. 2.

There exists a sample size $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ the matrix $(\mathbf{M}+\mathbf{L}(n)+\mathbf{I})^{-1}$ exists. 3. 3.

There exists a real $\gamma>0$ such that for all $j=1,\dots,p$ we have

[TABLE] 4. 4.

$\mathbf{R}(\bm{\theta},n)$* exists and is continuous in ${\bm{\theta}}\in\bm{\Theta}$ for all $n\in\mathbb{N}^{\ast}$ with $n\geq n^{\ast}$.*

Compared to Assumption D∗, Assumption D3 is weaker. First, the matrix $\mathbf{M}$ may be such that $\lVert\mathbf{M}\rVert_{F}\geq 1$, thereby relaxing the contraction mapping hypothesis on $\mathbf{a}({\bm{\theta}})$ given in Assumptions D, D∗ and D1. Second, conditions 3 and 5 of Assumption D∗ are removed.

Clearly, neither Assumption D2 nor Assumption D3 imply each other. The first part of Assumption D3 imposes a restrictive form for the bias function $\mathbf{a}({\bm{\theta}})$. If $\hat{\bm{\pi}}(\bm{\theta}_{0},n,\bm{\omega})$ is a consistent estimator of $\bm{\theta}_{0}$, then this restriction is automatically satisfied. Under these new conditions, the PT-unbiasedness of $\hat{\bm{\theta}}_{(n,H)}$ is established in Proposition 3 below.

Proposition 3:

Under Assumptions A, B and D3, $\hat{\bm{\theta}}_{(n,H)}$ is a PT-unbiased estimator, i.e. there exists a sample size $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ and for all $H\in\mathbb{N}^{\ast}$, we have \Big{\lVert}\mathbb{E}\left[\hat{\bm{\theta}}_{(n,H)}\right]-\bm{\theta}_{0}\Big{\rVert}_{2}=0.

Before giving the proof we introduce the following useful notation. Let $\mathbf{f}\,:\bm{\Theta}\to{\rm I\!R}^{p}$ and $\mathbf{F}({\bm{\theta}})~{}\equiv~{}\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{f}\left(\bm{\theta}\right)\in{\rm I\!R}^{p\times p}$ be its Jacobian matrix. Then, we define \mathbf{F}\big{(}{\bm{\theta}}^{(\mathbf{f})}\big{)} such that when using the multivariate mean value theorem between $\hat{\bm{\theta}}_{(n,H)}$ and ${{\bm{\theta}}}_{0}$ we obtain

[TABLE]

Therefore, ${\bm{\theta}}^{(\mathbf{f})}$ corresponds to a set of $p$ vectors lying in the segment $(1-\lambda)\hat{\bm{\theta}}_{(n,H)}+\lambda{{\bm{\theta}}}_{0}$ for $\lambda\in[0,1]$ (with respect to the function $\mathbf{f}$). Keeping the latter notation in mind (for a generic function $\mathbf{f}$), we provide the proof of Proposition 3.

Proof: From (B.1) we have that $\hat{\bm{\theta}}_{(n,H)}$ is such that for $(n,H)\in\mathbb{N}^{\ast}\times\mathbb{N}^{\ast}$

[TABLE]

Using Assumptions B and D3, we expand each side of the above equation as follows

[TABLE]

Therefore, we obtain

[TABLE]

By hypothesis, $n$ is such that the inverse of $\mathbf{B}\equiv\mathbf{I}+\mathbf{L}(n)+\mathbf{M}$ exists and we obtain

[TABLE]

By Assumptions A and D3, $\mathbf{r}(\hat{\bm{\theta}}_{(n,H)},n)$ is a bounded random variable on a compact set. Moreover, since $\mathbf{r}(\bm{\theta},n)=\mathcal{O}(n^{-\gamma})$ elementwise by Assumption D3, we have

[TABLE]

elementwise. Consequently, we deduce from (F.3) that

[TABLE]

elementwise.

The idea now is to re-evaluate $\mathbb{E}\left[\mathbf{r}(\hat{\bm{\theta}}_{(n,H)},n)-\mathbf{r}(\bm{\theta}_{0},n)\right]$ using the mean value theorem. This will allow us to make an induction that will show that for all $\delta\in\mathbb{N}$

[TABLE]

We will then use Lemma 1 to conclude the proof.

Applying the mean value theorem for vector-valued functions to $\mathbf{r}(\hat{\bm{\theta}}_{(n,H)},n)-\mathbf{r}(\bm{\theta}_{0},n)$ we have

[TABLE]

where the matrix $\mathbf{R}\left(\bm{\theta}^{(\mathbf{r})},n\right)$ is defined in (F.1). By Assumptions A and D3, $\mathbf{R}\left(\bm{\theta}^{(\mathbf{r})},n\right)$ is also a bounded random variable. Moreover, we have $\mathbb{E}\left[\mathbf{R}\left(\bm{\theta}^{(\mathbf{r})},n\right)\right]=\mathcal{O}(n^{-\gamma})$ elementwise since by Assumption D3 $\mathbf{r}(\bm{\theta},n)=\mathcal{O}(n^{-\gamma})$ elementwise, and

[TABLE]

for all $\bm{\theta}\in\bm{\Theta}$. For simplicity, we denote

[TABLE]

which implies that $\mathbf{\Delta}^{(\mathbf{r})}=\mathbf{R}\mathbf{\Delta}$. Moreover, a consequence of (F.5) is that

[TABLE]

for any $j=1,\dots,p$ and hence

[TABLE]

Now, for all $l=1,\dots,p$ we have

[TABLE]

Using Cauchy-Schwarz inequality, we have

[TABLE]

Since $\mathbf{R}_{l,m}$ and $\mathbf{\Delta}_{m}$ are bounded random variables on a compact set, $\mathbb{E}\left[\mathbf{R}_{l,m}\right]=\mathcal{O}\left(n^{-\gamma}\right)$ and $\mathbb{E}\left[\mathbf{\Delta}_{m}\right]=\mathcal{O}\left(pn^{-\gamma}\right)$ we have

[TABLE]

Therefore, we obtain

[TABLE]

and hence, using (F.6) we deduce that

[TABLE]

Considering this, and since $\mathbb{E}\left[\mathbf{\Delta}\right]=-\mathbf{B}^{-1}\mathbb{E}\left[\mathbf{\Delta}^{(\mathbf{r})}\right]$, we have, for all $j=1,\dots,p$

[TABLE]

and consequently,

[TABLE]

Since $\mathbb{E}\left[\mathbf{\Delta}\right]=-\mathbf{B}^{-1}\mathbb{E}\left[\mathbf{\Delta}^{(\mathbf{r})}\right]$ and $\mathbb{E}\left[\mathbf{\Delta}^{(\mathbf{r})}\right]=\mathbb{E}\left[\mathbf{R}\mathbf{\Delta}\right]$, one can repeat the same computations and deduce by induction that for all $\delta\in\mathbb{N}$

[TABLE]

Using equation (F.7) we see that

[TABLE]

which ends the proof by Lemma 1. ∎

Proposition 3 provides a strategy for proving the consistency of the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$. However, this requires slight modifications of the assumptions considered in Proposition 3. More in detail, we impose stronger requirements on $\alpha$ and $\gamma$ in Assumptions C and D3.

Assumption C1:

Preserving the same definitions and requirements of Assumption C, we additionally require that the constant $\alpha$ is such that

[TABLE]

Assumption D4:

Preserving the same definitions and requirements of Assumption D3, we additionally require that the constant $\gamma$ is such that

[TABLE]

Clearly, Assumption C is weaker than Assumption C1. Neither Assumption D4 nor Assumption D imply each other. Indeed, in the situations where $\alpha=\nicefrac{{1}}{{2}}$ and $\gamma=2$, the estimator $\hat{\bm{\theta}}_{(n,H)}$ is consistent for $\bm{\theta}_{0}$ under the conditions of Proposition 3 and if

[TABLE]

which is clearly a stronger requirement than (E.1) based on Proposition 2 (using Assumption D2 which is weaker than Assumption D). However, some requirements of Assumption D are stronger than the ones in Assumption D4. An important difference is that the function $\mathbf{a}(\bm{\theta})$ doesn’t need to be a contraction map in Assumption D4. In addition, most of the requirements on the matrix $\mathbf{L}(n)$ and the vector $\mathbf{c}(n)$ are removed. Using these new assumptions, we present the following corollary to Proposition 3.

Corollary 1:

Under Assumptions A, B, C1, and D4, $\hat{\bm{\theta}}_{(n,H)}$ is a consistent estimator of $\bm{\theta}_{0}$ for all $H\in\mathbb{N}^{\ast}$ in the sense of Proposition 2.

Proof: Fix $\varepsilon>0$ and $\delta>0$. We need to show that there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{\ast}$

[TABLE]

From Chebyshev’s inequality we have that

[TABLE]

Therefore, we only need to show that there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\geq n^{\ast}$

[TABLE]

Using the same notation as in the proof of Proposition 3, we have from (F.2) and for $n$ sufficiently large that

[TABLE]

where $\mathbf{v}$ is zero mean random vector of order $\mathcal{O}_{\rm p}\left(n^{-\alpha}\right)$ elementewise by Assumption C. Thus, by Assumption D4 we have that for all $j=1,\dots,p$

[TABLE]

Therefore we have,

[TABLE]

and thus since $\bm{\Theta}$ is compact by Assumption A, we have

[TABLE]

Using Assumptions C1 and D4 the last equality implies that there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{\ast}$, we obtain

[TABLE]

∎

Appendix G Asymptotic Distribution of the IB-Estimator

In order to study the asymptotic normality of the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ we modify some of the assumptions in order to incorporate additional requirements on the auxiliary estimator.

Before introducing the next assumption we define various quantities. First, we let

[TABLE]

where $\bm{\Sigma}({\bm{\theta}},n)$ is nonsingular. Then, using this definition, we introduce the following quantity

[TABLE]

where $\mathbf{u}\in{\rm I\!R}^{p}$ is such that $||\mathbf{u}||_{2}=1$. Clearly, from the definition of $\mathbf{v}(\bm{\theta},n,\bm{\omega})$ we have that $\mathbb{E}[Y({\bm{\theta}},n,\bm{\omega},\mathbf{u})]=0$. Therefore, without loss of generality we can always decompose $Y({\bm{\theta}},n,\bm{\omega},\mathbf{u})$ as follows:

[TABLE]

where $W({\bm{\theta}},n,\bm{\omega},\mathbf{u})$ and $Z({\bm{\theta}},n,\bm{\omega},\mathbf{u})$ are zero mean random variables and \big{(}\delta_{n}\big{)}\in\mathfrak{D} where

[TABLE]

In Assumption C∗ below, we restrict the behaviour of $\mathbf{v}(\bm{\theta},n,\bm{\omega})$ and in particular we require that it satisfies a specific Gaussian approximation based on the decomposition considered in (G.1).

Assumption C∗:

Preserving the requirements of Assumption C we add the following conditions on $\mathbf{v}(\bm{\theta},n,\bm{\omega})$. There exists a sample size $n^{\ast}\in\mathbb{N}^{*}$ such that:

For all $\bm{\theta}\in\bm{\Theta}$ and all $n\in\mathbb{N}^{*}$ with $n\geq n^{\ast}$, the matrix $\bm{\Sigma}({\bm{\theta}},n)$ exists, is nonsingular and is such that $\bm{\Sigma}_{j,l}({\bm{\theta}},n)=\mathcal{O}(1)$ for all $j,l=1,\dots,p$. 2. 2.

For all $\bm{\omega}\in\bm{\Omega}$ and all $n\in\mathbb{N}^{*}$ with $n\geq n^{\ast}$, the Jacobian matrix $\mathbf{V}(\bm{\theta},n,\bm{\omega})\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{v}\left(\bm{\theta},n,\bm{\omega}\right)$ exists and is continuous in ${\bm{\theta}}\in\bm{\Theta}$. 3. 3.

Considering the decomposition in (G.1), then for all $\mathbf{u}\in{\rm I\!R}^{p}$ such that $||\mathbf{u}||_{2}=1$, there exist sequences \big{(}\delta_{n}\big{)}\in\mathfrak{D}, as well as a random variable $W({\bm{\theta}},n,\bm{\omega},\mathbf{u})=\mathcal{O}_{\rm p}(1)$, such that $Z({\bm{\theta}},n,\bm{\omega},\mathbf{u})$ is a standard normal random variable.

The first requirement of Assumption C∗ is quite mild and commonly assumed as it simply requires that $\bm{\Sigma}({\bm{\theta}},n)$ is a suitable covariance matrix. In addition, it implies that $\mathbf{v}(\bm{\theta},n,\bm{\omega})=\mathcal{O}_{\rm p}\left(n^{-\nicefrac{{1}}{{2}}}\right)$, elementwise. Similarly to the discussion following Assumption D∗, our second requirement ensures that $\mathbf{V}_{j,l}(\hat{\bm{\theta}}_{(n,H)},n,\bm{\omega})$ is a bounded random variable for all $j,l=1,\dots,p$ and all $(n,\bm{\omega})\in\mathbb{N}^{*}\times\bm{\Omega}$. The third condition of Assumption C∗ describes how “close” the distribution of $\mathbf{v}(\bm{\theta},n,\bm{\omega})$ is to a multivariate normal distribution. Such an assumption is quite strong and may not always be satisfied. In the case where $p\to c<\infty$ as $n\to\infty$, our requirement on the distribution of $\mathbf{v}(\bm{\theta},n,\bm{\omega})$ can simply be expressed as

[TABLE]

where $\bm{\Sigma}({\bm{\theta}})\equiv\displaystyle\lim_{n\to\infty}\bm{\Sigma}({\bm{\theta}},n)$. Finally, Assumption C∗ clearly implies Assumption C but not necessarily Assumption C1.

In Assumption D5 below, we impose an additional requirement on the bias of the auxiliary estimator.

Assumption D5:

There exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{\ast}$ the Jacobian matrices $\mathbf{A}(\bm{\theta})\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{a}\left(\bm{\theta}\right)$ and $\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}\left(\bm{\theta},n\right)$ exist and are continuous in ${\bm{\theta}}\in\bm{\Theta}$. Moreover, we require that for all $\mathbf{u}\in{\rm I\!R}^{p}$ with $||\mathbf{u}||_{2}=1$, there exists a sequence $(\delta_{n})\in\mathfrak{D}$ such that

[TABLE]

where

[TABLE]

[TABLE]

while $\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})$ and $\mathbf{R}({\bm{\theta}}^{(\mathbf{r})})$ are defined using notation in (F.1).

Assumption D5 allows us to quantify how “far” the matrices $\mathbf{B}(\bm{\theta}_{0},n)$ and \mathbf{B}\big{(}\hat{\bm{\theta}}_{(n,H)},n\big{)} are from each other. In the case where $p\to c<\infty$ as $n\to\infty$, the consistency of $\hat{\bm{\theta}}_{(n,H)}$ and the continuity of $\mathbf{B}(\bm{\theta})\equiv\displaystyle\lim_{n\to\infty}\mathbf{B}(\bm{\theta},n)$ in $\bm{\theta}_{0}$ would be sufficient so that Assumption D5 would not be needed to establish the asymptotic distribution of $\hat{\bm{\theta}}_{(n,H)}$. However, when $p\to\infty$ as $n\to\infty$ the requirement in (G.2) can be strong and difficult to verify for a specific model and auxiliary estimator. Having stated this, the following proposition defines the distribution of the IB-estimator.

Proposition 4:

Under Assumptions A, B, C∗ and D5, for all $\mathbf{u}\in{\rm I\!R}^{p}$ such that $||\mathbf{u}||_{2}=1$, there exist a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ and a sequence $\left(\delta_{n}\right)\in\mathfrak{D}$ such that for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have

[TABLE]

where $Z\sim\mathcal{N}(0,1)$.

Proof: Without loss of generality we can assume that the sample sizes $n^{*}\in\mathbb{N}^{*}$ of Assumption C∗ and D5 are the same. Let $n\in\mathbb{N}^{*}$ be such that $n\geq n^{*}$. By the definition of the IB-estimator we have

[TABLE]

Since the auxiliary estimator may be expressed, for any $({\bm{\theta}},n,\bm{\omega})\in\bm{\Theta}\times\mathbb{N}^{*}\times\bm{\Omega}$, as

[TABLE]

we deduce from (G.3), by using the mean value theorem on $\mathbf{a}({\bm{\theta}})$ and $\mathbf{r}({\bm{\theta}},n)$, that

[TABLE]

where $\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\equiv\mathbf{I}+\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})+\mathbf{L}(n)+\mathbf{R}\left({\bm{\theta}}^{(\mathbf{r})},n\right)$ and $\bm{\Delta}\equiv\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}$. Applying the mean value theorem on $\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,{\bm{\omega}_{h}}\right)$ and using the notation defined in (F.1), we obtain for all $h=1,\dots,H$

[TABLE]

Hence, let $\mathbf{u}\in{\rm I\!R}^{p}$ be such that $||\mathbf{u}||_{2}=1$, then by Assumption C∗ we have,

[TABLE]

and for all $h=1,\dots,H$

[TABLE]

where $(\delta_{n}^{(0)}),(\delta_{n}^{(1)}),\dots,(\delta_{n}^{(H)})\in\mathfrak{D}$ and

[TABLE]

Without loss of generality, we can suppose that $(\delta_{n})\equiv(\delta_{n}^{(0)})=(\delta_{n}^{(1)})=\dots=(\delta_{n}^{(H)})$ as one may simply modify $W(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})$ and $W(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})$ for all $h=1,\dots,H$. For simplicity, we write $\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\equiv\sqrt{n}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}$, $Z_{0}\equiv Z(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})$ and $Z_{h}\equiv Z(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})$. Multiplying the right hand side of (G.5) by $\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}$, we have

[TABLE]

The first three terms of (G.6) are considered separately in the following computation

[TABLE]

where $Z\sim\mathcal{N}\left(0,1\right)$. We now consider the last term of (G.6). Since $\bm{\Sigma}({\bm{\theta}},n)^{-\nicefrac{{1}}{{2}}}=\mathcal{O}(1)$ and $\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)=\mathcal{O}_{p}\left(n^{-\nicefrac{{1}}{{2}}}\right)$ elementwise by Assumption C∗ and $||\mathbf{u}||_{2}=1$ we have for all $h=1,\dots,H$

[TABLE]

where $\mathbf{1}_{p}\in{\rm I\!R}^{p}$ is a vector of ones. Since Assumptions A, B, C∗ and D5 imply that $\hat{\bm{\theta}}_{(n,H)}$ is consistent estimator of $\bm{\theta}_{0}$ (see Figure 4), there exists a sequence $(\delta^{\prime}_{n})\in\mathfrak{D}$ such that $\bm{\Delta}=\mathcal{O}_{\rm p}(\delta^{\prime}_{n})$ elementwise for all $n\in\mathbb{N}^{*}$. Without loss of generality, we may assume that $(\delta^{\prime}_{n})=(\delta_{n})$. Therefore, we have

[TABLE]

which leads to the computation for the last term of (G.6),

[TABLE]

Combining (G.5), (G.6) and (G.8), we obtain

[TABLE]

The proof is completed using (G.2) of Assumption D5. ∎

Remark D:

From the proof of Proposition 4, we have that

[TABLE]

without using (G.2) in Assumption D5. As a result, Assumption D5 is not strictly necessary for the Gaussian approximation of the distribution of $\hat{\bm{\theta}}_{(n,H)}$. However, while this assumption is strong, it is quite convenient to deliver a reasonable approximation of $\mathbf{B}(\bm{\theta}^{(\mathbf{a},\mathbf{r})},n)$ in (G.9) when $p$ is allowed to diverge, since its applicability is quite limited in practice.

Appendix H Proofs of the Main Results

The results and assumptions presented so far can now be combined to deliver the three theorems of this article. While Theorem 1 and 2 were stated in the main text, Theorem 3, which is stated below, presents the asymptotic distribution of the IB-estimator. Theorem 3 is based on Assumption D*∗∗*, which is the union of Assumption D∗ and Assumption D5. Assumption D*∗∗* is provided below for clarity of exposition.

Assumption D*∗∗*:

The bias function $\mathbf{d}\left(\bm{\theta},n\right)$ is such that:

The asymptotic bias function $\mathbf{a}(\bm{\theta})$ can be written as

[TABLE]

where $\mathbf{M}\in{\rm I\!R}^{p\times p}$ with $||\mathbf{M}||_{F}<1$ and $\mathbf{s}\in{\rm I\!R}^{p}$. 2. 2.

There exists a sample size $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ the matrix $(\mathbf{M}+\mathbf{L}(n)+\mathbf{I})^{-1}$ exists. 3. 3.

There exist real $\beta,\gamma>0$ such that for all $\bm{\theta}\in\bm{\Theta}$ and any $j,l=1,\dots,p$, we have

[TABLE] 4. 4.

Defining $c_{n}\equiv\displaystyle{\max_{j=1,\dots,p}}\mathbf{c}_{j}(n)$ for all $n\in\mathbb{N}^{*}$, we require that the sequence $\left\{c_{n}\right\}_{n\in\mathbb{N}^{*}}$ is such that

[TABLE] 5. 5.

There exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ the Jacobian matrices $\mathbf{A}(\bm{\theta})\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{a}\left(\bm{\theta}\right)$ and $\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}\left(\bm{\theta},n\right)$ exist and are continuous in ${\bm{\theta}}\in\bm{\Theta}$. Moreover, we require that for all $\mathbf{u}\in{\rm I\!R}^{p}$ with $||\mathbf{u}||_{2}=1$, there exists a sequence $(\delta_{n})\in\mathfrak{D}$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

Using this new assumption, we now restate Theorem 1 and 2, and state Theorem 3, followed by their proof.

Theorem 1:

Under Assumptions A, B, C and D, for all $H\in\mathbb{N}^{\ast}$,

There exist a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}$ with $n\geq n^{*}$, $\left\{\hat{\bm{\theta}}_{(n,H)}\right\}=\widehat{\bm{\Theta}}_{(n,H)}$, i.e. the set $\widehat{\bm{\Theta}}_{(n,H)}$ is a singleton. 2. 2.

There exist a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}$ with $n\geq n^{*}$, the sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ has the following limit

[TABLE]

Moreover, there exists a real $\epsilon\in(0,\,1)$ such that for any $k\in\mathbb{N}^{\ast}$

[TABLE] 3. 3.

$\hat{\bm{\theta}}_{(n,H)}$* is a consistent estimator of $\bm{\theta}_{0}$, in that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have:*

[TABLE]

Theorem 2:

Under Assumptions A, B, C and D∗, and for all $H\in\mathbb{N}^{\ast}$, the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ is consistent and PT-unibiased. Therefore, the IB-estimator satisfies:

$\hat{\bm{\theta}}_{(n,H)}$* is a consistent estimator of $\bm{\theta}_{0}$, in that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have:*

[TABLE] 2. 2.

$\hat{\bm{\theta}}_{(n,H)}$* is PT-unbiased, that is there exists a $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$, we have \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0.*

Theorem 3:

Under Assumptions A, B, C∗ and D*∗∗*, and for all $H\in\mathbb{N}^{\ast}$, the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ satisfies:

$\hat{\bm{\theta}}_{(n,H)}$* is a consistent estimator of $\bm{\theta}_{0}$, in that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such for all $n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have:*

[TABLE] 2. 2.

$\hat{\bm{\theta}}_{(n,H)}$* is PT-unbiased, that is there exists a $n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$, we have \big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0.* 3. 3.

For all $\mathbf{u}\in{\rm I\!R}^{p}$ such that $||\mathbf{u}||_{2}=1$, there exist a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ and a sequence $\left(\delta_{n}\right)\in\mathfrak{D}$ such that for all for all $n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ we have

[TABLE]

where $Z\sim\mathcal{N}(0,1)$.

Proof: The figure below delivers a proof of Theorems 1, 2 and 3.

∎

Remark E:

Other results can be derived from Figure 4. For example, one can state a theorem ensuring that the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ is unbiased and asymptotically normal (and therefore consistent) without the guarantee that the IB sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ is convergent. This can be of interest if another method is used to compute the IB-estimator. Another case is when one is only interested in the convergence of the IB sequence $\left\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\}_{k\in\mathbb{N}}$ and the consistency of the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$. Assumption A is stronger than necessary for such a result. The reason we use this assumption to prove Theorem 1 is guided by the rationale that we want a restricted number of general assumptions for the main Theorems 1, 2 and 3.