Efficient Bias Correction for Cross-section and Panel Data

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

84,517 characters

Efficient Bias Correction for Cross-section and Panel Data



\title{Efficient Bias Correction for Cross-section and Panel Data\footnote{T. Rothenberg provided helpful comments.}}
\author{Jinyong Hahn\thanks{Department of Economics, UCLA, Los Angeles, CA 90095-1477. Email: [email removed]}\\UCLA
\and David W. Hughes\thanks{Boston College, Department of Economics, Chestnut Hill, MA 02467. Email: [email removed].}\\Boston College
\and Guido Kuersteiner\thanks{University of Maryland, Department of Economics, College Park, MD 20742-7211. Email:
[email removed]. Financial support from NSF grants SES-0095132 and
SES-0523186 is gratefully acknowledged. }\\University of Maryland
\and Whitney K. Newey\thanks{MIT, Department of Economics, Cambridge, MA 02139. Email: [email removed]. Financial support from NSF Grant 1757140 is gratefully acknowledged.}\\MIT and NBER}
\maketitle

\begin{abstract}
Bias correction can often improve the finite sample performance of estimators. We show that the choice of bias correction method has no effect on the higher-order variance of semiparametrically efficient parametric estimators, so long as the estimate of the bias is asymptotically linear. It is also shown that bootstrap, jackknife, and analytical bias estimates are asymptotically linear for  estimators with higher-order expansions of a standard form. In particular, we find that for a variety of estimators the straightforward bootstrap bias correction gives the same higher-order variance as more complicated analytical or jackknife bias corrections. In contrast, bias corrections that do not estimate the bias at the parametric rate, such as the split-sample jackknife, result in larger higher-order variances in the i.i.d.\ setting we focus on. For both a cross-sectional MLE and a panel model with individual fixed effects, we show that the split-sample jackknife has a higher-order variance term that is twice as large as that of the `leave-one-out' jackknife.

\end{abstract}

\newpage


\section{Introduction}

Asymptotic bias corrections can be useful for centering estimators nearer to the
truth. One approach is to use analytical corrections such as the standard
textbook expansion for functions of sample means and the more complicated
formulas required for other estimators. Alternatively, we may use jackknife and
bootstrap bias corrections. To help choose a bias correction method it would
be useful to know which, if any, is preferable on asymptotic efficiency
grounds. Although the bias correction does not affect the first-order asymptotic
variance, it can affect the higher-order variance.  We show that the method of bias correction does not affect the higher-order variance of any parametric estimator that is efficient in a semiparametric model,
as long as the bias estimator is asymptotically linear. Thus, one can choose a bias correction
for an efficient estimator based on computational convenience, or some other
criteria, without affecting its higher-order efficiency. We give a formal
expansion showing this property for a parametric estimator in a general semiparametric model, i.e. a model with a parametric component in which some other components are left unspecified. We also prove
that the bootstrap, jackknife and analytical bias estimates are asymptotically linear, when the estimator of the parameters of interest has a standard form of stochastic expansion (which is known to exist for a large class of models). Derivations in the case of the MLE, show that the jackknife, bootstrap, and one type of analytical bias correction deliver estimators that have identical stochastic expansions to third order, and so have an even stronger equivalence property.

There are many implications of this higher-order efficiency result. One is that bias corrections that are not asymptotically linear may not have the same higher-order variance as those that are.
For example, split-sample jackknife bias corrections are not asymptotically linear in cross-section or panel data and have a larger higher-order variance than other bias-corrected estimators.
We find that the higher-order variance term is twice the size of that for the leave-one-out jackknife bias correction. On the other hand, the split-sample jackknife is useful in time series or panel data when the observations are not independent over time because the leave-one-out panel jackknife does not work in this case.

Another implication of the result is that it allows researchers to choose the bias correction method that is computationally convenient.
For example, Newey and Smith (2004) showed that the empirical likelihood estimator is higher-order efficient in moment condition models when certain analytical bias corrections are used.
Asymptotic linearity of the bootstrap means that calculation of the bias formula can be avoided by using the bootstrap bias correction instead.
As another example, Cattaneo, Jansson, and Ma (2019)
showed that a jackknife bias correction can be important when a first step regression with many regressors is plugged into a second step regression. Although our current results do not include asymptotics in which the number of regressors increases with the sample size, we conjecture that the bootstrap bias correction has a similar properties as the jackknife bias correction when the number of regressors increases slowly enough.

The higher-order variance concept that we consider is the $O(
n^{-1})$ variance of a third-order stochastic expansion of the
estimator. Its use for comparison of estimators was pioneered by Nagar (1959).
As shown in Pfanzagl and Wefelmeyer (1978) and Ghosh, Sinea, and Wieand (1980),
and discussed in Rothenberg (1984), under appropriate regularity conditions
rankings based on this higher-order variance correspond to rankings based on
the variance of an Edgeworth approximation. Thus, the bias and variance of leading terms in a stochastic expansion are also the leading
terms of an expansion of the bias and variance of an approximating
distribution. Furthermore, as noted by Rothenberg (1984), Akahira and
Takeuchi (1981) have shown that all well-behaved asymptotically efficient
estimators of a parametric likelihood model necessarily have the same skewness
and kurtosis to an order $n^{-1}$ approximation and that to compare the
dispersion of second-order approximations to the distribution of efficient estimators it suffices to compare their higher-order variances.
This motivates our focus here on the higher-order variance of bias-corrected estimators, which serves to quantify their higher-order
efficiency.
In having this focus we follow much of the more recent literature, such as Rilstone, Srivastava, and Ullah (1996) and Newey and Smith (2004).

We also derive asymptotic higher-order variance expressions for panel data
models with unobserved, individual fixed
effects as well as for cross-section models. These estimators are known to suffer from asymptotic bias under asymptotic sequences
in which $n$ and $T$ grow at the same rate. Two common methods of bias correction in this setting are the
`leave-one-out' panel jackknife (Hahn and Newey, 2004), and the split-sample jackknife (Dhaene and Jochmans, 2015).
Both deliver estimates that are asymptotically normal and centered at the truth when
$n$ and $T$ grow at the same rate, with equal first-order (asymptotic) variances.
Our analysis makes it possible to compare the two bias corrections in terms of their higher-order variance. We find that with i.i.d.\ data the split-sample correction has a higher-order variance that is twice the size of the `leave-one-out' jackknife. Although we focus on the maximum likelihood setting for panel data, the results are applicable to a broader set of moment condition estimators under suitable assumptions on the moment functions.
Numerical comparisons in recent papers confirm that the difference in higher-order variance can be meaningful in practice, with the split-sample jackknife having larger dispersion and lower coverage than analytical or leave-one-out jackknife corrections in a variety of settings; see for example,
Alexander and Breunig (2016), Fernandez-Val and Weidner (2018), and Czarnowske and Stammann (2019). Our own simulations of a panel probit model with one
common parameter and individual fixed effects also support this result.
This comparison is also true of estimates of a marginal effect parameter.


\subsection{Related literature}

Higher-order efficiency of the MLE was analyzed by Pfanzagl and Wefelmeyer (1978) in terms of risk functions or Akahira and Takeuchi (1981) in terms of concentration probabilities. Ghosh (1994), and Taniguchi and Kakizawa (2000) for time series models, contain surveys of this literature. A common theme is that a higher-order bias-corrected version of the MLE is higher-order efficient in the case of higher-order squared risk. Similar results obtain for median bias-corrected MLE's in the case of concentration probabilities. The bias correction in this literature is typically of a known parametric form that only depends on the estimated parameters. A plug-in estimator is then a regular estimator for the bias term. Amari (1982) obtains similar results for curved exponential families using differential geometry that characterizes the MLE in terms of tangent spaces. Akahira (1983, 1989) shows that when instead of using parametric bias correction one relies on the jackknife the same higher-order efficiency results remain true for the jackknifed MLE. We add to this literature by analyzing the effects of bias correction for first-order, but not necessarily higher-order, semiparametrically efficient estimators.  We demonstrate that any bias-corrected semiparametrically efficient estimator using an asymptotically linear bias estimator has a higher-order variance that does not depend on the nature of the bias estimate. This result applies to cases where the higher-order bias may not be known in closed form. We show that efficient bias correction may be based on sample averages, bootstrap or jackknife methods as long as the bias estimates are asymptotically linear.

The use of a jackknife bias estimator goes back to Quenouille (1949, 1956) and
Tukey (1958). Bootstrap bias estimation was discussed by Parr (1983), Shao
(1988a,b), Hall (1992), and Horowitz (1998) in the context of nonlinear
transformations of OLS estimators of linear models and nonlinear functions of
the mean. Akahira (1983) showed that the jackknife bias-corrected maximum
likelihood estimator is higher-order efficient. Our work extends the
literature on the bootstrap and jackknife bias-corrected estimators by
analyzing any semiparametrically efficient parametric estimator rather than nonlinear transformations
of linear estimators as in Shao (1988a,b).

In the panel data literature, methods to control for unobserved heterogeneity are well established
(early literature includes Rasch (1960, 1961) and Andersen (1970)); see for example Chamberlain (1984), Arellano and Honor\'{e} (2001), and Arellano and Hahn (2010) for reviews. Because of the incidental parameters problem, the
best that can be achieved in a fixed-$T$ setting is partial identification in general;
this is especially true for policy relevant parameters such as average marginal
effects (see Chernozhukov, Newey, Hahn and Fern\'{a}ndez-Val (2013)). Even
under sequences in which $T$ grows at the same rate as $n$, fixed effects estimators may be
asymptotically biased, as discussed in Hahn and Kuersteiner (2002) and Hahn
and Newey (2004). Given the typical size of panel data sets, in which $n$ is
much larger than $T$, it is desirable to find estimators that have biases of
order $O(T^{-2})$ or smaller, rather than the typical $O(T^{-1})$ of fixed effects estimators.

For a static model, Hahn and Newey (2004) show that a `leave-one-out'
jackknife estimator is asymptotically normal and centered at the truth when
$n$ and $T$ grow at the same rate. Other styles of jackknife bias correction
are also possible. For example, Dhaene and Jochmans (2015) suggest a
split-sample bias correction that, in its simplest form, is constructed by
splitting the sample into two half-panels of length $T/2$. It should be understood that the split-sample jackknife provides valid bias corrections with autocorrelated data where the `leave-one out' jackknife does not. Thus the split-sample jackknife is preferred to the leave-one-out jackknife with autocorrelated data. Our results show that the leave-one-out jackknife is preferred to the split-sample jackknife in i.i.d.\ data in the sense that the higher-order variance is smaller and small sample performance is better.

The remainder of the paper is set out as follows. In Section \ref{sec:semipar_result} we discuss the higher-order bias and variance of parametric estimators, and provide our main results on the higher-order efficiency of bias corrections for semiparametric efficient estimators and asymptotic linearity of analytical, bootstrap and jackknife bias corrections. In Section \ref{sec:mle_crosssec}, we provide expressions for the higher-order expansions and variances of various bias-corrections in a cross-sectional MLE. We extend the results to panel settings
in Section \ref{section-panel-model}, by deriving the higher-order variances of the leave-one-out and split-sample jackknifes in a model with individual fixed effects.
Section \ref{section-MC} provides Monte Carlo evidence to support the theory in the panel setting. Section \ref{section-summary} concludes.









\section{Higher-order efficiency of bias-corrected estimators} \label{sec:semipar_result}

\subsection{Higher-order bias and variance}
We begin with a discussion of higher-order expansions for parametric estimators, and make precise our definition of higher-order variance, following closely the exposition in Rothenberg (1984). We focus on semiparametric models in which the
data $\left\{  Z_{i}\right\}  _{i=1}^{n}$ are i.i.d.\ with some
distribution $F_{0}$ contained in a class of distributions that is the model.\footnote{By `semiparametric model', we refer to a model which has a parametric component, but leaves the functional form of some other components unspecified. See Newey (1990) for further discussion.} Let $\widehat{\theta}$ denote an estimator of a parameter
$\theta_{0}.$ When $\widehat{\theta}$ is an estimator based on parametric moment conditions,
such as GMM, and the moment conditions are smooth enough in the parameter
$\theta$ and possibly other parameters, there will be a stochastic expansion of
the form
\begin{align} \label{eq:expan}
\sqrt{n}\left(  \widehat{\theta}-\theta_{0}\right)   &  =A_{1}+\frac{A_{2}}
{\sqrt{n}}+\frac{A_{3}}{n}+o_{p}\left( n^{-1} \right)  ,
\\
A_{1}  &  =\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(Z_{i}),\quad E[\psi
(Z_{i})]=0,\nonumber
\end{align}
where $A_{2}$ and $A_{3}$ are second and third-order products of sample averages of mean-zero random variables, each multiplied by $\sqrt{n}$; see
for example, Rilstone, Srivastava, and Ullah (1996). The terms $A_{1},$ $A_{2},$
and $A_{3}$ are bounded in probability so that equation (\ref{eq:expan})
is an expansion where the stochastic order of terms is smaller as one moves to
the right in the expansion.

We define the higher-order bias and variance of estimators using the first and second moments of the leading three terms in this expansion. We also compare higher-order efficiency of bias-corrected estimators by comparing their higher-order variance, as discussed in the introduction. The higher-order bias of $\sqrt{n}(\widehat{\theta}-$ $\theta_{0})$ is given by
$E[A_{2}]/\sqrt{n}.$ This follows from the fact that $E[A_{1}]=0$, while the expectation of $A_{3}/n$
is generally of smaller order than $E[A_{2}]/\sqrt{n}$. Dividing through
by $\sqrt{n},$ the higher-order bias of $\widehat{\theta}$ is
\[
Bias(\widehat{\theta})\approx\frac{b_{0}}{n},\text{ }b_{0}=E[A_{2}].
\]
In general $b_{0}=E[B(Z)],$ where the function $B(z)$ captures the `own observation' term in $A_{2},$ i.e.\ where the observation indices coincide in the
product of sample means that make up $A_{2}$. The term $B(z)$ may depend on the distribution of the
data through the parameter $\theta$ or in other ways. We let this dependence
be implicit for notational convenience.

Similarly, the higher-order variance of $\sqrt{n}(\widehat{\theta}-\theta_{0})$ can be obtained from the variance of the first three expansion terms, up to order $1/n.$ This gives
\begin{align*}
&Var\big(\sqrt{n}(\widehat{\theta}-\theta_{0})\big) \approx Var(A_1) + \frac{\upsilon}{n} \\
&\upsilon:=Var(A_{2})+2\sqrt{n}E[A_{1}
A_{2}]+2E[A_{1}A_{3}].
\end{align*}
where $Var(A_1)$ is the asymptotic (first-order) variance of the estimator and we refer to $\upsilon$ as the higher-order variance.




\medskip
{\bf Example 1: Normally distributed data.}
In order to demonstrate ideas, we will follow a simple example. Suppose that we observe a sample of observations $Z_{i}\sim N\left(
\sqrt{\theta},1\right) $. The MLE provides an efficient estimate for $\theta$, and is given by
$\widehat{\theta}=\left(  \frac{1}{n}\sum_{i=1}^{n}Z_{i}\right)  ^{2}$. It can be shown that the asymptotic expansion for this estimator is of the form in (\ref{eq:expan}), with
\begin{align} \label{eq:exp_example}
    A_1 = 2\sqrt{\theta} \frac{1}{\sqrt{n}}\sum_i( Z_i - \sqrt{\theta}),\qquad
    A_2 = \big(\frac{1}{\sqrt{n}} \sum_i (Z_i - \sqrt{\theta}) \big)^2,
\end{align}
and $A_3=0$.\footnote{See Section \ref{sec-4.2-proofs} in the Supplementary Appendix for details of this example.} In this case, since the estimator is quadratic, the expansion is exact. We can conclude that the bias of the estimator is $E[A_{2}]/n = 1/n$. Since $E[A_1A_2]=0$ in this example, the higher-order variance is
\[
Var\big(\sqrt{n}(\widehat{\theta}-\theta)\big) = Var(A_1) + Var(A_2)/n = 4\theta + 2/n.
\]

\subsection{The effect of bias correction}
A bias-corrected estimator can be formed by subtracting off an estimator $\widehat
{b}$ of $b_{0}$ from $\widehat{\theta}$
\[
\tilde{\theta}=\widehat{\theta}-\frac{\widehat{b}}{n}
\]
The focus of this paper is on the effect of the choice of $\widehat{b}$ on the
variance of $\tilde{\theta}.$ Generally $\widehat{b}$ has no effect on the
asymptotic variance (i.e. the first-order variance) of $\tilde{\theta}$, as long as $\widehat{b}$ is bounded in
probability, because then $\sqrt{n}(\widehat{b}/n)=\widehat{b}/\sqrt{n}=o_{p}(1).$ The
choice of $\widehat{b}$ can have an effect on the higher-order variance of
$\tilde{\theta}.$ To describe this effect note that the asymptotic expansion
for $\widehat{\theta}$ implies that, by adding and subtracting $b_{0}/\sqrt{n},$
\begin{align}\label{eq:exp_bias_corr}
\sqrt{n}(\tilde{\theta}-\theta_{0})  &  =A_{1}+\frac{A_{2}-b_{0}}{\sqrt{n}
}+\frac{A_{3}-\sqrt{n}(\widehat
{b}-b_{0})}{n}+o_{p}(n^{-1})\\
&  =A_{1}+\frac{\tilde{A}_{2}}{\sqrt{n}
}+\frac{\tilde{A}_{3}}{n}+o_{p}(n^{-1}),
\nonumber
\end{align}
where $\tilde{A}_{2}=A_2-b_0$, and $\tilde{A}_{3}=A_3-
 \sqrt{n}(\widehat
{b}-b_{0})$.
This is again an expansion whose terms are of decreasing stochastic order, so long as $\widehat{b}$ is
$\sqrt{n}$-consistent so that the third term is smaller order than the second term. Importantly, the second-order term  $\tilde{A}_{2}/\sqrt{n}$ has expectation zero, so that the bias-corrected estimator $\tilde{\theta}$ is higher-order unbiased.

To analyze the higher-order variance of the bias-corrected estimator it is
helpful to be more specific about $\widehat{b}.$ For now we assume $\widehat{b}$
is asymptotically equivalent to a sample average, i.e.\ that there exists a $\phi(z)$
with $E[\phi(Z)]=0$, $Var(\phi(Z))<\infty,$ such that
\[
\sqrt{n}(\widehat{b}-b_{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\phi(Z_{i})+o_{p}(1).
\]
Here we are assuming that $\widehat{b}$ is asymptotically linear with influence
function $\phi(z).$ Bias corrections that are based on directly estimating
$b_{0}=E[B(Z)]$ will often have this property and we find that other bias
corrections, like the jackknife and bootstrap also have this property. In this
case the expansion in equation (\ref{eq:exp_bias_corr}) continues to hold with
$\tilde{A}_{3}=A_{3}-\Delta$ for $\Delta:=\sum_{i=1}^{n}\phi(Z_{i})/\sqrt{n}.$
The variance of the third-order approximation for $\tilde{\theta}$ is $Var(A_1) + \tilde{\upsilon}/n$, with higher-order variance
\begin{align} \label{eq:higher-var}
\tilde{\upsilon}&:=Var(\tilde{A}_{2})+2\sqrt{n}E[A_{1}\tilde{A}_{2}]+2E[A_{1}\tilde{A}_{3}] \\
&= Var(A_{2}-b_0)+2\sqrt{n}E[A_{1}(A_{2}-b_0)]+2E[A_{1}(A_3 - \Delta)] \nonumber\\
&= \upsilon -2E[A_{1}\Delta] \nonumber.
\end{align}
Higher-order efficiency of the bias-corrected estimator refers to the size of $\tilde{\upsilon}$. One bias-corrected estimator is higher-order more efficient than another if it has smaller $\tilde{\upsilon}$.

The contribution of the bias correction to the higher-order variance $\tilde{\upsilon}$ is through the term
\begin{align} \label{eq:psi_phi_cov}
-2E[A_{1}\Delta]=-2E[\psi(Z)\phi(Z)].
\end{align}
Thus the bias correction $\widehat{b}$ affects the higher-order variance only
through the covariance of its influence function with that of $\widehat{\theta}$.
In the following theorem, we show that when $\tilde{\theta}$ is an efficient
estimator, its  higher-order variance does not depend on the choice of $\widehat{b}$, by demonstrating that the covariance in (\ref{eq:psi_phi_cov}) is the same for any $\phi(Z)$. We adopt the notation and terminology of Newey (1990, pp 104-106) for the statement and proof of this result. Here $\mathcal{T}$ denotes the tangent set for the semiparametric model, which is the mean-square closure of the set of all scores for regular parametric submodels. A parameter is differentiable if there is a random variable $d(Z)$ such that the derivative of the parameter with respect to parameters of any regular parametric submodel exists and equals the expected product of $d(Z)$ with the submodel score.  Under additional regularity conditions, the influence function of any asymptotically linear estimator will be equal to such a $d(Z)$.

\begin{theorem} \label{thm:semip_equiv}
If $\mathcal{T}$ is linear,  $\theta_0$ and $b_0$ are differentiable parameters of the semiparametric model,  $\widehat{\theta}$ is asymptotically linear and efficient with influence function $\psi(Z)$, and $\widehat{b}$ is asymptotically linear with influence function $d(Z)$, then the higher-order variance $\tilde{\upsilon}$ does not depend on $\widehat{b}$.
\end{theorem}

\begin{proof}

It follows from asymptotic linearity and efficiency of $\widehat{\theta}$ that its influence function $\psi(Z)$ is an element of $\mathcal{T}$. By Theorem 3.1 of Newey (1990) the efficient influence function $\delta(Z)$ for $b_0$ is equal to the projection of $d(Z)$ on $\mathcal{T}$. This allows us to decompose $d(Z)$ as,
\[
d(Z)=\delta(Z)+U(Z),\qquad E[U(Z)t(Z)]=0\text{ for all }
t\in\mathcal{T}\text{.}
\]
As shown above in (\ref{eq:psi_phi_cov}), the higher-order variance of $\tilde{\theta} = \widehat{\theta} - \widehat{b}/n$, depends on $\widehat{b}$ only through the covariance of $\psi(Z)$ with $d(Z)$. By the above decomposition, since $\psi(Z)$ is in the tangent set, this covariance is
\[
E[\psi(Z)d(Z)]=E[\psi(Z)\delta(Z)]+E[\psi(Z)U(Z)]=E[\psi(Z)\delta(Z)],
\]
Since the efficient influence function $\delta(Z)$ is unique, it follows
that $E[\psi(Z)d(Z)]$ does not vary with $d(Z)$, and hence $\tilde{\upsilon}$ does not vary with $\widehat{b}$.

\end{proof}

The proof uses geometry associated with semiparametric models, and relies on the efficiency of $\widehat{\theta}$.
There is also some relatively simple intuition for this result. Consider any two semiparametric estimators $\widehat{b}_{1}$ and $\widehat{b}_{2}$ of $b_{0}.$ The random variable $\widehat{b}_{1}
-\widehat{b}_{2}$ is a semiparametric estimator of $0$. If the asymptotic
covariance of $\widehat{\theta}$ with $\widehat{b}_{1}-\widehat{b}_{2}$ were nonzero then
for some fixed constant $C$ the asymptotic variance of $\bar{\theta}=\widehat{\theta
}+C(\widehat{b}_{1}-\widehat{b}_{2})$ would be less than the asymptotic variance of
$\widehat{\theta}$, so that $\widehat{\theta}$ could not be semiparametrically efficient.
That is, semiparametric efficiency of $\widehat{\theta}$ implies zero asymptotic
covariance of $\widehat{\theta}$ with $\widehat{b}_{1}-\widehat{b}_{2}$, which is
equivalent to the asymptotic covariance of $\widehat{\theta}$ with $\widehat{b}_{1}$
being equal to the asymptotic covariance of $\widehat{\theta}$ with $\widehat{b}_{2}.$
The asymptotic covariance of $\widehat{\theta}$ with any asymptotically linear
$\widehat{b}$ is $E[\psi(Z)\phi(Z)],$ so semiparametric efficiency of $\widehat
{\theta}$ implies that this covariance does not depend on $\phi(Z).$
This intuition extends the Hausman (1978) result that the covariance of any estimator of a parameter of interest with an efficient estimator of that parameter is equal to the variance of the efficient estimator, i.e.\ the covariance does not depend on the estimator. The intuition and result here show that the covariance of an efficient estimator of a parameter of interest with an estimator of any object does not depend on the estimator of that object.


\subsection{Asymptotically linear bias correction \label{subsec:2_3}}

There are many examples of bias corrections to semiparametric estimators that
have different influence functions. Consider a parametric model of the
conditional pdf of an outcome variable $Y$ given regressors $X$, where the pdf of
$Y$ conditional on $X$ has a parametric form $f(y|x,\beta)$ and $\theta$ is
some function of $\beta.$ This is a familiar semiparametric model where the
marginal distribution of $X$ is unspecified. A special case is a parametric
likelihood model where $X$ is a constant and $f(y|x,\beta)$ specifies the unconditional pdf
of observation $Y$. Since generally $b_{0}=E[B(Z)]$, for $Z=(Y,X)$ and for some function
$B(z)$ that may also depend on $\beta$, one could devise a number of estimates for $b_{0}$.
For example, an analytical bias correction could be obtained by plugging in estimates $\widehat{\beta}$ and replacing the expectation operators with
sample averages. Alternative forms of analytical bias correction could make use
of the structure imposed by $f(y|x,\beta)$ to integrate over $y$ in $B(z)$, or use restrictions implied by $\int f(y|x,\beta)dy=1$ to estimate $b_{0}$ (e.g. applying the information equality).

One could also use a nonparametric method like the bootstrap or jackknife to
estimate $b_{0}$. The bootstrap bias correction estimates the bias by $\frac{1}{n}\widehat{b}_B\equiv E^{\ast}\left[  \widehat{\theta}^{\ast}\right]
-\widehat{\theta}$, where $\widehat{\theta}^{\ast}$ are estimates obtained using bootstrap samples from the empirical distribution of $Z$. Alternatively, the jackknife estimates the bias by taking the difference between $\widehat{\theta}$ and the average of all `leave-one-out' estimates, i.e.\ estimates formed by excluding a single observation. The jackknife bias estimator is given by
$\frac{1}{n} \widehat{b}_J =(n-1)\big(\frac{1}{n}\sum_{i=1}^{n}
\widehat{\theta}_{\left(  i\right)  } - \widehat{\theta} \big)$
where $\widehat{\theta}_{(i)}$ is the estimate of the parameter that excludes observation $i$. In Section \ref{sec:mle_crosssec} we describe these bias correction techniques in detail in the context of MLE.

\medskip
{\bf Example 1, continued:}
Continuing our earlier example, we can compute the jackknife and analytical bias corrections for $\widehat{\theta}$. From the asymptotic expansion in (\ref{eq:exp_example}), we have that
\[
    E[A_2] = E[(Z_i - \sqrt{\theta})^2]
\]
which suggests the analytical bias estimate $\widehat{b}_a = \frac{1}{n}\sum_i (Z_i  - \bar{Z})^2$. Alternatively, using the known variance, we could simply use $b_a=1$. To construct a jackknife bias estimator, we use
\[
    \frac{1}{n} \widehat{b}_J = (n-1)\Big( \frac{1}{n}\sum_i \big(\frac{1}{n-1} \sum_{j\ne i} Z_j \big)^2 - \big(\frac{1}{n} \sum_{i} Z_i \big)^2 \Big)
    = \frac{1}{n(n-1)}\sum_i (Z_i  - \bar{Z})^2
\]
which is the same as the first analytical bias estimate, up to the factor $n/(n-1)$. It is straightforward to show that both bias estimates can be written as
\[
    \sqrt{n}(\widehat{b}-1) = \frac{1}{n}\sum_i \big((Z_i  - \sqrt{\theta})^2 -1\big) + o_p(1)
\]
so that they are both asymptotically linear, in this case with the same influence function $d(Z)= (Z  - \sqrt{\theta})^2 -1$. Given that the third moment of $Z$ is equal to zero, since it is normally distributed, we have $E[A_1(\widehat{b}-1)]=0$, so that the higher-order variance of the bias-corrected estimators that use either of the analytical or the jackknife bias estimates, are the same as that of $\widehat{\theta}$, i.e. $4\theta + 2/n$.
\medskip


Theorem \ref{thm:semip_equiv} implies that the various analytical bias corrections as well as both nonparametric bias correction methods
lead to the same higher-order variance as any other method when the $\widehat{b}$
is asymptotically linear. In the example above, both the analytical and jackknife bias estimates are asymptotically linear. This is not a coincidence. It turns out that for a fairly general class of models, the analytical, jackknife, and bootstrap bias estimates are asymptotically linear, so that Theorem \ref{thm:semip_equiv} can be applied.\footnote{The theorem is stated for a scalar parameter $\theta$, but straightforwardly extends to a vector-valued parameter, at the cost of more complicated notation.}

\begin{theorem} \label{thm:asym_lin}
Let $\widehat{\theta}$ be an estimator with a stochastic expansion
\begin{equation} \label{eq:expansion_asym_lin}
\widehat{\theta}-\theta_{0} = \frac{1}{\sqrt{n}}A_{1}+\frac{1}{n}A_{2}+\frac{1}{n^{3/2}}A_{3}
 +\frac{1}{n^2}A_{4}+\frac{1}{n^{5/2}}A_{5}+o_{p}(n^{5/2})
\end{equation}
where $n^{-k/2}A_{k}$ is a $k$-th order V-statistic, i.e.\
\begin{equation}
\frac{1}{n^{k/2}}A_{k}=\frac{1}{n^k}\sum_{i_1=1}^{n}\cdots \sum_{i_k=1}^{n}
g_{k,1}(z_{i_1},\theta_0)\cdots g_{k,k}(z_{i_k},\theta_0)
\end{equation}
and $g_{k,j}$ are functions of the data that are continuously differentiable in $\theta$ at $\theta_0$, mean zero, i.e.\ $E[g_{j,k}(z_{i},\theta_0)]=0$, and with $E[g_{j,k}(z_i,\theta_0)^{10}]\leq C <\infty$.

Then, the jackknife, bootstrap, and (sample average) analytical\footnote{See Section \ref{sec:mle_crosssec} for exact definitions of these three bias corrections, $\widehat{\theta}_J$, $\widehat{\theta}_B$, and $\tilde{\theta}_a$.} bias estimates of the higher-order
bias $b_{0}=E[A_2]=E[g_{2,1}(z_{i},\theta_0)g_{2,2}(z_{i},\theta_0)]$, are asymptotically linear.
\end{theorem}

\begin{proof}
See Appendix \ref{app:main}
\end{proof}

The V-statistic structure assumed in Theorem \ref{thm:asym_lin} is shown to hold for the case of MLE under standard regularity conditions on the likelihood function (see Supplementary Appendix I \ref{Proof-Main-Results}). In Section 4 we derive such expansions for a maximum likelihood estimator. The derivation is based on an expansion of the first order condition of the MLE and so applies to any other estimator that can be characterized by a moment equation satisfying similar regularity conditions. Therefore, the structure can be shown to appear in the expansion of any M-estimator, i.e. an estimator that maximizes the sample average of some function, which need not be a log-likelihood.

An important semiparametric model is a model of unconditional moment
restrictions that motivates GMM estimation. Newey and Smith (2004) derive the asymptotic expansion for generalized empirical likelihood (GEL) estimators up to third order and examination of that expansion shows the same V-statistic structure. Under sufficient regularity conditions, the structure should also exist in the fifth-order expansion used in Theorem \ref{thm:asym_lin}, so that the equivalence result of Theorem \ref{thm:semip_equiv} would imply that jackknife and bootstrap bias corrections lead to the same higher-order variance for GEL estimators. This extends the result in Newey and Smith (2004), who found that
averaging over an efficient estimator of the distribution of the data to
obtain $\widehat{b}$ does not affect the higher-order variance of bias-corrected
GMM and generalized empirical likelihood (GEL) estimators.


Another interesting example is a nonparametric model where no restrictions are placed
on the distribution of the data. In this case $b_{0}$ will also be a
nonparametric object and there will exist only one influence function
corresponding to $b_{0}$; see Van der Vaart (1991) and Newey (1994). Here, all
asymptotically linear estimators of $b_{0}$ will have the same influence
function so that Theorem \ref{thm:semip_equiv} holds trivially. In particular, when general
misspecification is allowed for in debiasing, so that the distribution of
the data is unrestricted, there will only be one influence function for
$b_{0}$. Then any bias correction method, such as the analytical, bootstrap, or jackknife,
must have the same higher-order variance, as long as they are asymptotically linear.

We emphasize that Theorem \ref{thm:semip_equiv} only applies to comparisons between bias
corrections for a given semiparametrically efficient estimator. The stochastic
expansion term $A_{2}$ may differ across semiparametrically efficient estimators so that their
higher-order variances need not be the same.
For example, Newey and Smith (2004) showed that GMM will tend to have larger second-order bias than GEL when there are many instrumental variables, so that $A_{2}$ is different for GMM
and GEL. Theorem \ref{thm:semip_equiv} has nothing to say about how the higher-order variance of
bias-corrected GMM compares with that of bias-corrected GEL. It only implies that the form
of $\widehat{b}$ does not affect the higher-order variance of bias-corrected GMM
and separately for bias-corrected GEL.

Note that while Theorem \ref{thm:asym_lin} implies that bootstrap, jackknife, and analytical bias estimates are asymptotically linear in general, this is not true of all bias estimates. An example of this is the split-sample jackknife, which estimates the bias using $\frac{1}{n}\widehat{b}_{SS} = \frac{1}{2}(\widehat{\theta}_{(1)} + \widehat{\theta}_{(2)}) - \widehat{\theta}$, where $\widehat{\theta}_{(1)}$ and $\widehat{\theta}_{(2)}$ are estimates using two separate halves of the data set.

{\bf Example 1, continued:}
Let $N=2m$, and define $\xi^{(1)} = \frac{1}{\sqrt{m}} \sum_{i\leq m} (Z_i-\sqrt{\theta})$ and $\xi^{(2)} = \frac{1}{\sqrt{m}} \sum_{i>m} (Z_i-\sqrt{\theta})$.  Some algebra gives the expression
\begin{align*}
    \widehat{b}_{SS} &= \frac{1}{2} \big( \xi^{(1)} - \xi^{(2)} \big)^2
\end{align*}
Since $\xi^{(1)}$ and $\xi^{(2)}$ are independent standard normal variables, we have
$E[\widehat{b}_{SS}]=1$, so that the split-sample jackknife gives an unbiased estimator for $b_0$. However, $\widehat{b}_{SS}=O_p(1)$, so that the bias estimator is also inconsistent, and hence not asymptotically linear. In this case, the higher-order variance of the bias-corrected estimator can be shown to be $Var\big(\sqrt{n}(\widehat{\theta}_{SS})\big) = 4\theta + 4/n$, which is larger than that of the analytical and jackknife estimators. In Section \ref{sec:SS_var}, we give a formal result that shows that this larger higher-order variance is true in general for MLE.








\section{Bias-corrected MLE\label{sec:mle_crosssec}}

In this section, we demonstrate the higher-order equivalence of bias corrections in a fully parametric model by deriving higher-order expansions and variance expressions for several bias-corrected maximum likelihood estimators. We consider the bootstrap bias correction, jackknife bias correction,
as well as three different versions of analytical bias correction, including the bias-corrected MLE of Pfanzagl and Wefelmeyer (1978), which was shown to be third-order optimal.
We show that all of these corrections lead to estimators with the same higher-order variance. Further, we show that the bootstrap, jackknife, and a particular analytical bias-correction, result in identical third-order asymptotic expansions.

To describe the parametric model, let $\left\{  Z_{i}\right\}  _{i=1}^{n}$ be an i.i.d.\ sample
$Z_{i}\sim f\left(  z,\theta_{0}\right)  $, such that $f(z,\theta)$ satisfies
sufficient smoothness conditions.\footnote{The results in Section
\ref{sec:mle_crosssec} are predicated on regularity conditions,
including Conditions \ref{HH}, \ref{M} and \ref{EC}, which are presented in
Appendix \ref{sec:reg-conds}.} The density $f(z,\theta)$ is a member of
a parametric family of distributions $P_{\theta}$ indexed by $\theta\in\Theta$
with $\Theta\in\mathbb{R}$ a compact set.\footnote{For simplicity of notation,
we will assume that $p=1$, where $p=\dim\left(  \theta\right)  $. The result
is expected to hold for any finite $p$.} We consider properties of the
MLE $\widehat{\theta}$, where
\[
\widehat{\theta}\equiv\arg\sup_{\theta\in\Theta}n^{-1}\sum_{i=1}^{n}\log f\left(
Z_{i},\theta\right)  .
\]


We adopt the following set of notations in this section of the paper. We let $\ell\left(
\cdot,\theta\right)  \equiv\left.  \partial\log f\left(  \cdot,\theta\right)
\right/  \partial\theta$, $\ell^{\theta}\left(  \cdot,\theta\right)
\equiv\left.  \partial^{2}\log f\left(  \cdot,\theta\right)  \right/
\partial\theta^{2}$, $\ell^{\theta\theta}\left(  \cdot,\theta\right)
\equiv\left.  \partial^{3}\log f\left(  \cdot,\theta\right)  \right/
\partial\theta^{3}$, etc. We also define $\mathcal{I}\equiv-E\left[
\ell^{\theta}(Z_{i},\theta_{0})\right]  $, $\mathcal{Q}_{1}\left(
\theta\right)  \equiv E\left[  \ell^{\theta\theta}(Z_{i},\theta)\right]  $,
and $\mathcal{Q}_{2}\left(  \theta\right)  \equiv E\left[  \ell^{\theta
\theta\theta}(Z_{i},\theta)\right]  $. Furthermore, we let $U_{i}\left(
\theta\right)  \equiv\ell\left(  Z_{i},\theta\right)  $, $V_{i}\left(
\theta\right)  \equiv\ell^{\theta}\left(  Z_{i},\theta\right)  -E\left[
\ell^{\theta}\left(  Z_{i},\theta\right)  \right]  $, $W_{i}(\theta) \equiv\ell
^{\theta\theta}\left(  Z_{i},\theta\right)  -E\left[  \ell^{\theta\theta}\left(
Z_{i},\theta\right)  \right]  $. Finally, we define $U\left(  \theta\right)  \equiv
n^{-1/2}\sum_{i=1}^{n}U_{i}\left(  \theta\right)  $, $V\left(  \theta\right)
\equiv n^{-1/2}\sum_{i=1}^{n}V_{i}\left(  \theta\right)  $, and $W\left(
\theta\right)  \equiv n^{-1/2}\sum_{i=1}^{n}W_{i}\left(  \theta\right) $.
We next describe the bias corrections and derive higher-order properties of the bias-corrected estimators.

\subsection{The bootstrap and jackknife bias corrections}

The bootstrap estimator constructs an estimate of the higher-order bias as the difference between the average of bootstrap replicates $\widehat{\theta}^{\ast}$ and the MLE $\widehat{\theta}$. In particular, we first obtain bootstrapped estimates $\widehat{\theta}^{\ast}$ by sampling
$Z_{1}^{\ast},...,Z_{n}^{\ast}$ identically and independently from the
empirical distribution $\widehat{F}(z) = \frac{1}{n}\sum 1\{ Z_i \leq z\}$. Let $E^{\ast}$ be the expectation
operator with respect to $\widehat{F}.$ The idea behind the bootstrap bias
correction is to estimate $E\left[  \widehat{\theta}\right]  -\theta_{0},$ if
it exists, by $\frac{1}{n}\widehat{b}_B\equiv E^{\ast}\left[  \widehat{\theta}^{\ast}\right]
-\widehat{\theta}$. This in turn will allow us to construct the bias-corrected estimate
\[
\widehat{\theta}_{B}\equiv\widehat{\theta}-\frac{\widehat{b}_B}{n}=2\widehat{\theta}-E^{\ast
}\left[  \widehat{\theta}^{\ast}\right]  .
\]

An alternative nonparametric bias correction is a jackknife bias-corrected estimator. The jackknife estimates the higher-order bias by taking the difference between the MLE and the average of all `leave-one-out' estimates, i.e.\ estimates formed by excluding a single observation. The jackknife estimate is given by
\[
\widehat{\theta}_{J} = \widehat{\theta} - \frac{\widehat{b}_J}{n}
=n\widehat{\theta}-\frac{n-1}{n}\sum_{i=1}^{n}
\widehat{\theta}_{\left(  i\right)  }
\]
where $\widehat{\theta}_{(i)}$ is the estimate of the parameter that excludes observation $i$. The jackknife uses the bias estimate $\frac{1}{n}\widehat{b}_J = (n-1)\left(\frac{1}{n}\sum_{i=1}^{n}\widehat{\theta}_{\left(  i\right)} - \widehat{\theta} \right)$.


The following proposition establishes the higher-order properties of the
bootstrap and jackknife bias-corrected MLEs.

\begin{proposition}
\label{boot-jack-higher}Let $\tilde{b}$ be either the bootstrap bias estimate $\widehat{b}_B$, or the jackknife bias estimate $\widehat{b}_J$. Then, under regularity conditions stated in Appendix \ref{sec:reg-conds}, $\tilde{b}$ satisfies
\begin{align*}
    \sqrt{n}(\tilde{b} - b_0) &= \frac{1}{\sqrt{n}}\sum_{i=1}^n \mathbb{B}_i + o_p(1)
\end{align*}
where
\begin{align*}
    \mathbb{B}_i &= \left( \frac{1}{2}\mathcal{I}^{-3}\mathcal{Q}_{2} + \frac{3}{2}\mathcal{I}^{-4}\mathcal{Q}_{1}^2 + 3\mathcal{I}^{-4}\mathcal{Q}_{1}E[U_iV_i] + \mathcal{I}^{-3}\left(E[U_iW_i] + E[(\ell^\theta)^2] \right) \right) U_i \\
    &+ \left( \frac{3}{2}\mathcal{I}^{-3}\mathcal{Q}_{1} + 2\mathcal{I}^{-3}E[U_iV_i]  \right) V_i
    + \frac{1}{2}\mathcal{I}^{-2} W_i \\
    &+ \frac{1}{2}\mathcal{I}^{-3}\mathcal{Q}_{1}\left(   \ell(Z_i,\theta_0)^2 - E[\ell^2]\right)
    + \mathcal{I}^{-2} \left( \ell(Z_i,\theta_0) \ell^\theta(Z_i,\theta_0) - E[\ell \ell^\theta] \right)
\end{align*}
and hence, for $\tilde{\theta}$ either the bootstrap or jackknife bias-corrected MLE,
\[
\sqrt{n}\left(  \tilde{\theta}-\theta_{0}\right)  =A_{1}+\frac{1}
{\sqrt{n}}\left(  A_{2}-b(\theta_{0})\right)  +\frac{1}{n}(A_{3}-
\mathbb{B})+o_{p}\left( n^{-1} \right)  ,
\]
where $\mathbb{B} = \sum_i \mathbb{B}_i / \sqrt{n}$.
\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}\footnote{The Supplementary Appendix is available at \href{http://arxiv.org/abs/2207.09943}{http://arxiv.org/abs/2207.09943}}.
\end{proof}


Here we see that the jackknife and bootstrap bias estimates are $\sqrt{n}$-consistent and asymptotically linear (note that $E[\mathbb{B}_i]=0$), which from Theorem  \ref{thm:semip_equiv} implies that their higher-order variances are the same as any other estimator that uses an asymptotically linear estimator of the bias. In fact, the jackknife and bootstrap share an even stronger equivalence, in the sense that their bias estimates have identical influence functions, which implies that their asymptotic expansions are identical up to the third order.

\subsection{Analytical bias corrections}

We next introduce and discuss three forms of analytical bias correction. An analytical expression for the higher-order bias is given by
\begin{align}
b\left(  \theta_{0}\right)  \equiv E \left[  A_{2}\right]   &
=\frac{1}{2\mathcal{I}^{3}}E\left[  \ell^{\theta\theta}\right]  E\left[
\ell^{2}\right]  +\frac{1}{\mathcal{I}^{2}}E\left[  \ell\ell^{\theta}\right]
\label{bias-no-info}
\end{align}

The first bias correction is based on the bias formula (\ref{bias-no-info}), replacing expectations with sample averages. This gives the estimator $\widetilde{\theta}_{a}\equiv\widehat{\theta}-\frac{\widetilde{b}\left(
\widehat{\theta}\right)  }{n}$, where
\begin{align} \label{analy-no-info}
\widetilde{b}\left(\widehat{\theta}\right) =
-\frac{\left(  \frac{1}{n}\sum_{i}\ell^{\theta\theta}(  Z_{i},\widehat{\theta
})  \right)  \left(  \frac{1}{n}\sum_{i} \ell(  Z_{i}
,\widehat{\theta})   ^{2}\right)  }{2\left(  \frac{1}{n}\sum_{i}
\ell^{\theta}(  Z_{i},\widehat{\theta})  \right)  ^{3}}
+\frac{\left(  \frac{1}{n}\sum_{i}\ell(  Z_{i},\widehat{\theta})
\ell^{\theta}(  Z_{i},\widehat{\theta})  \right)  }{\left(
\frac{1}{n}\sum_{i}\ell^{\theta}(  Z_{i},\widehat{\theta})  \right)
^{2}}
\end{align}

It can be shown that this form of analytical bias estimate shares the same influence function as the bootstrap and jackknife bias estimators. It follows that it has an identical third-order asymptotic expansion.

\begin{proposition}
\label{HigherOrderEfficiency3} Let the regularity conditions stated in Appendix \ref{sec:reg-conds} hold. The estimator $\widetilde{\theta}_a$ has the higher-order expansion
\[
\sqrt{n}\left(  \widetilde{\theta}_{a}-\theta_{0}\right)  = A_{1}+\frac
{1}{\sqrt{n}}\left(  A_{2}-b\left(  \theta_{0}\right)  \right)  +\frac{1}
{n}\left(  A_{3}-\mathbb{B}\right)  +o_{p}\left(  n^{-1}\right)  .
\]

\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}.
\end{proof}

Alternative analytical corrections can be constructed by imposing the information equality, $E[\ell^2]=\mathcal{I}$ in the characterization of the bias, that is, using the bias formula
\begin{align}
b\left(  \theta_{0}\right)  \equiv E \left[  A_{2}\right]
  =\mathcal{I}^{-2} \left( \frac{1}{2}E\left[  \ell^{\theta\theta}\right]  + E\left[  \ell\ell^{\theta}\right] \right). \label{bias-info}
\end{align}

An analytical bias correction that uses the information equality is given by
\[
\widehat{\theta}_{a}\equiv\widehat{\theta}-\frac{\widehat{b}\left(
\widehat{\theta}\right)  }{n}\equiv\widehat{\theta}-\frac{1}{n}\left(
\frac{\left(  n^{-1}\sum_{i}\ell^{\theta\theta}(  Z_{i},\widehat{\theta
})  \right)  }{2\left(  n^{-1}\sum_{i}\ell^{\theta}(
Z_{i},\widehat{\theta})  \right)  ^{2}}+\frac{\left(  n^{-1}\sum_{i}
\ell(  Z_{i},\widehat{\theta})  \ell^{\theta}(  Z_{i}
,\widehat{\theta})  \right)  }{\left(  n^{-1}\sum_{i}\ell^{\theta
}(  Z_{i},\widehat{\theta})  \right)  ^{2}}\right)  .
\]

It is also possible to construct a bias correction that computes the expectations in (\ref{bias-info}) in their integral form. This is the analytical bias-corrected estimator of Pfanzagl and Wefelmeyer (1978), which they showed is
higher-order efficient.
\[
\widehat{\theta}_{c}\equiv\widehat{\theta}-\frac{b\left(  \widehat{\theta
}\right)  }{n}=\widehat{\theta}-\frac{1}{n}\left(  \frac{\int\ell
^{\theta\theta}(  z,\widehat{\theta})  f(  z,\widehat{\theta
})  dz}{2\left(  \int\ell^{\theta}(  z,\widehat{\theta})
f(  z,\widehat{\theta})  dz\right)  ^{2}}+\frac{\int\ell(  z,\widehat{\theta})  \ell^{\theta}(  z,\theta
)  f(  z,\widehat{\theta})  dz}{\left(  \int\ell^{\theta
}(  z,\widehat{\theta})  f(  z,\widehat{\theta})
dz\right)  ^{2}}\right)  .
\]

The analytical bias-corrected estimators that impose the information equality do not have identical higher-order expressions to the bias-corrected estimators that do not rely on the information equality. Instead, they differ in their third expansion terms, as can be seen in the following proposition.

\begin{proposition}
\label{HigherOrderEfficiency}Let the regularity conditions stated in Appendix \ref{sec:reg-conds} hold. The analytical  bias estimates satisfy
\begin{align*}
    \sqrt{n}(\widehat{b}(\widehat{\theta})  - b_0) &= \frac{1}{\sqrt{n}}\sum_{i=1}^n \mathbb{A}_i + o_p(1) \\
    \sqrt{n}(b(\widehat{\theta}) - b_0) &= \frac{1}{\sqrt{n}}\sum_{i=1}^n \mathbb{C}_i + o_p(1)
\end{align*}
where $\mathbb{A}_i$ and $\mathbb{C}_i$ are as defined in Supplementary Appendix I \ref{Proof-Main-Results}. Hence, the analytical bias-corrected estimators $\widehat{\theta}_{c}$ and $\widehat{\theta}_{a}$ have higher-order expansions
\begin{align*}
\sqrt{n}\left(  \widehat{\theta}_{c}-\theta_{0}\right)
&  =A_{1}+\frac{1}{\sqrt{n}}\left(  A_{2}-b\left(  \theta_{0}\right)  \right)  +\frac{1}
{n}\left(  A_{3}-\mathbb{C}\right)  +o_{p}\left(  n^{-1}\right)  ,\\
\sqrt{n}\left(  \widehat{\theta}_{a}-\theta_{0}\right)
&  =A_{1}+\frac{1}{\sqrt{n}}\left(  A_{2}-b\left(  \theta_{0}\right)  \right)  +\frac{1}
{n}\left(  A_{3}-\mathbb{A}\right)  +o_{p}\left(  n^{-1}\right)  ,
\end{align*}
where $\mathbb{A} = \sum_i \mathbb{A}_i /\sqrt{n}$ and $\mathbb{C} = \sum_i \mathbb{C}_i /\sqrt{n}$.
\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}.
\end{proof}



Propositions \ref{boot-jack-higher},  \ref{HigherOrderEfficiency3}, and
\ref{HigherOrderEfficiency} imply that the
higher-order variances of $\widehat{\theta}_{B}$, $\widehat{\theta}_{J}$,
$\widetilde{\theta}_{a}$ are equal to $\operatorname{Var}\left(  A_{2}\right)
+2n^{1/2}E\left[  A_{1}A_{2}\right]  +2E\left[  A_{1}\left(  A_{3}
-\mathbb{B}\right)  \right]  $,
while the higher-order variances of $\widehat{\theta}_{a}$, and
$\widehat{\theta}_{c}$ are $\operatorname{Var}\left(  A_{2}\right)
+2n^{1/2}E\left[  A_{1}A_{2}\right]  +2E\left[  A_{1}\left(  A_{3}
-\mathbb{A}\right)  \right]  $ and $\operatorname{Var}\left(  A_{2}\right)
+2n^{1/2}E\left[  A_{1}A_{2}\right]  +2E\left[  A_{1}\left(  A_{3}
-\mathbb{C}\right)  \right]  $, respectively.
Therefore, it is natural to
conjecture that the higher-order variances of $\widehat{\theta}_{B}$,
$\widehat{\theta}_{J}$, and $\widetilde{\theta}_{a}$ are different from that
of $\widehat{\theta}_{c}$ or $\widehat{\theta}_{a}$.
However, the propositions also state that each estimator uses an estimate of the higher-order bias $b_0$ that is asymptotically linear, and hence a corollary of Theorem \ref{thm:semip_equiv} is that the higher-order variances of the bias-corrected estimators are all equal.

\begin{corollary}\label{ET1T3}
Let the regularity conditions stated in Appendix \ref{sec:reg-conds} hold. Then,
\[
    E\left[  \mathbb{B}A_{1}\right]  =E\left[
\mathbb{A}A_{1}\right]  =E\left[  \mathbb{C}A_{1}\right]
\]
\end{corollary}


\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}.
\end{proof}

By Corollary \ref{ET1T3}, and the expression for the higher-order variance (\ref{eq:higher-var}), we conclude that all five estimators are
higher-order efficient. This result follows because (i) they have identical higher-order variance; and (ii) $\widehat{\theta}_{c}$ was shown to be higher-order efficient by Pfanzagl and Wefelmeyer (1978).

The next proposition provides an expression for the higher-order variance of these five bias-corrected estimators.

\begin{proposition} \label{higher_var}
Let the regularity conditions stated in Appendix \ref{sec:reg-conds} hold, and define
\begin{align*}
    X_{i} &= \mathcal{I}^{-1}U_i \\
    Y_{i} &= \frac{1}{2}\mathcal{I}^{-2}\mathcal{Q}_1U_i + \mathcal{I}^{-1}V_i \\
    \Upsilon &= \left( E[X_{i}^2]E[Y_{i}^2] + E[X_{1}Y_{i}]^2\right)
\end{align*}
The higher-order variance of $\widehat{\theta}_B$, $\widehat{\theta}_J$, $\widetilde{\theta}_a$, $\widehat{\theta}_a$, and $\widehat{\theta}_c$ is
\begin{align} \label{eq:higher-var-expr}
    Var(A_1) + \frac{\tilde{\upsilon}}{n} =
    \mathcal{I}^{-1} + \frac{1}{n}\Upsilon
\end{align}
\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}.
\end{proof}

The above higher-order variance expression is useful since the higher-order variance formula in (\ref{eq:higher-var}) includes two covariance terms $2E\left[  A_{1}
A_{3}\right]  +2n^{1/2}E\left[  A_{1}A_{2}\right]  $, and as such, it is not
obvious that (\ref{eq:higher-var}) is positive. On the other hand, because
$\Upsilon>0$, the expression in the proposition consists of
positive terms and eliminates this concern. According to Theorem \ref{thm:semip_equiv},
(\ref{eq:higher-var-expr}) is also the higher-order variance of any other asymptotically linear bias correction.










\subsection{Inefficient bias correction} \label{sec:SS_var}

Our discussion above shows that many common methods of bias correction use estimates of the bias that are asymptotically linear, and hence estimate the bias at the $n^{-1/2}$-rate. These bias-corrected estimators have equivalent higher-order variances.
We now show by an example that the equivalence
result does not hold in general if this rate assumption is violated. We showed
in Example 1 (see the end of Section
\ref{subsec:2_3}) that the split-sample jackknife provided an inconsistent estimate of the higher-order bias, so that its
higher-order variance was larger than that of the leave-one-out jackknife. In the Appendix we show that this result holds more generally, for estimators with the expansion structure used in Theorem 2. Here we derive expressions for the higher-order variance in the MLE case.\footnote{Derivation of the expressions used below are available in Section \ref{HOV-SS} of Supplementary Appendix I.}

Recall that the split-sample jackknife estimator is given by \[\widehat{\theta}_{SS}\equiv2\widehat{\theta}-\left(  \widehat{\theta}_{(1)}+\widehat{\theta}_{(2)}\right)  /2,
\]
where $\widehat{\theta}_{(  1)  }$ and $\widehat{\theta}_{(  2)}$ denote the MLEs based on separate halves of the sample. The implicit bias estimate used by the split-sample jackknife is given by $\frac{1}{n}\widehat{b}_{SS}=(  \widehat{\theta}_{(1)}+\widehat{\theta}_{(2)})  /2 - \widehat{\theta}$.

\begin{proposition} \label{higher_var_ss}
Let the regularity conditions stated in Appendix \ref{sec:reg-conds} hold. The bias estimate used by the split-sample estimator can be written as
\begin{align*}
    \widehat{b}_{SS} - b_0 &= \frac{1}{2}\left( X_{(1)}Y_{(1)} + X_{(2)}Y_{(2)} - 2b_0 -\left( X_{(1)}Y_{(2)} + X_{(2)}Y_{(1)}  \right) \right) \\
    &= O_p(1)
\end{align*}
where, for $n=2m$, $X_{(1)} = \frac{1}{\sqrt{m}}\sum_{i=1}^m X_{i}$ is the scaled sum over $X_{i}$ in the first half of the sample, and $X_{(2)}$, $Y_{(1)}$, and $Y_{(2)}$ are defined similarly (for $X_i$ and $Y_i$ as defined in Proposition \ref{higher_var}).

The higher-order variance of $\widehat{\theta}_{SS}$ is
\begin{align*}
    Var(A_1) + \frac{\tilde{\upsilon}_{SS}}{n} =
    \mathcal{I}^{-1} + \frac{2}{n}\Upsilon
\end{align*}
\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{Proof-Main-Results}.
\end{proof}

That is, the higher-order variance of $\widehat{\theta}_{SS}$ is strictly larger than that of $\widehat{\theta}_{J}$. As shown in Proposition \ref{higher_var_ss}, the split-sample bias correction provides an inconsistent estimate of the bias term $b_0$, which results in a bias correction that impacts the second-order term $A_2$ in the expansion, and hence affects the higher-order variance of the estimator.










\section{Bias correction for panel data\label{section-panel-model}}

We now extend some of these ideas to the panel data setting by comparing the jackknife and split-sample bias correction methods for a model with individual fixed effects. We begin with a description of the fixed effects maximum likelihood estimator.
Let $Z_{it}$, for $i=1,\dots,n$ and $t=1,\dots,T$, be a vector of observed
data. Denote $\theta$ a $p\times1$ parameter vector and $\alpha_{i}$ a scalar
unobserved individual effect.\footnote{As before, we will assume that $p=1$
for notational simplicity, but results should be expected to hold for any finite
$p>1$. The analysis in this section assumes that time effects are not present.} The data have density function $f(z|\theta,\alpha)$ with respect to
some measure, and so (treating the $\alpha_{i}$ as parameters to be estimated)
we may estimate $\theta$ via maximum likelihood. Assuming that the $Z_{it}$
are independent across both $i$ and $t$, the MLE solves
\[
\widehat{\theta}_{T}\equiv\arg\max_{\theta}\sum_{i=1}^{n}\sum_{t=1}^{T}\ln
f(Z_{it}|\theta,\widehat{\alpha}_{i}(\theta)),\quad\widehat{\alpha}_{i}(\theta
)\equiv\arg\max_{\alpha}\sum_{t=1}^{T}\ln f(Z_{it}|\theta,\alpha)
\]

In contrast to the higher-order bias discussed earlier for the cross-sectional model, the fixed effects panel data estimator suffers from the well-known incidental parameters problem (Neyman and Scott, 1948). For fixed $T$, the probability limit of the estimator $\theta_{T}\equiv
\operatorname*{plim}_{n\rightarrow\infty}\widehat{\theta}_{T}$ generally differs from
$\theta_{0}$.\footnote{$\theta_{T}$ is
given by $\operatorname*{argmax}_{\theta}\lim_{n\rightarrow\infty}\frac{1}
{n}\sum_{i=1}^{n}E\left[  \sum_{t=1}^{T}\ln f(z_{it}|\theta,\widehat{\alpha}
_{i}(\theta))\right]  $} Even when the number of time periods grow, so that $n/T\rightarrow\rho$ (which will be assumed in this
section), the estimator $\widehat{\theta}$ remains asymptotically biased, i.e.\
$\sqrt{nT}(\widehat{\theta}-\theta_{0})\Rightarrow N(\sqrt{\rho}\mathbf{B}
,\Omega)$ for some bias term $\mathbf{B}$. The bias is of order $O(T^{-1})$,
and can be substantial if $T$ is not sufficiently large.

It is useful to think about $\theta_{0}$ and $\alpha_{i}$ as solutions to a
set of moment equations given by the score functions
\[
0=\sum_{i=1}^{n}E\left[  \frac{\partial}{\partial\theta}\ln f(z_{it}
|\theta_{0},\alpha_{i})\right]  ,\quad0=E\left[  \frac{\partial}
{\partial\alpha_{i}}\ln f(z_{it}|\theta_{0},\alpha_{i})\right]
\]
As earlier, we can expand these first order
conditions to produce asymptotic expansions. We may also consider other quantities of interest that can be
defined via some moment condition, for example an average effect parameter
$\mu_{0}$, defined as the solution to
\[
0=\mu_{0}-\frac{1}{n}\sum_{i=1}^{n}E\left[  m(z_{it},\theta_{0},\alpha
_{i})\right]
\]
where $m$ is some function of interest, for example, the partial derivative of a conditional expectation function. Stacking these moment conditions, we can define the common parameter to be
$(\theta,\mu)$, so that the results presented below will also apply to these types of parameters.

We use the following notation for this section. Let $u_{it}(\theta
,\alpha)\equiv\frac{\partial}{\partial\theta}\ln f(z_{it}|\theta,\alpha)$ and
$V_{it}(\theta,\alpha)\equiv\frac{\partial}{\partial\alpha_{i}}\ln
f(z_{it}|\theta,\alpha)$ be the score functions. When evaluating functions at
the true value of parameters, arguments will be dropped, e.g. $u_{it}
=u_{it}(\theta_{0},\alpha_{i})$. Further, let $U_{it}(\theta,\alpha
)=u_{it}(\theta,\alpha)-\delta V_{it}(\theta,\alpha)$, for $\delta
=E[u_{it}V_{it}]/E[V_{it}^{2}]$, be the efficient score for $\theta$. All
expectations are taken with respect to the distribution for an individual $i$,
that is $E[h(z_{it})]=\int h(z)f(z|\theta,\alpha_{i})dz$. Also define $\mathcal{I}_n = \frac{1}{n}\sum_{i=1}^n E[U_{it}^2]$. As in the cross
sectional case, we denote partial derivatives of these functions with
superscripts, e.g. $\partial U_{it}/\partial\theta=U_{it}^{\theta}$.




\subsection{Higher-order comparison of jackknife estimators\label{jackknife-discussion}}

In this section we derive the higher-order properties of both the leave-one-out jackknife and split-sample jackknife estimators for the panel data model.

It has be shown previously (e.g. Hahn and Kuersteiner (2002), Hahn and Newey (2004)) that the
MLE $\widehat{\theta}$ has an expansion of the form
\begin{equation}
\sqrt{nT}(\widehat{\theta}-\theta_{0})= A_1 +\frac{\sqrt{n}}{\sqrt{T}}A_2+\frac{\sqrt{n}}{T}A_3+O_{p}(T^{-1}) \label{eq:mle_exp}
\end{equation}
where the expansion terms $A_1,A_2,A_3$ are each $O_p(1)$.\footnote{In Supplementary Appendix IV, we provide an even
higher-order expansion than is available in Hahn and Newey (2004).} This expansion is similar in style to the asymptotic expansion for the cross-sectional model developed earlier, except that the asymptotic order of the leading terms in the expansion are in terms of $T^{-1/2}$ due to the presence of the individual fixed effects. Here the asymptotic bias is
given by $\sqrt{\frac{n}{T}}\mathbf{B} = \sqrt{\frac{n}{T}} E[A_2]$.

The `leave-one-out' jackknife estimator is
\[
\widehat{\theta}_{J}
=T\widehat{\theta}-(T-1)\frac{1}{T}\sum_{t=1}^{T}\widehat{\theta}_{(t)},
\]
where $\widehat{\theta}_{(t)}$ is the estimator formed from the subsample that
excludes time period $t$.\footnote{Jackknife estimators that drop $k$ time periods, rather than just one, could also be used. However, averaging over all $\binom{T}{k}$ leave-$k$-out estimates would be computationally demanding, and so we do not pursue this idea here.}
Dhaene and Jochmans (2015) propose
the use of split-panel jackknifes that only make use of subpanels that contain  consecutive time periods. The split-sample jackknife estimator is
\[
\widehat{\theta}_{SS}=2\widehat{\theta}-\bar{\theta}_{SS},
\]
where $\bar{\theta}_{SS}=\frac{1}{2}(\widehat{\theta}_{1}+\widehat
{\theta}_{2})$, with $\widehat{\theta}_{1}$ being the estimate using observations from the first half of
time periods, and $\widehat{\theta}_{2}$ the estimate that uses the second half of
time periods. Other choices of split-sample jackknife are available; however, the results in
Dhaene and Jochmans (2015) show that non-overlapping sub-panels in general
have lower asymptotic variance, and that among the non-overlapping options,
splitting in two leads to the smallest inflation of higher-order bias. Hence,
in this paper we focus on this half sample version, and simply refer to it as
the split-sample jackknife.

Both the jackknife and split-sample corrections have no impact on the first-order term in (\ref{eq:mle_exp})  so that both  estimators have the same asymptotic variance, equal to $\lim_{n\rightarrow\infty}
\operatorname{Var}(A_1)$ $=\lim_{n\rightarrow\infty}\mathcal{I}
_{n}^{-1}$. To compute the higher-order variance for the estimators, we take the variance
of the first three expansions terms, retaining terms up to $O(T^{-1})$,
similar to what was done in Section \ref{sec:mle_crosssec}. The
following proposition establishes the higher-order variance expressions for the two
estimators, from which we can conclude that the higher-order variance of the split-sample
estimator is larger than that of the jackknife.

\begin{proposition}
\label{thm:variance} Let the regularity conditions stated in Appendix \ref{sec:reg-conds-panel} hold. The higher-order variances of the jackknife and split-sample bias-corrected estimators are given by
\[
\operatorname{Var}(\widehat{\theta}_{J}) \approx \operatorname{Var}(A_1)+\frac{1}{T-1}\tilde{\upsilon},\quad \operatorname{Var}(\widehat{\theta}_{SS})\approx \operatorname{Var}(A_1)
+\frac{2}{T}\tilde{\upsilon},
\]
where
\begin{align*}
\operatorname{Var}(A_1)  &  =\mathcal{I}_{n}^{-1},\\
\tilde{\upsilon}  &  =\mathcal{I}_{n}^{-2}\frac{1}{n}\sum_{i}
\frac{\frac{1}{2}E[U_{it}^{\alpha\alpha}]^{2}+2E[U_{it}^{\alpha\alpha}
]E[V_{it}U_{it}^{\alpha}]
+E[V_{it}^{2}]E[(U_{it}^{\alpha})^{2}]+E[V_{it}U_{it}^{\alpha}]^{2}}{E[V_{it}^{2}
]^{2}}.
\end{align*}

\end{proposition}

\begin{proof}
See Supplementary Appendix I \ref{proof-for-panel-models}.
\end{proof}

\subsection{Accuracy in estimating the bias}

In Section \ref{sec:SS_var}, it was demonstrated that while the jackknife uses a consistent estimate of the bias term, the split-sample jackknife uses an unbiased, but inconsistent estimate. A similar situation arises in the panel model. For the `leave-one-out' jackknife, the implicit bias estimate is equal to $\frac{1}{T}
\widehat{b}_{J}=(T-1)\left(  \frac{1}{T}\sum_{t=1}^{T}\widehat{\theta}_{(t)}
-\widehat{\theta}\right) $, while the split-sample jackknife uses the bias
estimate $\frac{1}{T}\widehat{b}_{SS}=(\bar{\theta}_{SS}-\widehat{\theta})$.
The following
proposition establishes the accuracy of $\widehat{b}_{J}$ and $\widehat{b}_{SS}$
as estimators for $\mathbf{B}=\lim_{n\rightarrow\infty}E[A_2]$.

\begin{proposition}
\label{Thm:bias_est} Let the regularity conditions stated in Appendix \ref{sec:reg-conds-panel} hold. Let $\frac{1}{T}\widehat{b}_{J}=(T-1)(\frac{1}{T}\sum_{t=1}^{T}
\widehat{\theta}_{(t)}-\widehat{\theta})$ and $\frac{1}{T}\widehat{b}_{SS}=(\bar{\theta}
_{SS}-\widehat{\theta})$ be the jackknife and split-sample estimators for the
bias term $\mathbf{B}$. Then,
\begin{align*}
\sqrt{nT}\frac{1}{T}(\widehat{b}_{J} -\mathbf{B})=O_{p}(T^{-1}) \\
\sqrt{nT}\frac{1}{T}(\widehat{b}_{SS}-\mathbf{B})=O_{p}(T^{-1/2}).
\end{align*}
\end{proposition}

\begin{proof}
See Supplementary Appendix  I \ref{proof-for-panel-models}.
\end{proof}

Similarly to the cross-sectional analysis, the jackknife estimates the bias term at a faster rate than the split-sample estimator.
In contrast, in the panel setting the split-sample bias estimate is in fact consistent; it is an average over $n$ unbiased, but inconsistent, estimates of individual-level bias terms.
Nonetheless, Proposition \ref{Thm:bias_est} shows that the jackknife bias
correction affects the third-order, $O_{p}(T^{-1})$, part of the expansion,
while the split-sample bias correction appears as a second-order,
$O_{p}(T^{-1/2})$, term. This implies that the jackknife bias estimate will
only impact the higher-order variance through its covariance with the
first-order term $A_1$, i.e.\ through the term $\operatorname{Cov}\left(
\sqrt{n}A_{1},\sqrt{nT}(\widehat{b}_{J}-\mathbf{B})\right)$. In contrast,
the split-sample bias estimate appears in the higher-order variance both
through its covariance with $A_1$, as well as through its own variance, $\operatorname{Var}\left(\sqrt{n}(\widehat{b}_{SS}-\mathbf{B})\right)$.


It should be noted that the results in Propositions \ref{thm:variance} and \ref{Thm:bias_est} are derived under i.i.d sampling over both individuals and time. In this setting, the results mirror those for the cross sectional analysis. In settings where serial correlation exists, the leave-one-out style jackknife does not provide a valid estimate of the bias, while the  split-sample estimator continues to produce a consistent estimate of the bias, albeit at a slower rate. This is an important advantage of  the method in Dhaene and Jochmans (2015). Although a general result on higher-order properties under serial
correlation is not yet available, it may be reasonable to speculate that the
convergence rate of the estimator of the bias will continue to play an important role. This would suggest that analytical bias corrections involving parametric estimators, when available, are likely to be superior to nonparametric estimators that require estimates of long-run variances.

\subsection{Example 2: panel data\label{sec:example}}

A simple example may help to highlight the results. Suppose that $z_{it}\sim
N\left(  \alpha_{i},\theta\right)  $ are independent across $i$ and $t$. This model was studied by Neyman and Scott (1948). The MLE is given by $\widehat{\theta}=\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}(z_{it}
-\bar{z}_{i})^{2}$, and a standard calculation gives that $E[\widehat{\theta
}]=\frac{T-1}{T}\theta$ so that the MLE has a bias of $\frac{1}{T}
\mathbf{B}=-\frac{1}{T}\theta$. It can be shown that the jackknife
bias-corrected estimator has the form
\[
\tilde{\theta}_{J}=\frac{1}{n(T-1)}\sum_{i=1}^{n}\sum_{t=1}^{T}(z_{it}-\bar
{z}_{i})^{2}
\]
while the split-sample estimator is
\[
\tilde{\theta}_{SS}=2\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}(z_{it}-\bar
{z}_{i})^{2}-\frac{1}{2}\left(  \frac{1}{nM}\sum_{i=1}^{n}\sum_{t=1}
^{M}(z_{it}-\bar{z}_{i,1})^{2}+\frac{1}{nM}\sum_{i=1}^{n}\sum_{t=M+1}
^{T}(z_{it}-\bar{z}_{i,2})^{2}\right)
\]
where $T=2M$, and $\bar{z}_{i,1}$ and $\bar{z}_{i,2}$ are the sample means in
the first and second halves of the sample time period. In this simple model,
the formula given in Proposition \ref{thm:variance} gives $V_{J}
=\frac{2T}{T-1}\theta^{2}$ and $V_{SS}=\frac{2T+4}{T}\theta^{2}$, which we
can easily confirm are also the exact finite sample variances of
the two estimators in this case. Considering the estimation of the bias
itself, we can see that the jackknife bias-correction estimates the bias as
$\frac{1}{T}\widehat{b}_{J}=-\frac{1}{T}\frac{1}{n(T-1)}\sum_{i=1}^{n}
\sum_{t=1}^{T}(z_{it}-\bar{z}_{i})^{2}$, while the split-sample correction
estimates the bias using $\frac{1}{T}\widehat{b}_{SS}=-\frac{1}{2n}\sum
_{i=1}^{n}  \left(  (\bar{z}_{i,1}-\bar{z}_{i})^{2}+(\bar{z}_{i,2}
-\bar{z}_{i})^{2}\right)   $. Both are unbiased estimators of the bias
term in this case, i.e.\ $E[\widehat{b}_{J}]=E[\widehat{b}_{SS}]=-\theta$. It
is straightforward to show $\operatorname{Var}\left(  \sqrt{\frac{n}{T}}
(\widehat{b}_{J}-\mathbf{B})\right)  =\frac{2\theta^{2}}{T(T-1)}$, whereas
$\operatorname{Var}\left(  \sqrt{\frac{n}{T}}(\widehat{b}_{SS}-\mathbf{B}
)\right)  =\frac{2\theta^{2}}{T}$, so that the variance of the split-sample bias
estimate is larger by a factor $T$, as predicted by Proposition \ref{Thm:bias_est}.









\subsection{Extension to time series data}

In Proposition \ref{higher_var}, we noted that the higher-order variance of the bias
corrected MLE can be expressed in terms of $\Upsilon$ in (\ref{eq:higher-var-expr}). This
characterization is useful because it is intuitively positive, which is not
obvious from the definition in (\ref{eq:higher-var}). A natural question is
whether an analytical bias correction in a general time series environment
would lead to a higher-order variance of the form $E\left[  A_{1}^{2}\right]
+\frac{1}{n}\left(  \Upsilon+o\left(  n^{-1}\right)  \right)  $, where
$\Upsilon$ is appropriately replaced by the long-run
counterpart of
\[
E\left[  \left(  \frac{1}{\sqrt{T}}\sum_{t=1}^{T}X_{t}\right)  ^{2}\right]
E\left[  \left(  \frac{1}{\sqrt{T}}\sum_{t=1}^{T}Y_{t}\right)  ^{2}\right]
+\left(  E\left[  \left(  \frac{1}{\sqrt{T}}\sum_{t=1}^{T}X_{t}\right)
\left(  \frac{1}{\sqrt{T}}\sum_{t=1}^{T}Y_{t}\right)  \right]  \right)
^{2}.
\]
Although we are unable to answer the question in its general form, it can be
shown to be valid for a strictly stationary AR(1) model with normally
distributed errors. To be more precise, we go through the formal expansion of
the MLE of the AR(1) model $y_{t}=\theta y_{t-1}+\varepsilon_{t}$, where
$\varepsilon_{t}\sim N\left(  0,\sigma^{2}\right)  $ and $y_{0}\sim N\left(
0,\left.  \sigma^{2}\right/  \left(  1-\theta^{2}\right)  \right)  $.\footnote{The result is available in Supplementary Appendix III.} The
higher-order bias based on the asymptotic expansion for the MLE
yields
$\lim_{T\rightarrow\infty}E\left[  A_{2}\right]  =-2\theta$, so the analytical
bias-corrected estimator takes the form $\widehat{\theta}+\left.
2\widehat{\theta}\right/  T$. This means that the higher-order variance of
the bias-corrected estimator is $\operatorname{Var}\left(  A_{1}\right)
+\frac{1}{T}\operatorname{Var}\left(  A_{2}\right)  +\frac{2}{\sqrt{T}
}E\left[  A_{1}A_{2}\right]  +\frac{2}{T}E\left[  A_{1}\left(  A_{3}
+2A_{1}\right)  \right]  $, which is not obviously positive. We show that this
higher-order variance is equal to $\operatorname{Var}\left(  A_{1}\right)
+\frac{1}{T}\left(  \operatorname{Var}\left(  \mathcal{X}\right)
\operatorname{Var}\left(  \mathcal{Y}\right)  +\left(  E\left[
\mathcal{X}\mathcal{Y}\right]  \right)  ^{2}\right)  $ for
$\mathcal{X}=\left(  1-\theta^{2}\right)  \frac{1}{\sqrt{T}}\sum_{t=1}
^{T}\frac{y_{t-1}\varepsilon_{t}}{\sigma^{2}}$ and $\mathcal{Y}=-\left(
1-\theta^{2}\right)  \frac{1}{\sqrt{T}}\sum_{t=1}^{T}\left(  \frac{y_{t-1}
^{2}}{\sigma^{2}}-\frac{1}{1-\theta^{2}}\right)$. Analogously to the
result in Proposition \ref{higher_var_ss}, we show that the split-sample
jackknife bias-corrected estimator for the AR(1) model has higher-order variance
that is given by $\operatorname{Var}\left(  A_{1}\right)
+\frac{2}{T}\left(  \operatorname{Var}\left(  \mathcal{X}\right)
\operatorname{Var}\left(  \mathcal{Y}\right)  +\left(  E\left[
\mathcal{X}\mathcal{Y}\right]  \right)  ^{2}\right)$, and so
the higher-order part of its variance is larger than that of the
analytical bias-corrected estimator by a factor of two. See
Propositions \ref{prop-AR1-higher-var} and \ref{prop-AR1-higher-var-SS} in Supplementary Appendix III.
These results are
admittedly specific to the AR(1) model. We conjecture that they carry over to
more general parametric time series models as long as the bias has a closed
form parametric expression that can be estimated at parametic rates. On the
other hand, estimators with non-parametrically estimated bias corrections that
typically involve estimated long run variances may not share the same
efficiency properties. We conjecture the same is true for non-parametric
block bootstrap procedures and other subsampling techniques used to estimate
the higher-order bias.









\section{Monte Carlo analysis\label{section-MC}}

To highlight the relevance of the results in a more practical setting,
we conduct two Monte Carlo exercises. In the cross-sectional setting we estimate a marginal treatment effect (MTE) model as in Heckman and Vytlacil (2005). The simulation design follows that used in Catteneo et.\ al.\ (2019) and is a simplified model with no covariates. The treatment $T_i$ is assigned according to $T_i = 1\{ P(Z_i) \geq V_i \}$, where $P(Z_i)$ is a propensity score (which is a function of observed instrumental variables $Z_i$), and $V_i \sim \text{Uniform}[0,1]$ is an unobservable shock that is correlated with potential outcomes and generates selection. Potential outcomes under treated and control states are generated according to
\begin{align*}
    Y_i(0) &= U_{0i}, \quad U_{0i} \vert Z_i, V_i \sim \text{Uniform}[-1,1],\\
    Y_i(1) &= 0.5 + U_{1i}, \quad U_{1i} \vert Z_i, V_i  \sim \text{Uniform}[-0.5,1.5-2V_i].
\end{align*}
The observed outcome variable is given by $Y_i = T_i Y_i(1) + (1-T_i)Y_i(0)$. We create a set of 19 potential instruments, $Z_{j,i}\sim \text{Uniform}[0,1]$ for $j=1,\dots,19$, to be used in estimation, although the true propensity score only depends on the first four of these (in addition to the constant) $P(Z_i)  = 0.1 + Z_{1,i}+ Z_{2,i}+ Z_{3,i}+ Z_{4,i}$.

The MTE function is defined as $MTE(v) = E[Y_1 - Y_0 \vert V=v]$, and measures the average treatment effect for individuals with a given level of unobserved resistance to treatment $V_i=v$. Many objects of interest can be represented as weighted averages of the MTE function, although for the purposes of this simulation we follow Catteneo et.\ al.\ (2019) and focus on $\theta=MTE(0.5)$. The MTE function can be identified as $MTE(v)=\partial E[Y_i \vert P(Z_i)=p]/\partial p \vert_{p=v}$.

To estimate the marginal treatment effect parameter, the propensity score $P(z)$ is first estimated via regression of the treatment dummy variable $T_i$ on $Z_i(k)=(1, Z_{1,i},\dots,Z_{k-1,i})$, where $k$ ranges from 5 to 20. The second step then regresses the outcome $Y_i$ on a quadratic in the propensity score, $Y_i = \beta_1 + \beta_2\widehat{P}(z) + \beta_3\widehat{P}(z)^2$. The scalar parameter of interest $\theta$ is the derivative of this function at $p=0.5$, i.e. $\widehat{\theta} = \widehat{\beta}_2+\widehat{\beta}_3$.

\begin{table}[!h]
\caption{Simulation of marginal treatment effect estimates}
\label{tab:sims_mte}
\centering \small
\begin{threeparttable}

\begin{tabular}{cccccccc}
\hline
 &  & \multicolumn{3}{c}{Conventional} & \multicolumn{3}{c}{Bootstrap} \tabularnewline
n & k & Bias & SD & RMSE & Bias & SD & RMSE \tabularnewline
\hline
1000 & 5  & 0.441 & 4.922 & 4.941 & 0.099 & 5.135 & 5.135 \tabularnewline
     & 10 & 1.059 & 4.712 & 4.828 & 0.184 & 5.381 & 5.383 \tabularnewline
     & 15 & 1.570 & 4.460 & 4.727 & 0.326 & 5.447 & 5.456 \tabularnewline
     & 20 & 1.952 & 4.231 & 4.659 & 0.425 & 5.435 & 5.450 \tabularnewline
\hline
2000 & 5  & 0.317 & 4.730 & 4.740 & 0.074 & 4.830 & 4.830 \tabularnewline
     & 10 & 0.798 & 4.636 & 4.703 & 0.122 & 4.964 & 4.964 \tabularnewline
     & 15 & 1.226 & 4.505 & 4.668 & 0.191 & 5.021 & 5.024 \tabularnewline
     & 20 & 1.604 & 4.421 & 4.702 & 0.264 & 5.105 & 5.110 \tabularnewline
\hline
3000 & 5  & 0.155 & 4.756 & 4.757 & -0.047 & 4.822 & 4.821 \tabularnewline
     & 10 & 0.582 & 4.666 & 4.701 & 0.008 & 4.881 & 4.880 \tabularnewline
     & 15 & 0.971 & 4.574 & 4.675 & 0.069 & 4.924 & 4.923 \tabularnewline
     & 20 & 1.300 & 4.517 & 4.700 & 0.103 & 4.995 & 4.994 \tabularnewline
\hline
&  & \multicolumn{3}{c}{Jackknife} & \multicolumn{3}{c}{Split-Sample} \tabularnewline
n & k & Bias & SD & RMSE & Bias & SD & RMSE \tabularnewline
\hline
1000 & 5  & 0.105 & 5.146 & 5.146 & 0.092 & 5.426 & 5.426 \tabularnewline
     & 10 & 0.254 & 5.320 & 5.325 & 0.243 & 5.683 & 5.686 \tabularnewline
     & 15 & 0.511 & 5.274 & 5.297 & 0.459 & 5.611 & 5.628 \tabularnewline
     & 20 & 0.713 & 5.167 & 5.215 & 0.716 & 5.563 & 5.607 \tabularnewline
\hline
2000 & 5  & 0.075 & 4.840 & 4.839 & 0.087 & 5.012 & 5.012 \tabularnewline
     & 10 & 0.154 & 4.953 & 4.954 & 0.175 & 5.185 & 5.187 \tabularnewline
     & 15 & 0.275 & 4.972 & 4.979 & 0.260 & 5.204 & 5.210 \tabularnewline
     & 20 & 0.423 & 5.014 & 5.031 & 0.395 & 5.267 & 5.281 \tabularnewline
\hline
3000 & 5  & -0.044 & 4.827 & 4.826 & -0.068 & 4.924 & 4.923 \tabularnewline
     & 10 & 0.025 & 4.879 & 4.877 & -0.024 & 5.006 & 5.005 \tabularnewline
     & 15 & 0.122 & 4.899 & 4.899 & 0.068 & 5.093 & 5.092 \tabularnewline
     & 20 & 0.202 & 4.951 & 4.953 & 0.243 & 5.180 & 5.184 \tabularnewline
\hline
\end{tabular}

\begin{tablenotes} \footnotesize
\item[*] Results of estimators over 2000 simulations. \emph{Conventional} denotes the standard two-step estimator described in the text; \emph{Jackknife} denotes the leave-one-out jackknife bias-corrected estimator; \emph{Bootstrap} denotes the bootstrap bias-corrected estimator based on $n/2$ bootstrap draws.
\end{tablenotes}
\end{threeparttable}
\end{table}

We investigate the performance of the bootstrap, jackknife, and split-sample bias corrections. As predicted by the theory, the standard deviation of the bootstrap and jackknife bias-corrected estimators are very similar across the different simulation settings. This is particularly true as the sample size grows large; for $n=3000$ the jackknife and bootstrap bias corrections are very close in terms of standard deviation. The split-sample estimator has larger standard deviation than the other two bias corrections, as expected given its larger higher-order variance. As the sample size grows, the difference decreases; this is expected given that the estimators all have the same asymptotic variance.

As a panel data example, we estimate a probit model with strictly exogenous covariates and individual fixed effects.
\begin{align*}
y_{it}  &  = 1\{\theta_{0}'x_{it}+\alpha_{i}+\varepsilon_{it}>0\},\quad
\varepsilon_{it}\sim N(0,1)
\end{align*}
The simulation is calibrated to the female labor force participation application of Fernandez-Val (2009), and is the same as that used in Fernandez-Val and Weidner (2018). Here the outcome is an indicator for participation in the labor force, and the covariates include three measures of fertility, the number of children aged 0-2, 3-5, and 6-17 years, as well as the log of husband's income, and a quadratic in age. We focus on the coefficients on the three fertility variables. Below we report the results of simulations drawn from a sample of $n=500$ individuals and $T=\{4,8\}$ time periods.\footnote{Fernandez-Val and Weidner (2018) report results using the $n=664$ and $T=9$, which matches the sample size in the PSID data set, and find similar results.}

Table \ref{tab:sims_beta} reports the results for the biased MLE fixed effects estimator as well as three bias corrections: the analytical correction, leave-one-out jackknife, and the split-sample jackknife. We report the bias, standard deviation and root mean-squared error as a percentage of the true coefficient values, and the rejection rate for a test with 5\% significance. The MLE fixed effects estimator has large bias and rejection rates as large as 54\%. As is evident from the theory, the size of the bias is decreasing with the number of time periods. Both the jackknife and analytical bias corrections lead to significant reductions in the bias with no cost in precision; in fact, both bias corrected estimators have smaller standard deviations than the MLE. In contrast, while the split-sample jackknife also reduces bias (although to a lesser degree than the other corrections), it has substantially larger variance and mean-squared error. It is evident that, even for $T=8$, the impact of higher-order differences in the bias corrections remains important for the finite sample properties of the estimator.

\begin{table}[t]
\caption{Simulation of probit model with individual fixed effects}
\label{tab:sims_beta}
\centering \small
\begin{threeparttable}
\begin{tabularx}{0.93\textwidth}{@{}
                       l l S[table-format=4.1,table-column-width=1cm]
                         S[table-format=3.1,table-column-width=1cm]
                         S[table-format=3.1,table-column-width=1cm]
                         S[table-format=1.2,table-column-width=1cm] l
                         S[table-format=3.1,table-column-width=1cm]
                         S[table-format=3.1,table-column-width=1cm]
                         S[table-format=3.1,table-column-width=1cm]
                         S[table-format=1.2,table-column-width=1cm] @{}}
\toprule
&  & \multicolumn{4}{c}{MLE} & & \multicolumn{4}{c}{Analytical}\\
& & {Bias} & {SD} & {RMSE} & {Rej 5\%} & & {Bias} & {SD} & {RMSE} & {Rej 5\%} \\
\hline
T=4 & Ages 0-2 & -41.9 & 24.7 & 48.6 & 0.54 & & -8.2 & 19.0 & 20.7 & 0.06\\
& Ages 3-5 & -42.4 & 47.8 & 63.9 & 0.24 & & -8.2 & 37.4 & 38.3 & 0.05\\
& Ages 6-17 & -42.7 & 132.1 & 138.8 & 0.12 & & -1.1 & 102.9 & 102.9 & 0.03\\ \\[-0.8em]

T=8 & Ages 0-2 & -16.9 & 11.2 & 20.3 & 0.36 & & -3.9 & 10.0 & 10.7 & 0.05\\
& Ages 3-5 & -16.8 & 18.6 & 25.1 & 0.18 & & -3.7 & 16.6 & 17.0 & 0.05\\
& Ages 6-17 & -18.6 & 50.7 & 54.0 & 0.09 & & -4.8 & 45.2 & 45.5 & 0.05\\
\midrule
&  & \multicolumn{4}{c}{Jackknife} & & \multicolumn{4}{c}{Split-sample}\\
&  & {Bias} & {SD} & {RMSE} & {Rej 5\%} & & {Bias} & {SD} & {RMSE} & {Rej 5\%} \\
\hline
T=4 & Ages 0-2 & 21.2 & 16.7 & 27.0 & 0.14 & & 37.7 & 49.5 & 62.2 & 0.53\\
& Ages 3-5 & 20.9 & 30.6 & 37.1 & 0.03 & & 35.2 & 99.5 & 105.5 & 0.46\\
& Ages 6-17 & 22.4 & 88.8 & 91.6 & 0.02 & & 31.2 & 303.7 & 305.3 & 0.44\\  \\[-0.8em]
T=8 & Ages 0-2 & 4.3 & 9.3 & 10.2 & 0.04 & & 7.1 & 16.3 & 17.8 & 0.25\\
& Ages 3-5 & 4.4 & 15.6 & 16.2 & 0.04 & & 7.1 & 30.6 & 31.4 & 0.27\\
& Ages 6-17 & 3.1 & 42.2 & 42.3 & 0.03 & & 8.6 & 82.1 & 82.5 & 0.28\\
\bottomrule
\end{tabularx}

\begin{tablenotes} \footnotesize
\item[*] Results of estimators over 1000 simulations. Bias, SD, and RMSE are percentages of the true parameter values.
\end{tablenotes}
\end{threeparttable}
\end{table}


\section{Summary\label{section-summary}}

We show that the choice of bias correction method does not affect the higher-order
variance of any parametric estimator that is semiparametrically efficient, as long as the bias
estimator is asymptotically linear, i.e. asymptotically equivalent to a sample average. We give a formal expansion
showing this property in a general semiparametric model. We also prove that the bootstrap, jackknife and a version of the analytical bias estimates are asymptotically linear, when the estimator of the parameters of
interest has a standard form of stochastic expansion (which is known to exist for a large class
of models). The result implies that a researcher
may choose a bias correction for an efficient estimator based on computational convenience, or
some other criteria, without affecting its higher-order efficiency.

We have verified this result using derivations of the asymptotic expansion and higher-order variance for maximum likelihood
estimation of a parametric model, using a bootstrap, jackknife, and three forms of analytical bias corrections. Furthermore, we found that the third-order stochastic expansion of the bootstrap, jackknife, and one type of analytical bias-corrected MLE are identical, and hence have an even stronger higher-order equivalence property.

These results show that the higher-order efficiency of bias-corrected efficient estimators does
not depend on the form of the bias correction, as long as the estimate of the bias term is asymptotically linear (and hence $\sqrt{n}$-consistent). Thus, in practice one might use
whatever bias correction method is most convenient. An important caveat
is that the split-sample jackknife estimator does not estimate the bias at the $\sqrt{n}$-rate, and so is not asymptotically linear, and we show the resultant higher-order variance to be
strictly larger in an i.i.d.\ setting, suggesting the importance of the accuracy in estimating the bias.

We generalized the result to the analysis of a panel data model with fixed effects, and established that
the split-sample bias-corrected estimator has larger higher-order variance
than the jackknife estimator, confirming the importance of the accuracy in
estimation of the bias even in panel settings. In non-i.i.d.\ settings, the
standard jackknife cannot be used. Comparison of the split-sample correction
with alternatives, such as the analytical correction given in Hahn and
Kuersteiner (2002, 2011), is a topic that we leave for future research.