# Robust Inference in Panel Data Models:
Some Effects of Heteroskedasticity and Leveraged Data in Small Samples

Annalivia Polselli Institute of Analytics and Data Science (IADS) and Centre
for Micro-Social Change (MiSoc), University of Essex. Email:
<EMAIL_ADDRESS>This is a preliminary revised version of the
first chapter of my doctoral thesis, available at
https://repository.essex.ac.uk/34481/. This project was supported by the South
East Network for Social Sciences (SeNSS) as part of the Doctoral Training
Partnership (award ES/P00072X/1/2128216).

###### Abstract

With the violation of the assumption of homoskedasticity, least squares
estimators of the variance become inefficient and statistical inference
conducted with invalid standard errors leads to misleading rejection rates.
Despite a vast cross-sectional literature on the downward bias of robust
standard errors, the problem is not extensively covered in the panel data
framework. We investigate the consequences of the simultaneous presence of
small sample size, heteroskedasticity and data points that exhibit extreme
values in the covariates (‘good leverage points’) on the statistical
inference. Focusing on one-way linear panel data models, we examine asymptotic
and finite sample properties of a battery of heteroskedasticity-consistent
estimators using Monte Carlo simulations. We also propose a hybrid estimator
of the variance-covariance matrix. Results show that conventional standard
errors are always dominated by more conservative estimators of the variance,
especially in small samples. In addition, all types of HC standard errors have
excellent performances in terms of size and power tests under
homoskedasticity.

JEL codes: C13, C15, C23.
Keywords: cluster-robust standard errors, jackknife methods, test size, power
of test.

##  1 Introduction

When the assumption of homoskedasticity is violated and the disturbances show
non-constant variance (within the cross-sectional timention, or time
dimension, or both), least squares (LS) estimators are no longer efficient.
Consequently, standard errors based on the incorrect assumption of
homoskedastic disturbances lead to misleading statistical inferences. A common
practice is to account for heteroskedasticity with robust standard errors when
estimating the model. The Eicker-Huber-White (EHW) estimator (Eicker,, 1967;
Huber et al.,, 1967; White,, 1980) has become the norm to account for any
degree of heteroskedasticity in the cross-sectional environment. Its
counterpart for the panel data is the Arellano,’s (1987) formula. The presence
of data points that exhibit extreme values in the covariates – i.e., _good
leverage points_ – makes the EHW estimator systematically downward biased
leading to liberal statistical inferences (Long and Ervin,, 2000; Godfrey,,
2006; Hayes and Cai,, 2007; MacKinnon,, 2013; Şimşek and Orhan,, 2016). The
bias is severe when the cross-sectional sample size is sufficiently small
(e.g., with less than 250 units in the sample), and persists even in large
samples (MacKinnon and White,, 1985; Chesher and Jewitt,, 1987; Silva,, 2001;
Verardi and Croux,, 2009). While much discussion has involved the cross-
sectional framework, little has been investigated for panel data, despite
similar issues of Arellano,’s (1987) standard errors.111To the best of our
knowledge, there are only two available studies for panel data. Kezdi, (2003)
compares the finite sample properties of a series of estimators of the
variance-covariance matrix with an without serial correlation in the error
term in large-N and small-T panels. Hansen, (2007) derive the asymptotic
properties of the conventional estimator of the variance-covariance matrix and
studies its finite sample behaviour under heteroskedasticity in the cases
where both $N,T$ jointly go to infinity, and where either $N$ or $T$ goes to
infinity holding the other dimension fixed. Extensions of a class of HC-based
estimators to linear panel data mode ls has been conducted by Cattaneo et al.,
(2018) in high dimensional literature..

In this paper, we investigate the consequences of the simultaneous presence of
small sample size, good leveraged data, and heteroskedastic disturbances on
the validity of the statistical inference in linear panel data models. We
formalise panel versions of MacKinnon and White,’s (1985) and Davidson et
al.,’s (1993) estimators, and propose a new hybrid estimator, $PHC6$, that
penalises only units with high leverage in the covariates. We derive the
asymptotic distributions of this battery of estimators, and analyse their
finite sample properties with Monte Carlo (MC) simulations in terms of
proportional bias, rejection probability (or empirical size), root mean
squared error, and adjusted power. The analysis is conducted across different
panel sample sizes and degrees of heteroskedasticity. Units are randomly
contaminated with good leverage points. While we treat homoskedasticity as a
special case, heteroskedasticity is assumed to be a core component of the
correct regression specification.

We find that under heteroskedasticity and with good leveraged data test
statistics obtained with Arellano,’s (1987) standard errors are, as expected,
over-sized, upward biased, and with low power, especially when the panel size
is smaller than 2,500 observations. Test statistics calculated with PHC6
formula mimic the behaviour of those based on jackknife standard errors in
terms of bias, empirical size and adjusted power test, converging to the same
rates as the sample size increases. The panel version of MacKinnon and
White,’s (1985) estimator shows similar patterns but with different
magnitudes. Under homoskedasticity and with good leveraged data, all
estimators have good performances in terms of proportional bias, rejection
probabilities, and adjusted power, suggesting that the heteroskedasticity
correction should be used. A similar result was found in MacKinnon and White,
(1985) and Long and Ervin, (2000) for cross-sectional models who claimed that
jackknife-type standard errors might enhance inference even with small degrees
of heteroskedasticity.

We focus on small sample sizes for a two reasons. First, the cross-sectional
HC literature has extensively discussed the finite sample bias of the EHW
estimator in the presence of leverage points, and we want to document the
behaviour of Arellano,’s (1987) estimator under the same circumstances.
Second, the nature of the research and/or data availability may force the
investigator to deal with a reduced number of observations in the dataset.

Despite the remarkable methodological contribution in the cross-sectional HC
literature, HC-type estimators222HC-type estimators include: HC2 by Horn et
al., (1975), HC3 by MacKinnon and White, (1985), HC$jk$ by Davidson et al.,
(1993), HC4 by Cribari-Neto, (2004), HC5 by Cribari-Neto et al., (2007), and
HC4m by Cribari-Neto and da Silva, (2011). have not found much application in
practice, although by construction they alleviate the effect of leveraged data
being less sensitive to anomalous cases (Hinkley,, 1977). This study
contributes to the HC literature by creating a link between cross-sectional
and panel HC estimators of the sampling variance. We provide the formulae and
derive the distribution of a selected group of variance-covariance estimators
to panel data. We document the downward bias of conventional robust standard
errors under certain circumstances and provide alternative solutions to obtain
more reliable statistical inferences. This study provides simulation evidence
that these estimators outperform the conventional cluster-robust standard
errors under specific circumstances and should be used in linear panel data
models.

The rest of the paper is structured as follows. Section 2 introduces the
static linear panel data model and its assumptions, and the asymptotic
properties of the _within-group_ estimator. In Section 3, we discuss the
estimation of the variance-covariance matrix, formalise HC estimators for
panel data and propose a new estimator. Section 4 shows the MC simulation
design and discusses the simulation results. In Section 5, we examine the
performances of the four estimators in terms of their proportional bias,
empirical size, adjusted power, and mean squared errors. Section 6 concludes.

##  2 The Model and Estimator

###  2.1 Model and Assumptions

Consider the static linear panel regression model with one-way error component

$y_{it}=\mathbf{x}_{it}^{\prime}\bm{\beta}+\alpha_{i}+u_{it},\hskip
5.0pti\in\mathcal{I}=\\{1,\dots,N\\}\hskip 3.0pt\text{and}\hskip
3.0ptt\in\mathcal{T}=\\{1,\dots,T\\}$ (1)

where $y_{it}$ is the response variable for the cross-sectional unit $i$ at
time period $t$; $\mathbf{x}_{it}$ is a $k\times 1$ vector of time-varying
inputs, $\bm{\beta}$ is a $k\times 1$ vector of parameters of interest;
$\alpha_{i}$ is the individual-specific unobserved heterogeneity (or _fixed
effects_); and $u_{it}$ is a stochastic error component.

Stacking observations for $t$, model (1) at the level of the observation
becomes

$\mathbf{y}_{i}=\mathbf{X}_{i}\bm{\beta}+\bm{\alpha}_{i}+\mathbf{u}_{i},\,\,\,\text{for
all}\,\,\,i=1,\dots,N,$ (2)

where $\mathbf{y}_{i}$ is $T\times 1$ vector of outcomes; $\mathbf{X}_{i}$ is
a $T\times k$ matrix of time-varying regressors;
$\bm{\alpha}_{i}=\alpha_{i}\bm{\iota}$ is a $T\times 1$ vector of individual
fixed effects, and $\bm{\iota}$ is a vector of ones of order $T$; and
$\mathbf{u}_{i}$ is a $T\times 1$ vector of one-way error component. The fixed
effects $\bm{\alpha}_{i}$ in Equation (2) are removed to consistently
estimating the parameter of interest $\bm{\beta}$ by applying an appropriate
transformation of the original data, i.e., the _time-demeaning_ or _first-
differencing_ procedure, because it might be the case that
$\mathbb{E}(\bm{\alpha}_{i}|\mathbf{X}_{i})=h(\mathbf{X}_{i})$. For the rest
of the discussion, we focus on the first approach when applied to Equation
(2). The _time-demeaning_ data transformation delivers a consistent estimator
of $\bm{\beta}$ even when the regressor is correlated with the unobserved
heterogeneity $\alpha_{i}$, but is less efficient than the First-Difference
(FD) transformation with errors that are not identically distributed.

The estimating equation becomes

$\widetilde{\mathbf{y}}_{i}=\widetilde{\mathbf{X}}_{i}\bm{\beta}+\widetilde{\mathbf{u}}_{i},\,\,\,\text{for
all}\,\,\,i=1,\dots,N,$ (3)

where
$\widetilde{\mathbf{y}}_{i}=(\mathbf{I}_{T}-\,T^{-1}\bm{\iota}\bm{\iota}^{\prime})\mathbf{y}_{i}$
is $T\times 1$;
$\widetilde{\mathbf{X}}_{i}=(\mathbf{I}_{T}-\,T^{-1}\bm{\iota}\bm{\iota}^{\prime})\mathbf{X}_{i}$
is $T\times k$; and
$\widetilde{\mathbf{u}}_{i}=(\mathbf{I}_{T}-\,T^{-1}\bm{\iota}\bm{\iota}^{\prime})\mathbf{u}_{i}$
is $T\times 1$. Note that
$(\mathbf{I}_{T}-\,T^{-1}\bm{\iota}\bm{\iota}^{\prime})\bm{\alpha}_{i}=\bm{0}$
as $T^{-1}\bm{\iota}\bm{\iota}^{\prime}\bm{\alpha}_{i}=\bm{\alpha}_{i}$. The
within-group estimator is the Pooled OLS estimator of Equation (3).

The model assumptions are as follows

  1. asm.1

(_data-generating process_):

    1. i

(_independent variables_): $\\{\mathbf{X}_{i}\\}$ is an independent and
identically distributed (_iid_) sequence of random variables, for all
$\,i=1,\dots,N$;

    2. ii

(_disturbances_): $\\{\mathbf{u}_{i}\\}$ is an independent but not identically
distributed (_inid_) sequence of random error terms, for all $\,i=1,\dots,N$.

  2. asm.2

(_on the relation of $\widetilde{\mathbf{X}}_{i}$ and
$\widetilde{\mathbf{u}}_{i}$_):

    1. i

(_strong exogeneity_):
$\mathbb{E}\bigl{(}\widetilde{\mathbf{u}}_{i}|\widetilde{\mathbf{X}}_{i}\bigr{)}=0$,
for all $i=1,\dots,N$;

    2. ii

(_heteroskedasticity_):
$\overline{\bm{\Sigma}}_{N}=N^{-1}\sum_{i=1}^{N}\bm{\Sigma}_{i}\to\bm{\Sigma}$,
where the matrix of the heteroskedastic disturbances
$\bm{\Sigma}_{i}=\mathbb{E}\bigl{(}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}|\widetilde{\mathbf{X}}_{i}\bigr{)}=\mathrm{diag}\bigl{\\{}\sigma_{it}^{2}\bigr{\\}}$
is symmetric of dimension $T$, finite, positive definite, and diagonal.

The above model assumptions have the following implications. asm.1.asm.1i
guarantees that the sequence of random variables
$\\{\mathbf{X}_{i}^{\prime}\mathbf{X}_{i}\\}$ is _iid_ [prop 3.3 in White,
(1984, p.30)]. asm.1.asm.1ii imposes cross-sectional independence and,
together with asm.1.asm.1i, implies that
$\\{\mathbf{X}_{i}^{\prime}\mathbf{u}_{i}\\}$ is an _inid_ sequence of random
vectors [prop 3.10 in White, (1984, p.34)]. Assumption asm.1 and its
implications remain unaltered after any data transformation.

The strict exogeneity assumption asm.2.asm.2i rules out feedback effects and
implies contemporaneous exogeneity, i.e.,
$\mathbb{E}\bigl{(}\widetilde{u}_{it}|\widetilde{\mathbf{X}}_{i}\bigr{)}=\mathbb{E}\bigl{(}\widetilde{u}_{it}|\widetilde{\mathbf{x}}_{it}\bigr{)}=0$,
and is a crucial assumption to prove consistency of the _within-group_
estimator. The projection analog of asm.2.asm.2i is the strong exogeneity
condition, i.e,
$\mathbb{E}\bigl{(}\widetilde{\mathbf{x}}_{is}\widetilde{u}_{it})=0\Leftrightarrow\mathbb{E}\bigl{(}\widetilde{u}_{it}|\widetilde{\mathbf{X}}_{i}\bigr{)}\mbox{=0}$,
for all $s\in\mathcal{T}$ and $s\neq t$. Because the exogeneity of the non-
demeaned variables might not be strong enough to guarantee that exogeneity is
preserved _after_ the transformation333This occurs because the regressor is
correlated with
$T^{-1}\bm{\iota}\bm{\iota}^{\prime}\widetilde{\mathbf{u}}_{i}$ since it
includes the whole history., i.e.,
$\mathbb{E}\bigl{(}\mathbf{X}_{i}^{\prime}\mathbf{u}_{i}\bigr{)}\mbox{=0}\not\Rightarrow\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\bigr{)}\mbox{=0}$
(Cameron and Trivedi,, 2005, p.707). Assumption asm.2.asm.2ii allows the
conditional error variance to vary across observations and time periods, and
imposes serial uncorrelation over time dimension,
$\mathbb{E}\bigl{(}\widetilde{u}_{it}\widetilde{u}_{is}|\widetilde{\mathbf{X}}_{i}\bigr{)}=0$
with $(t,s)\in\mathcal{T}$ and $t\neq s$.

The assumptions for the existence and optimality properties of the estimator
of the true population parameter $\bm{\beta}$ are

  1. asm.3

(_rank condition_): $\mathbf{S}_{N}\equiv
N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
is a finite symmetric matrix with full column rank $k$.

  2. asm.4

(_moment conditions_):

    1. i

$\mathbb{E}\big{\|}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\big{\|}<\infty$
for
$\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\in\mathbb{R}^{k\times
k}$;

    2. ii

$\mathrm{sup}_{i}\,\mathbb{E}\big{\|}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\big{\|}^{2+\delta}<\infty$
for some $\delta>0$, $\forall i$ and
$\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\in\mathbb{R}^{k}$,

where $\|\cdot\|$ denotes the Euclidean norm.

  3. asm.5

(_average variance-covariance matrix convergence_):
$\overline{\mathbf{V}}_{N}=N^{-1}\sum_{i=1}^{N}\,\mathbf{V}_{i}\to\mathbf{V}$,
where
$\mathbf{V}_{i}=\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}\bigr{)}$
and $\mathbf{V}$ is a finite positive definite $k\times k$ matrix.

The full column rank condition in asm.3 implies non-singularity of the matrix
$\mathbf{S}_{N}$ and, hence, no perfect multicollinearity that guarantees the
invertibility of the matrix. The limiting matrix
$\mathbf{S}_{XX}\equiv\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}$
possesses the properties of $\mathbf{S}_{N}$ by the _Weak Law of Large
Numbers_ (_WLLN_) [thm 6.6]. Another implication of asm.3 is that the matrix
of regressors $\widetilde{\mathbf{X}}_{i}$ is full column rank. Assumption
asm.4 defines the finiteness and boundedness of moments in terms of the
Euclidean norm. Assumption asm.5 ensures that the average variance-covariance
matrix converges to a finite quantity, satisfying one of the conditions of the
_Multivariate Central Limit Theorem_ (_MCLT_) for _inid_ processes [thm 6.16
in Hansen, (2019, p.189)].

No restrictions are placed on influential data points – such as, high leverage
points that possess extreme values in the covariates – but we possibly allow
for their presence. We consider a framework where the panel is small, that is,
the time period length is smaller than the number of units $N$ such that $T\ll
N$. Under this notation $T$ is the full set of time information, and the total
number of observations in the sample is given by $n=N\cdot T$ with balanced
data sets.

This set of assumptions and their implications remain valid under any
monotonic data transformation due to the _Continuous Mapping Theorem_ (_CMT_)
[thm 6.19 in Hansen, (2019, p.192)]. Later in this work, we consider the
_within-group_ transformation of the data.

###  2.2 Asymptotic Properties of the Estimator

Under asm.1–asm.4.asm.4i, the _within-group_ estimator of the true population
parameter $\bm{\beta}$ exists with form
$\widehat{\bm{\beta}}_{N}=\Bigl{(}N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\Bigr{)}^{-1}\,\,N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{y}}_{i}$,
and is consistent, i.e.,

$\widehat{\bm{\beta}}_{N}-\bm{\beta}=\Biggl{(}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\Biggr{)}^{-1}\,\,\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\overset{p}{\to}\bm{0}\,\,\,\mathrm{as}\,\,\,N\to\infty.$
(4)

The consistency of the _within-group_ estimator under the aforementioned
assumptions is a known result (as reference, see Hansen,, 2019, pp. 612–613).
By the previously discussed implication of asm.1.asm.1i and prop 3.3 in White,
(1984, p.30),
$\bigl{\\{}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\bigr{\\}}$
is an _iid_ sequence of random variables with finite moments given asm.4. The
elements of the sequence satisfy the _Weak Law of Large Numbers_ (_WLLN_) [thm
6.6 in Hansen, (2019, p.182)] such that
$\mathbf{S}_{N}\overset{p}{\to}\mathbf{S}_{XX}<\infty$. Because both matrices
are invertible by asm.3, then thm 6.19 [_Continuous Mapping Theorem (CMT)_ in
Hansen, (2019, p.192)] yields the result
$\mathbf{S}_{N}^{-1}\overset{p}{\to}\mathbf{S}_{XX}^{-1}$.

Now, we show that the second component in (4) converges in probability to
zero. We know that the sequence
$\bigl{\\{}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\bigr{\\}}$
is _inid_ as an implication of asm.1 [prop 3.10 in White, (1984, p.33)]. Then,
the Chebyshev inequality is

$\mathrm{Pr}\Biggl{(}\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Bigg{\|}\geq\epsilon\Biggr{)}\leq\frac{\mathbb{E}\Big{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Big{\|}^{2}}{\epsilon^{2}},$
(5)

where the numerator in (5) can be expanded as follows

$\begin{split}\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Bigg{\|}^{2}&=\mathrm{tr}\Biggl{\\{}\Biggl{(}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Biggr{)}\,\Biggl{(}\frac{1}{N}\sum_{j=1}^{N}\widetilde{\mathbf{u}}_{j}^{\prime}\widetilde{\mathbf{X}}_{j}\Biggr{)}\Biggr{\\}}\\\
&=\frac{1}{N^{2}}\,\mathrm{tr}\Biggl{\\{}\sum_{i}\sum_{j}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{j}^{\prime}\widetilde{\mathbf{X}}_{j}\Biggr{\\}}.\end{split}$
(6)

By the aforementioned implication of asm.1.asm.1ii and under asm.2.asm.2ii the
conditional error variance is

$\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{j}^{\prime}\widetilde{\mathbf{X}}_{j}\big{|}\widetilde{\mathbf{X}}_{i}\widetilde{\mathbf{X}}_{j}\bigr{)}=\begin{cases}\bm{0}&\forall
i\neq j\\\
\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}&\forall
i=j\end{cases}$ (7)

Applying the expected value operator to (6), and using result (7) jointly with
the _Law of Iterated Expectations_ (LIE), the above equality becomes as
follows

$\begin{split}\mathbb{E}\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Bigg{\|}^{2}&=\frac{1}{N^{2}}\,\mathrm{tr}\Biggl{\\{}\sum_{i}\sum_{j}\,\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{j}^{\prime}\widetilde{\mathbf{X}}_{j}\bigr{)}\Biggr{\\}}\\\
&=\frac{1}{N^{2}}\,\mathrm{tr}\Biggl{\\{}\sum_{i}\sum_{j}\,\mathbb{E}\Bigl{[}\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{j}^{\prime}\widetilde{\mathbf{X}}_{j}\big{|}\widetilde{\mathbf{X}}_{i},\widetilde{\mathbf{X}}_{j}\bigr{)}\Bigr{]}\Biggr{\\}}\\\
&=\frac{1}{N^{2}}\,\mathrm{tr}\Bigl{\\{}\sum_{i}\,\mathbb{E}\Bigl{[}\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}\Bigr{]}\Bigr{\\}}\\\
&=\frac{1}{N}\,\mathrm{tr}\biggl{\\{}\frac{1}{N}\sum_{i}\,\mathbb{E}\Bigl{[}\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}\Bigr{]}\biggr{\\}}\\\
&=\frac{1}{N}\mathrm{tr}\bigl{\\{}\overline{\mathbf{V}}_{N}\bigr{\\}}\to
0,\hskip 5.69054pt\mathrm{as}\,\,N\to\infty\end{split}$ (8)

since assumption (asm.5) implies that
$\mathrm{tr}\bigl{\\{}\overline{\mathbf{V}}_{N}\bigr{\\}}\to\mathrm{tr}\bigl{\\{}\mathbf{V}\bigr{\\}}$,
which is finite.

As a result, the right-hand side of Equality (8) converges in probability to
zero. So does the left-hand side. Inequality (5) becomes

$\mathrm{Pr}\Biggl{(}\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\Bigg{\|}\geq\epsilon\Biggr{)}\to
0\,\,\,\,\mathrm{as}\,\,N\to\infty,$

and, hence,
$N^{-1}\sum_{i}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\overset{p}{\to}\bm{0}$.
By thm 6.19 [CMT in Hansen, (2019, p.192)] , the result follows
$\widehat{\bm{\beta}}_{N}-\bm{\beta}\overset{p}{\to}\mathbf{S}_{XX}^{-1}\cdot\bm{0}=\bm{0}$,
or alternatively $\widehat{\bm{\beta}}_{N}\overset{p}{\to}\bm{\beta}$. This
result holds for any monotonic transformation of the data.

Under asm.1–asm.5, the estimator has the known asymptotic distribution below

$\sqrt{N}\bigl{(}\widehat{\bm{\beta}}_{N}-\bm{\beta}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},\mathbf{S}_{XX}^{-1}\mathbf{V}\mathbf{S}_{XX}^{-1}\bigr{)}\hskip
11.38109pt\mathrm{as}\,\,N\to\infty\text{ and }T\textrm{ fixed}.$ (9)

A reference for this result is Hansen, (2019, pp. 624–625). The left-hand-side
of Equation (9) can be re-written as follows

$\small\sqrt{N}\bigl{(}\widehat{\bm{\beta}}_{N}-\bm{\beta}\bigr{)}=\Biggl{(}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\Biggr{)}^{-1}\,\,\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}.$

The sequence of random variables
$\bigl{\\{}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\bigr{\\}}$
is _iid_ as implication of asm.1.asm.1i and by prop 3.3 in White, (1984,
p.30). With analogous arguments as those used above to prove consistency,
$\mathbf{S}_{N}^{-1}\overset{p}{\to}\mathbf{S}_{XX}^{-1}$. Under assumptions
asm.1 and asm.2.asm.2i and by prop 3.10 in White, (1984, p.33), the sequence
$\bigl{\\{}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\bigr{\\}}\in\mathbb{R}^{k}$
is _inid_ with means
$\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\bigr{)}=0$
and variance matrices
$\mathbf{V}_{i}=\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}\bigr{)}$,
by LIE and asm.2.asm.2ii. The limit in probability asm.5 and assumption
asm.4.asm.4ii are the two conditions that satisfy the _Multivariate Central
Limit Theorem_ (_MCLT_) for _inid_ processes [thm 6.16 in Hansen, (2019,
p.189)]. Therefore,
$N^{-1/2}\,\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},\mathbf{V}\bigr{)}$
as $N\to\infty$. _Slutsky’s Theorem_ [thm 6.22.2 in Hansen, (2019, p.193)]
yields the result
$\sqrt{N}\bigl{(}\widehat{\bm{\beta}}_{N}-\bm{\beta}\bigr{)}\overset{d}{\to}\mathcal{N}\bigl{(}\bm{0},\mathbf{S}_{XX}^{-1}\mathbf{V}\mathbf{S}_{XX}^{-1}\bigr{)}$,
where
$\mathbf{S}_{XX}\equiv\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\bigr{)}$.
thm 6.19 [Hansen, (2019, p.192)] ensures that the above limits hold for any
monotonic transformation of the data, e.g., the within-group transformation.

##  3 Estimating the Asymptotic Variance

Given the above results under the model assumptions we made, the approximate
distribution of the estimator of $\bm{\beta}$ for large but finite samples is

$\widehat{\bm{\beta}}_{N}\mathbin{\overset{a}{\kern
0.0pt\sim}}\mathcal{N}(\bm{\beta},N^{-1}\mathbf{S}_{XX}^{-1}\mathbf{V}\mathbf{S}_{XX}^{-1}),$
(10)

where the limiting matrices $\mathbf{S}_{XX}$ and $\mathbf{V}$ need to be
estimated, and so does the average variance-covariance matrix
$\overline{\mathbf{V}}_{N}$. While $\mathbf{S}_{XX}$ is estimated by
$\mathbf{S}_{N}=N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$,
the estimation of the average variance-covariance matrix needs further
discussion. According to White, (1980), a computationally feasible practice
consists in estimating each expectation,
$\mathbf{V}_{i}=\mathbb{E}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i}\bigr{)}$,
individually, and a plausible estimator of $\overline{\mathbf{V}}_{N}$ would
be
$N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
if the error term were known. Because it is unobserved, a consistent estimator
of the variance-covariance matrix is in practice
$N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\widetilde{\mathbf{u}}}_{i}\widehat{\widetilde{\mathbf{u}}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$,
where
$\widehat{\widetilde{\mathbf{u}}}_{i}=\widetilde{\mathbf{y}}_{i}-\widetilde{\mathbf{X}}_{i}\widehat{\bm{\beta}}$.
Define $\widehat{\widetilde{\mathbf{u}}}_{i}=\widehat{\mathbf{u}}_{i}$ to
simplify the notation.

Using a generalised expression for regression residuals, the variance-
covariance matrix can be re-written as follows:
$\widehat{\mathbf{V}}_{N}=N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$,
where
$\widehat{\mathbf{v}}_{i}=\mathbf{M}_{i}^{-1}\widehat{\widetilde{\mathbf{u}}}_{i}$
are the transformed regression residuals with $\mathbf{M}_{i}$ being the
transformation matrix that differs across estimators of the variance-
covariance. When the transformed residuals equalise the residuals from the
regression, $\widehat{\mathbf{v}}_{i}=\mathbf{I}_{T}\widehat{\mathbf{u}}_{i}$,
the variance-covariance matrix takes the familiar _“sandwich-like”_ formula of
Arellano,’s (1987) estimator.

The variance-covariance matrix $\widehat{\mathbf{V}}_{N}$ with transformed
residuals is still a consistent estimator of the true variance. Let
$\widehat{\bm{\Sigma}}_{i}=\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}$,
from White,’s (1980) general result and under the above model assumptions and
thm 7.7 in Hansen, (2019, p.232), it follows that
$\big{\|}\widehat{\mathbf{V}}_{N}-\overline{\mathbf{V}}_{N}\big{\|}\overset{p}{\to}\bm{0}$
and, hence,
$\big{\|}N^{-1}\sum_{i}\widehat{\bm{\Sigma}}_{i}-\overline{\bm{\Sigma}}_{N}\big{\|}\overset{p}{\to}\bm{0}$,
for all $i=1,\dots,N$, as $N\to\infty$ and keeping $T$ fixed.

In the next sections, we review Arellano,’s (1987) well-known formula,
formalise MacKinnon and White,’s (1985) jackknife-type estimator for panel
data, provide a panel version of Davidson et al.,’s (1993) estimator, and
propose a new hybrid estimator, $PHC6$. The consistency of estimators with
transformed residuals is derived in Appendix D.

###  3.1 HC$k$-type Estimators

The well-known formula of Arellano,’s (1987) estimator (henceforth, PHC0) is

$\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{0}=c_{0}\,\mathbf{S}_{N}^{-1}\widehat{\mathbf{V}}_{N}^{0}\mathbf{S}_{N}^{-1}\,,$
(11)

where $c_{0}=\frac{n-1}{n-k}\cdot\frac{N}{N-1}$, and
$\widehat{\mathbf{V}}_{N}^{0}=N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
with $\mathbf{M}_{i}=\mathbf{I}_{T}$. The finite-sample correction
factor444Computationally, statistical software, like stata, use a finite-
sample modification of the conventional (i.e., Arellano,’s (1987)) variance-
covariance matrix multiplying
$N^{-1}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
by the correction factor $c=\frac{n-1}{n-k}\cdot\frac{N}{N-1}$, where
$n=N\cdot T$ for one-way clustering in panel data, otherwise cluster-robust
standard error turn out to be downward biased (Arellano,, 1987; Bertrand et
al.,, 2004; Cameron et al.,, 2011)., $c_{0}$, ensures that
$\widehat{\mathbf{V}}_{N}^{0}$ is consistent under asm.2.asm.2ii with fixed
$T$; the ratio $N/(N-1)$ is a computational necessary degree-of-freedom
correction to control for individual correlation (Stock and Watson,, 2008;
Cameron et al.,, 2011).

The estimator that resembles Davidson et al.,’s (1993) HC3 in the panel data
framework (PHC3) is as follows

$\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{3}=c_{3}\,\mathbf{S}_{N}^{-1}\widehat{\mathbf{V}}_{N}^{3}\mathbf{S}_{N}^{-1}\,,$
(12)

where $c_{3}=(N-1)N^{-1}$,
$\widehat{\mathbf{V}}_{N}^{3}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
with $\mathbf{M}_{i}=(\mathbf{I}_{T}-\mathbf{H}_{i})$ and the individual
leverage matrix555The individual leverage matrix is a $T\times T$ matrix
defined as follows
$\mathbf{H}_{i}=\begin{pmatrix}h_{i11}&h_{i12}&\dots&h_{i1T}\\\
h_{i21}&h_{i22}&\dots&h_{i21T}\\\ \vdots&\vdots&\vdots&\vdots\\\
h_{iT1}&h_{iT2}&\dots&h_{iTT}\end{pmatrix}\text{ for all }i=1,\dots,N$ (13)
with elements
$h_{its}=\widetilde{\mathbf{x}}_{it}^{\prime}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}})^{-1}\widetilde{\mathbf{x}}_{is}$
with $t,s=1,\dots,T$. ,
$\mathbf{H}_{i}=\widetilde{\mathbf{X}}_{i}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}$,
whose diagonal elements
$h_{itt}=\widetilde{\mathbf{x}}_{it}^{\prime}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}})^{-1}\widetilde{\mathbf{x}}_{it}$
lie in the $(0,1)$ interval but the off-diagonal elements may be negative.
Predicted residuals, $\widehat{\mathbf{v}}_{i}$, assign a penalty to LS
residuals based on the degree of leverage making the estimates of the variance
less sensitive to leverage points. This type of standard errors tend to be
asymptotically conservative as the number of covariates is allowed to grow as
fast as the sample size, despite being asymptotically valid (Cattaneo et al.,,
2018).

The estimator of the jackknife asymptotic variance for panel data models
(PHC$jk$) adapts MacKinnon and White,’s (1985) HC$jk$ estimator and has form

$\displaystyle\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{jk}$
$\displaystyle=\Biggl{(}\frac{N-1}{N}\biggl{)}\sum_{i=1}^{N}\Bigl{(}\widehat{\bm{\beta}}_{(i)}-\bar{\bm{\beta}}\Bigr{)}\Bigl{(}\widehat{\bm{\beta}}_{(i)}-\bar{\bm{\beta}}\Bigr{)}^{\prime}$
$\displaystyle=\Biggl{(}\frac{N-1}{N^{2}}\biggl{)}\mathbf{S}_{N}^{-1}\Bigl{\\{}\widehat{\mathbf{V}}_{N}^{3}-\bm{\mu}^{*}\bm{\mu}^{*}{{}^{\prime}}\Bigl{\\}}\mathbf{S}_{N}^{-1}\,,$
(14)

where the Leave-One-Out estimator is
$\widehat{\bm{\beta}}_{(i)}=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}$
with $\mathbf{M}_{i}=(\mathbf{I}_{T}-\mathbf{H}_{i})$,
$\bar{\bm{\beta}}=\frac{1}{N}\sum_{i=1}^{N}\widehat{\bm{\beta}}_{(i)}$, and
$\bm{\mu}^{*}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}$.
The jackknifed variance-covariance estimator with _fixed effects_ can be found
in Belotti and Peracchi, (2020).

In practice, the jackknife procedure consists in deleting the entire history
of each unit one at a time without replacement. Because the jackknife
resamples in such a way to construct “pseudo-data” on which the estimator of
interest is tested, this technique – as well as the bootstrap – is suitable
for the assessment of the variability of an estimate, e.g., the estimation of
standard errors (Efron,, 1982; Freedman and Peters,, 1984, Chapter 6). The
advantages of the jackknife procedure are double: it is an entirely data-
driven approach, and it is able to alleviate the impact of influential units
on inference (Cattaneo et al.,, 2019). The main drawback is that the jackknife
estimator becomes computationally infeasible for sufficiently large number of
groups.

PHC3 is a special case of Equation (3.1) when the contribution of the second
block is null as $N\to\infty$ and fixing $T$. The two estimators are
asymptotically equivalent and coincide in sufficiently large samples. The
derivation of (3.1) involves considerable algebraic manipulations (see
Appendix B).

###  3.2 A Hybrid Estimator: PHC6

We propose a hybrid estimator of the variance, PHC6, that nests PHC0 and PHC3
estimators using a threshold criterion from the decision rule of the penalty
factor in Cribari-Neto, (2004). PHC3 is chosen because Monte Carlo simulations
showed that Davidson et al.,’s (1993) HC3 possess the best final sample
properties in terms of lower bias, with rejection rates closer to the nominal
one (Long and Ervin,, 2000). The threshold criterion is designed to account
for the time period in which each unit has exerted the maximal leverage with
respect to the average leverage in the same period. PHC6 is designed to
deliver standard errors that are higher in magnitude than PHC0 with
contaminated observations but the same as PHC0 standard errors with no extreme
observations in the sample.

Before presenting the proposed estimator, we clarify beforehand the notation
we will be using. Let the $T\times 1$ vector

$\mathbf{h}_{i}=diag(\mathbf{H}_{i})=\begin{pmatrix}h_{i11}\\\ h_{i22}\\\
\vdots\\\ h_{iTT}\end{pmatrix}\text{ for all }i=1,\dots,N$

be the individual leverage vector constructed from the diagonal elements of
the individual leverage matrix $\mathbf{H}_{i}$ defined in (13), and let the
$T\times 1$ vector

$\mathbf{\overline{h}}_{tt}=\begin{pmatrix}\overline{h}_{11}\\\
\overline{h}_{22}\\\ \vdots\\\
\overline{h}_{TT}\end{pmatrix}=\begin{pmatrix}N^{-1}\sum_{i=i}^{N}h_{i11}\\\
N^{-1}\sum_{i=i}^{N}h_{i22}\\\ \vdots\\\
N^{-1}\sum_{i=i}^{N}h_{iTT}\end{pmatrix}$

be constructed from the average leverage at time $t$ across units. Then, let
$\overline{\mathbf{h}}$ be a $T\times 1$ vector with elements
$\big{(}\mathbf{\overline{h}}_{tt}\exp\circ\mathbf{j}\big{)}$, where the
expression $\exp\circ\mathbf{j}$ indicates the element-wise power of
$\mathbf{j}$ which is a $T\times 1$ vector of negative ones. The Hadamard
(element-wise) product, $\mathbf{h}_{i}\odot\overline{\mathbf{h}}$, is a
$T\times 1$ vector whose elements, $h_{itt}\overline{h}_{tt}^{-1}$, inform on
the relative leverage of unit $i$ at time $t$ with respect to the average
leverage at time $t$. Specifically, values of $h_{itt}\overline{h}_{tt}^{-1}$
above one signal that the relative leverage of unit $i$ at time $t$ exceeds
the average influence at time $t$. Units with values slightly grater than one
cannot automatically be flagged as highly influential because in the absence
of influential units at time $t$, the denominator may be very close to the
numerator, by construction and, hence, one cannot be chosen as cut-off value.
Conversely, high values of $h_{itt}\overline{h}_{tt}^{-1}$ indicate that unit
$i$ is exerting high leverage at time $t$ with respect to the mean influence
at time $t$.

The PHC6 estimator of the variance is defined as follows

$\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{6}=c_{6}\,\mathbf{S}_{N}^{-1}\widehat{\mathbf{V}}_{N}^{6}\mathbf{S}_{N}^{-1}\,,$
(15)

where the variance-covariance matrix is
$\widehat{\mathbf{V}}_{N}^{6}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i},$
and the matrix $\mathbf{M}_{i}$ has functional form

$\mathbf{M}_{i}=\begin{cases}\mathbf{I}_{T}&\text{if }h^{*}_{i}<2\\\
\mathbf{I}_{T}-\mathbf{H}_{i}&\text{otherwise}\end{cases}$ (16)

where
$h^{*}_{i}=max\big{\\{}h_{i11}/\overline{h}_{11},\dots,h_{iTT}/\overline{h}_{TT}\big{\\}}$
is the maximal individual leverage of unit $i$; and
$\overline{h}_{tt}=N^{-1}\sum_{i=i}^{N}h_{itt}$ is the average leverage at
time $t$, with $h_{itt}$ being the individual leverage of unit $i$ at time
$t$. The finite sample correction of PHC6 is

$c_{6}=\begin{cases}\frac{(NT-1)N}{(NT-k)(N-1)}&\text{if }h^{*}_{i}<2\\\
\frac{N-1}{N}&\text{otherwise}\end{cases}\\\ .$

According to the cut-off rule, residuals of units with maximal individual
relative leverage, $h^{*}_{i}=h_{itt}\overline{h}_{tt}^{-1}$, are discounted
by the penalty matrix $\mathbf{M}_{i}$. Unlike PHC3 that penalises both low
and high leverage points at the same rate, PHC6 discounts at the same
discounting rate as PHC3 only if the unit exerts high leverage. When the
individual relative leverage does not exceed the cutoff, no penalty is applied
and PHC6 coincides with Arellano,’s (1987) estimator. Conversely, when the
average level of leverage exceeds the cut-off value, PHC6 residuals are
penalised as in PHC3. In addition, PHC6 always weights for a final sample
correction.

The cut-off is set to be equal to 2 such that no penalty is assigned to fairly
influential units at time $t$. One is not chosen as a cutoff value because in
the absence of anomalous cases, the denominator $\overline{h}_{tt}$ would be
very close to the numerator $h_{itt}$ for some units with meaningless
individual leverage but above the mean average. This would drive the ratio to
exceed one.

##  4 Monte Carlo Simulation

In this section, we present the MC simulation design which illustrates the
behaviour of the four types of estimators of the variance in finite
samples666Monte Carlo simulations provide computational evidence of finite
sample properties of an estimator or a test when applied to fictitious data
(Hendry,, 1984; Kiviet et al.,, 2012)., when variables are contaminated with
anomalous data points. For simplicity, the simulation set up uses synthetic
balanced data set and does not allow for any correlation between the
individual-specific fixed effects and the regressor777This design leaves open
the possibility to estimate the regression equation consistently and
efficiently using the _random effects_ (RE) estimator. However, our objective
is not to analyse RE because its assumptions are unlikely to be satisfied in
practice. Also, we are not focusing on unbalanced datasets, whose discussion
is postponed to future analysis while addressing the issue of attrition in
panel data.. The data generating process for the Monte Carlo simulation is
designed to be closely related to: (i) Godfrey, (2006), Stock and Watson,
(2008), and MacKinnon, (2013) in terms of the form of heteroskedasticity,
number of regressors and the calibrated parameters; and (ii) Bramati and
Croux, (2007) for the contamination with cell-isolated good leverage points.
However, we depart from these settings by making some modifications to the
simulation designs.

The data generating process (DGP) of Monte Carlo simulations is as follows

$\displaystyle
y_{it}=\beta_{0}+\sum_{k=1}^{K}\beta_{j}x_{it,k}+\alpha_{i}+u_{it},\text{ for
all }i=1,\dots,N\text{ and }t=1,\dots,T_{i}$ (17) $\displaystyle
x_{k}\mathbin{\overset{}{\kern 0.0pt\sim}}\mathcal{N}(0,1)\,\,\textrm{for
$k=\\{1,2\\}$ except contaminated cases}$ (18) $\displaystyle
x_{k}=f(x_{1},x_{2})\,\,\textrm{for}\,\,k=\\{3,4,5\\}$ (19)
$\displaystyle\alpha_{i}\mathbin{\overset{}{\kern 0.0pt\sim}}{}U(0,1)$ (20)
$\displaystyle
u_{it}=\sigma_{it}\epsilon_{it}+\theta\epsilon_{it-1},\,\,\epsilon_{it}\mathbin{\overset{}{\kern
0.0pt\sim}}\mathcal{N}\bigl{(}0,1\bigr{)},\,\,u_{it}\mathbin{\overset{}{\kern
0.0pt\sim}}\mathcal{N}\bigl{(}0,\sigma_{it}^{2}\bigr{)}$ (21)
$\displaystyle\sigma_{it}^{2}=z(\gamma)\Biggl{(}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\Biggr{)}^{\gamma},\hskip
5.0pt\mathrm{with}\,\,z(\gamma)=\Biggl{[}\mathbb{E}\Biggl{(}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\Biggr{)}^{\gamma}\Biggr{]}^{-1}$
(22)

where the number of regressors in the model is $K=5$ and $K=J$; model
parameters are calibrated to be $\beta_{k}=1$, for $k=1,\dots,4$, and
$\beta_{5}=0$; $\theta=0$ because errors are conditionally serially
uncorrelated by assumption as in Stock and Watson, (2008); the degree of
heteroskedasticity assumes values of $\gamma=\\{0,2\\}$, where $\gamma=0$
stands for homoskedasticity and $\gamma\gg 1$ for severe heteroskedasticity.
The scaling factor, $z(\gamma)$, is chosen such that the average variance of
the error term is equal to one888The error term $u_{it}$ is intrinsically
heteroskedastic but not on average due to the presence of the scaling factor
$z(\gamma)$. The distribution of the random variable
$\mathrm{W}=\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}$ and $\mathrm{W}^{2}$ is
provided in Appendix F. The algebraic derivation of the means and variances
are shown..

The contamination of random variables with good leverage points is completely
random over the observations (i.e., _cell-isolated_ anomalous cases). Good
leverage points are obtained by randomly replacing 10% of the values999The
degree of contamination could have been set to be even more or less severe
according to the relevance the researcher attributes to the presence of
extreme observations in the sample.  of $x_{1}$ with extreme observations
drawn from a normal distribution with mean $\mu_{x_{1}}=5$ and standard
deviation $\sigma_{x_{1}}=25$. Because $x_{1}$ is contaminated, then the
variables generated from the former are directly affected by this source of
contamination. The remaining random variables – $x_{3},x_{4},x_{5}$ – are
either generated from the square or the product of $x_{1}$ and $x_{2}$ and,
hence, follow a $\chi^{2}_{(\nu_{1})}$ and a Gamma distribution, respectively.

The model is estimated including the set of aforementioned time-varying
covariates and individual specific fixed effects, $\alpha_{i}$. We estimate
model (17) using fixed effects (FE) by applying the _within-group_ (or time-
demeaning) transformation to simulated data. Then, we estimate the time-
demeaned regression specification using OLS101010We do not use the FGLS-FE to
estimate the estimating Equation (17) for three main reasons. First, when the
sample size is not sufficiently large there is an efficiency loss with respect
to the FE-OLS estimator. In this analysis, we are interested in investigating
the finite sample properties of the estimator, when $N$ is not very large.
Second, the FGLS-FE procedure requires to drop one of the time periods because
the variance matrix is not invertible, leading to the reduction of the
(already small) panel sample size (Cameron and Trivedi,, 2005, ch.21.6,
p.729). Third, FGLS-FE relies on the quality of the estimation of the variance
and on the knowledge of the form of heteroskedasticity. However, the form of
heteroskedasticity is always unknown from the data and the researcher has to
make assumptions on the relationship between the variance of the disturbances
and observables and unknown parameters (Cameron and Trivedi,, 2005, Chapter
21, pp. 720-721, 729). This is unpractical in many areas of application and
subjective to the researcher’s guess. To overcome this limit, an objective
criterion that has become a standard practice in applied works consists in
using conventional robust standard errors due to software facilities
(Verbeek,, 2008).. As in Hansen, (2007), the DGP for the simulations involves
only random effects (RE) model because with (20) we assumed that the
unobserved fixed effect is uncorrelated with the regressors. The model could
be estimated more efficiently with RE but FE models are commonly used in
empirical studies with panel data111111In future analysis we will re-assess
the current version of the Monte Carlo simulation allowing $\alpha_{i}$ to be
correlated with $\mathbf{x}_{i}$ to satisfy FE assumptions. Under this
simulation design, $\widehat{\bm{\beta}}$ estimated with FE remains consistent
but is less efficient than $\widehat{\bm{\beta}}$ estimated with RE..

Our Monte Carlo simulation involves 10,000 replications. The simulations are
run for a combination of cross-sectional units $N=\\{25,50,150,500\\}$ and
time periods $T=\\{2,5,10,20\\}$. Both cross-sectional units and time periods
can be grouped as small ($N=\\{25,50\\}$; $T=\\{2,5\\}$), moderately small
($N=\\{150\\}$; $T=\\{10\\}$), and moderately large ($N=\\{500\\}$;
$T=\\{20\\}$). The simulation is programmed in Stata16-mp and the main
procedure is implemented in mata.

##  5 Testing the Performance of HC Estimators

We examine the performance of each estimator in terms of proportional bias
(PB), rejection probability (RP, or empirical size), adjusted power test, and
root mean squared error (RMSE). Results are provided for a battery of
estimators by a combination of panel units, time periods, and degree of
heteroskedasticity, $\\{N,T,\gamma\\}$, where the number of units $N$ varies
in an interval from 25 to 500 units, time is fixed at $T=\\{2,5,10,20\\}$, and
the parameter that controls for the degree of heteroskedasticity is
$\gamma\in\\{0,2\\}$. This design is in accordance with the finite $T$
assumption in the model as time periods are fixed while the number of
observations increases.

Good leveraged data and heteroskedasticity make, as expected, test statistics
calculated with conventional robust standard errors over-sized, upward biased,
and with low power when the panel size $n<2,500$. The proposed PHC6 mimics the
behaviour of PHCjk in terms of PB, RP and power in all samples. PHC3 shows
similar patterns but with different magnitudes.

###  5.1 Rejection Probability and Probability Bias

RP (i.e., the size of a test) in a Monte Carlo exercise with $R$ runs is the
frequency at which a rejection of the true null hypothesis occurs on average.
A test statistics has a good size if rejects the null hypothesis approximately
around the chosen $\alpha\%$ of the simulations, when the model is generated
under the assumption that the null hypothesis is actually true.

The steps to obtain the empirical size in a two-sided single coefficient test
are as follows. First, for each combination of $\\{N,T\\}$ and each simulation
run $r=1,\dots,R$, compute the test statistics under the true null hypothesis,

$T_{N,T}^{0}(\widehat{\bm{\beta}}_{N,T,r})=\frac{\bigl{(}\widehat{\bm{\beta}}_{N,T,r}-\bm{\beta}^{0}\bigr{)}}{\sqrt{\widehat{\mathrm{AVar}(\widehat{\bm{\beta}}_{N,T,r})}}}\mathbin{\overset{a}{\kern
0.0pt\sim}}t_{(df_{r},\,\alpha/2)}.$

Second, set the indicator $\mathds{1}\\{\cdot\\}$ to turn on when the null
hypothesis is rejected according to the rule

$J_{N,T,r}^{0}(\widehat{\bm{\beta}})\equiv\mathds{1}\bigl{\\{}\big{|}T_{N,T}^{0}(\widehat{\bm{\beta}}_{N,T,r})\big{|}>t_{(df_{r},\,\alpha/2)}\bigr{\\}},$

where $t_{(df_{r},\,\alpha/2)}$ is the critical value from a student-t
distribution with $df_{r},$ degrees of freedom for a two-sided hypothesis
test121212With non-clustered inference $df_{r}=(NT-1)-(N+k-1)$ otherwise
$df_{r}=N-1$.. Third, count the total number of times a rejection has occurred
and average it out by the number of replications $R$; the empirical size
denotes the percentage of rejections in the Monte Carlo exercise as

$\bar{J}_{N,T,r}^{0}(\widehat{\bm{\beta}})\equiv\frac{1}{R}\sum_{r=1}^{R}J_{N,T,r}^{0}(\widehat{\bm{\beta}})=\alpha_{test}.$

For a two-sided test with $q$ linear restrictions, the coverage probability is
computed as follows. First, for each combination of $\\{N,T\\}$ and each
simulation run $r=1,\dots,R$, compute the Wald statistics under the true null
hypothesis, $H_{0}:\mathbf{R}\bm{\beta}-\mathbf{r}^{0}=\bm{0}$,

$W_{N,T,r}^{0}(\widehat{\bm{\beta}})=N(\mathbf{R}\widehat{\bm{\beta}}_{N,T,r}-\mathbf{r}^{0})^{\prime}\Bigl{\\{}\mathbf{R}\widehat{\mathrm{AVar}(\widehat{\bm{\beta}}_{N,T,r})}\mathbf{R}^{\prime}\Bigr{\\}}^{-1}(\mathbf{R}\widehat{\bm{\beta}}_{N,T,r}-\mathbf{r}^{0})\mathbin{\overset{a}{\kern
0.0pt\sim}}\chi^{2}(q),$

where $\mathbf{R}$ is a $q\times K$ matrix with $q\leq K$, and
$\mathbf{r}^{0}$ is a $q\times 1$ vector. Second, Mark as one every time a
rejection occurs according to the rule

$\tilde{J}_{N,T,r}^{0}(\widehat{\bm{\beta}})\equiv\mathds{1}\bigl{\\{}W_{N,T,r}^{0}(\widehat{\bm{\beta}})>cv_{\chi^{2}(q)}\bigr{\\}},$

where $cv_{\chi^{2}(q)}$ is the critical value from a $\chi^{2}$ distribution
with $q$ degrees of freedom for a two-sided hypothesis
test131313Alternatively, the F statistic can be computed from the Wald test
statistics as
$F_{N,T,r}^{0}(\widehat{\bm{\beta}})=W_{N,T,r}^{0}(\widehat{\bm{\beta}})/q\mathbin{\overset{a}{\kern
0.0pt\sim}}F_{\alpha}(q,df_{r})$ under the true null hypothesis, where $q$ are
the number of restrictions and degrees of freedom at the numerator, and
$df_{r}$ are the residual degrees of freedom or degrees of freedom at the
denominator.. Third, sum the cases when the null hypothesis has been rejected
according to the above rule, and divide the number by the total number of
simulation runs. The empirical size for a joint coefficient test is given by
the percentage of rejections in the overall Monte Carlo as follows

$\bar{\tilde{J}}^{0}_{N,T,r}(\widehat{\bm{\beta}})\equiv
R^{-1}\sum_{r=1}^{R}\tilde{J}^{0}_{N,T,r}(\widehat{\bm{\beta}})=\alpha_{test}.$

In the simulations, we test $H_{0}:\beta_{j}=1$ against $H_{1}:\beta_{j}\neq
1$ for $j=1$ while in a two-sided joint test we test
$H_{0}:\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=1$ against $H_{1}:\text{at
least one }\beta_{j}\neq 1$, for $j=1,\dots,4$. The closer the rejection
probability is to the nominal level of $5\%$, the better the estimator’s
performance in terms of empirical size (or type I error).

The proportional bias (PB) is a measure of the bias of the estimator of the
variance-covariance matrix computed as
$PB=1-SE(\hat{\beta}_{j})/SD(\hat{\beta}_{j})$, where _SE_ stands for standard
error and _SD_ for standard deviation. Positive (negative) values of PB
indicate by how much the standard error obtained using one of the four
formulae presented above underestimates (overestimates) the “true” standard
error.

In this section, we comment on the performance of each estimator taking into
account its ability to reject the true null hypothesis at 5% significance
level along with its accuracy. Tables 1 and 2 report the results of the Monte
Carlo simulations respectively, with and without heteroskedasticity. Each
table compares the PB, RP and RMSE141414Results for the RMSE are commented in
Section 5.2. of four alternative formulae of the variance-covariance matrix
(i.e., PHC0, PHC3, PHC6 and PHCjk). Results are grouped by different
combinations of sample size $N$ and time length $T$. Figures refer to the
slope parameter $\beta_{1}$, which is associated with the contaminated
variable $x_{it,1}$. The t-test statistics are at 5%-level.

Under heteroskedasticity, PHC0 standard errors considerably underestimate the
“true” variance (positive PB) on average by at least $30\%$ when $n\leq
2,500$. PHC6 mimics the behaviour of PHCjk in small and large samples,
overestimating the true variance (negative PB: min= 1.2% and max = 12.3%) for
$n\leq 300$ and slightly underestimating the true variance (positive PB: min=
4.9% and max = 10.6%) in the other cases. For $N=\\{25,50\\}$ and all $T$ the
PB of PHCjk is larger in absolute value than PB of PHC6 if the bias is
positive, and smaller otherwise. From $N\geq 150$ PHCjk and PHC6 produce the
same bias but PHC3 produces a smaller bias in absolute value when the
estimators over-estimate the variance.

Test statistics of PHC0 are largely over-sized (RP above 0.05) when
$N=\\{25,50\\}$ and all $T$ but approach the true $\alpha\%$-size when $n\geq
5,000$, despite the high positive PB. The most conservative estimators always
under-reject the null hypothesis (RP below 0.05), and as the cross-sectional
size increases (fixing the time dimension) the RP gradually converges to 5%
but their test statistics still remain slightly under-sized. However, looking
at the (positive/negative) distance from 0.05 PHC0 turns out to be more over-
sized then he other estimators when $n\leq 750$.

In general, a smaller PB in absolute value (signaling a good approximation of
the “true” variance) does not automatically imply that the empirical size is
the closest to the actual nominal significance level. The “true” standard
errors remain under-estimated (over-estimated) if the bias is positive,
despite producing test statistics that reject the null hypothesis with much
precision.

Under homoskedasticity, PHC0 always underestimates the true variance
(especially for $n\leq 1,500$). The PB reduces as the panel size increases but
only when $n=10,000$ it drops considerably. The other PHC estimators tend to
over-estimate the true variance (negative PB) but the magnitudes are smaller
in absolute value than the figures of PHC0. PHC6 has similar bias to PHCjk
while PHC3 is slightly more biased. Test-statistics of PHC0 are over-sized
(large RP) for $n\leq 1,500$ but show a convergence pattern to 5% as the
sample size increases. The test size of PHC6 and PHCjk is always closest to
the true $\alpha$-size followed by PHC3.

Tables 3 and 4 report the Wald test statistics and RP from the joint
coefficient test for the slope coefficients different from zero (i.e.,
$\beta_{i}$ for $i=1,\dots,4$) under heteroskedasticity and homoskedasticity,
respectively. Results for different combinations of $\\{N,T\\}$ are displayed.
The nominal level of significance is set at $\alpha=0.05$. The closer the
value of the rejection rate of the test statistic is to $\alpha=0.05$, the
better the estimator’s performance in terms of empirical size.

Under heteroskedasticity and good leveraged data points, the RPs of the four
estimators slowly converge to 5% as the sample size increases with exception
of PHC6. PHC6 is outperformed by PHC0 in terms of RP for $n\geq 2,500$.
Despite the upward distortion of all test statistics for $N\leq 500$ and all
$T$, PHC0 raw sizes are the largest in magnitude among the four estimators in
very small cross-sectional samples. The Wald statistics of the other two
conservative estimators are the lowest in magnitudes for $n\geq 300$. Similar
patterns are observed under the assumption of homoskedasticity. PHC0 performs
as well as the two conservative estimators only for large $N$ ($N\geq 150$
fixing $T$).

###  5.2 RMSE Assessment

An additional evaluation on the quality of the four estimators is done in
terms of the RMSE. For each estimator of the variance, the RMSE is computed as
the square root of the average deviation of the standard error from the
standard deviation of the estimated coefficient of $\beta_{j}$. In formulae,

$\mathrm{RMSE}_{j}^{s}=\frac{1}{R}\sum_{r=1}^{R}\sqrt{(\hat{\sigma}^{s}(\hat{\beta_{j}})_{r}-\sigma(\hat{\beta_{j}})_{r})^{2}}$
(23)

where $\hat{\sigma}^{s}(\hat{\beta_{j}})_{r}$ is the standard error of
$\hat{\beta_{j}}$ in the $r$th run of the simulation computed using one of the
HC formulae, and $\sigma(\hat{\beta_{j}})_{r}$ is the standard deviation of
the estimated coefficient $\beta_{j}$. A good quality estimator has its RMSE
close to zero. Because the RMSE and PB are constructed from the same
quantities, $\hat{\sigma}^{s}(\hat{\beta_{j}})$ and $\sigma(\hat{\beta_{j}})$,
they are linked one to the other. The larger the proportional bias in absolute
value, the larger the RMSE of the estimator is in magnitude.

Results are presented in Table 1 and 2 for different combinations of cross-
sectional units and time length, and under different degrees of
heteroskedacity. Under heteroskedasticity, the RMSE of PHC0 estimator is much
higher than those of the other three estimators for all combinations of $N$
and $T$. The RMSE of the three conservative estimators gradually converges to
zero in large samples, displaying similar values in small samples. Under
homoskedasticity, the RMSE of all estimators are always very close to zero for
different combinations of panel sample size. The only exception is for $n\leq
100$ when PHC6 has the smallest RMSE.

###  5.3 Adjusted Power Test

The power of the test is the average frequency at which the false null
hypothesis is rejected in a simulation. In a two-sided single coefficient
test, the adjusted power for the false null hypothesis is obtained through the
steps below. First, for each combination of $\\{N,T\\}$ and for each
simulation run $r=1,\dots,R$, compute the test statistics under the false null
hypothesis as

$T_{N,T}^{1}(\widehat{\bm{\beta}}_{N,T,r})=\frac{\bigl{(}\widehat{\bm{\beta}}_{N,T,r}-\bm{\beta}^{1}\bigr{)}}{\sqrt{\widehat{\mathrm{AVar}(\widehat{\bm{\beta}}_{N,T,r})}}}\mathbin{\overset{a}{\kern
0.0pt\sim}}t_{(df_{r},\,\alpha/2)}.$

Second, the indicator $\mathds{1}\\{\cdot\\}$ turns on every time that the
rejection rule holds

$J_{N,T,r}^{1}(\widehat{\bm{\beta}})\equiv\mathds{1}\bigl{\\{}T_{N,T}^{1}(\widehat{\bm{\beta}}_{N,T,r})<\textbf{t}_{\alpha/2}^{0}\textrm{
or
}T_{N,T}^{1}(\widehat{\bm{\beta}}_{N,T,r})>\textbf{t}_{1-\alpha/2}^{0}\bigr{\\}},$

where $\textbf{t}_{\alpha/2}^{0}$ and $\textbf{t}_{1-\alpha/2}^{0}$ are values
lying respectively at the $(\alpha/2)^{th}$ and $(1-\alpha/2)^{th}$
percentiles of $T_{N,T}^{0}(\widehat{\bm{\beta}}_{N,T,r})$, and used as
critical values151515We cannot use conventional critical values from the
t-distribution because size-unadjusted power curves make any comparison
between estimators meaningless.. The empirical critical values differ due to
the asymmetric distribution of the test statistics. Third, count the total
number of rejections in the simulation and divide by the number of runs; the
adjusted power of a test is

$\bar{J}^{1}_{N,T,r}(\widehat{\bm{\beta}})\equiv
R^{-1}\sum_{r=1}^{R}J_{N,T,r}^{1}(\widehat{\bm{\beta}})=1-\theta_{test}.$

Similarly, for a two-sided test with $q$ linear restrictions the adjusted
power of a test is conducted as follows. First, for each combination of
$\\{N,T\\}$ and for each simulation run $r=1,\dots,R$, compute the Wald
statistics under the true null hypothesis,
$H_{0}:\mathbf{R}\bm{\beta}-\mathbf{r}^{1}=\bm{0}$,

$W_{N,T,r}^{1}(\widehat{\bm{\beta}})=N(\mathbf{R}\widehat{\bm{\beta}}_{N,T,r}-\mathbf{r}^{1})^{\prime}\Bigl{\\{}\mathbf{R}\widehat{\mathrm{AVar}(\widehat{\bm{\beta}}_{N,T,r})}\mathbf{R}^{\prime}\Bigr{\\}}^{-1}(\mathbf{R}\widehat{\bm{\beta}}_{N,T,r}-\mathbf{r}^{1})\mathbin{\overset{a}{\kern
0.0pt\sim}}\chi^{2}(q),$

where $\mathbf{r}^{1}$ is a $q\times 1$ vector. Second, define the F
statistics
$F_{N,T,r}^{1}(\widehat{\bm{\beta}})=W_{N,T,r}^{1}(\widehat{\bm{\beta}})/q$
under the false null hypothesis for replication run $r$, and sample
combination $\\{N,T\\}$. The rejection rule is defined as

$\tilde{J}^{1}_{N,T,r}(\widehat{\bm{\beta}})\equiv\mathds{1}\bigl{\\{}F_{N,T,r}^{1}(\widehat{\bm{\beta}})>F_{\alpha}^{0}\bigr{\\}},$

where $F^{0}_{\alpha}$ is the value lying at the $\alpha^{th}$ quantile of
distribution of $F_{N,T,r}^{0}(\widehat{\bm{\beta}})$ derived under the true
null hypothesis, and used as empirical critical in the rejection rule. Third,
the percentage of rejections that occur in the Monte Carlo exercise is the
adjusted power of a test,

$\bar{\tilde{J}}^{1}_{N,T,r}(\widehat{\bm{\beta}})\equiv
R^{-1}\sum_{r=1}^{R}\tilde{J}^{1}_{N,T,r}(\widehat{\bm{\beta}})=1-\theta_{test}.$

In the simulations, we test $H_{0}:\beta_{j}=1$ against $H_{1}:\beta_{j}\neq
001$ for $j=\\{1,2\\}$ for two-sided single coefficient tests, where
$\beta^{1}$ is a value taken from a narrow interval around the true
$\beta_{j}$. For two-sided joint tests we test
$H_{0}:\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=1$ against $H_{1}:\text{at
least one }\beta_{j}\neq 1$, for $j=1,\dots,4$.

Figures 1 and 2 plot size-adjusted power curves of a battery of HC estimators
for different panel sample sizes and degree of heteroskedasticity for
$\beta_{1}$. The vertices of all power curves correspond to the nominal size
of the test statistics, $\alpha=0.05$. It is common practice to adjust the
power for the empirical size because the empirical distributions of test
statistics may depend on the nature of the specific regressor and, therefore,
any comparison across estimators turns out to be meaningless without size-
adjustment. Precisely, in the absence of any size-adjustment the most liberal
estimator would tend to have greater power than the most conservative
estimator because the former is more likely to over-reject the null hypothesis
in favour of the alternative, while the opposite is true for the latter.
Unlike the test size, simulation results for the test power do not differ
considerably in terms of the overall pattern, but they do in terms of
magnitudes.

Under heteroskedasticity, simulation results show that PHC0 does not have as
good power performance as PHC3, PHC$jk$ and PHC6 in small samples
($N=\\{25,50\\}$ and especially with $T=2$). In fact, its rejection
probabilities at a given parameter value are lower than those of the other
three estimators. Fixing T and letting N change, the power performance of PHC0
does not improve. Rejection probabilities remain the lowest and slowly
converge to one, even when the distance from the true value of $\beta$
increases. Conversely, we do not observe such a remarkable loss in power when
we let $T$ increase and fixing $N$ as the difference with other estimators in
the rejection probabilities at a given parameter value becomes negligible or
vanishes completely.

Under the assumption of homoskedasticity, PHC0 has better power than PHC3,
PHC$jk$, and PHC6 with $N=25$ for all $T$. This result is in stark contrast
with PHC0 poor test size (i.e., RP) described above due to the usual trade-off
between type I and type II error. When $T=2$, all power curves show a lower
convergence to one.

Figures 3 and 4 show the adjusted power curves for the joint coefficient test.
From the graphs we observe that all power curves are well-behaved under
homoskedasticity with rejection rates approaching one quite rapidly as the
tested parameter values depart from the true value, and with the increase in
the sample size. This cannot be said under heteroskedasticity and, especially,
when the panel sample size is small (small $N$ and small $T$) because test
statistics of all estimators have low rejection power, especially PHC0 test
statistics when $N=\\{25,50\\}$ and $T=\\{2,5\\}$.

Overall, the four estimators have similar asymptotic behaviour with or without
heteroskedasticity. This can be explained by the sensitivity of the test of
hypothesis to sample size. In fact, as the sample size increases the
probability of rejecting the false null hypothesis (i.e., the power of the
test) increases as well, by construction. The opposite happens to the size of
a test instead.

##  6 Conclusion

In this chapter, we investigated the effects of the simultaneous presence of a
small sample size, heteroskedasticity, and good leveraged data on the validity
of conventional statistical inference in linear panel data models with fixed
effects. We documented their detrimental effects on the statistical inference
calculated with robust standard errors. More conservative estimators of the
sampling variance produce test statistics that have unbiased empirical sizes
and higher power under these circumstance.

We formalised a panel version of MacKinnon and White,’s (1985) and Davidson et
al.,’s (1993) estimators, and proposed a new hybrid estimator, $PHC6$. We
derived the finite sample properties and the asymptotic distributions of the
discussed HC estimators. With MC simulations we compared the performances of
four types of standard errors, computed with Arellano,’s (1987) and three
types of jackknife-like formulae, in terms of empirical size and power. We
documented the downward bias of conventional robust standard errors under
specific circumstances, suggesting alternatives to obtain more reliable
statistical inference.

The main findings can be summarised as follows. Under heteroskedasticity, more
conservative standard errors should be used in the presence of leverage points
because their test statistics possess a low proportional bias, small size
distortions, and have higher power. Conversely, conventional standard errors
and the proposed formula, PHC6, should be preferred with homoskedasticity
because the other conservative estimators excessively under-reject the true
null hypothesis. Under homoskedasticity cluster-robust formulae should always
be used. A similar result was found in MacKinnon and White, (1985) and Long
and Ervin, (2000) for cross-sectional models. The cross-sectional dimension
matters for the finite sample properties of the estimators but not the size of
$N$ relative to $T$. However, conventional cluster-robust standard errors
remain upward biased even when their empirical size is correct, and even in
larger samples.

## References

  * Arellano, (1987) Arellano, M. (1987).  Practitioners’ corner: Computing robust standard errors for within-groups estimators.  Oxford bulletin of Economics and Statistics, 49(4):431–434.
  * Banerjee and Frees, (1997) Banerjee, M. and Frees, E. W. (1997).  Influence diagnostics for linear longitudinal models.  Journal of the American Statistical Association, 92(439):999–1005.
  * Belotti and Peracchi, (2020) Belotti, F. and Peracchi, F. (2020).  Fast leave-one-out methods for inference, model selection, and diagnostic checking.  The Stata Journal, 20(4):785–804.
  * Bertrand et al., (2004) Bertrand, M., Duflo, E., and Mullainathan, S. (2004).  How much should we trust differences-in-differences estimates?  The Quarterly journal of economics, 119(1):249–275.
  * Bramati and Croux, (2007) Bramati, M. C. and Croux, C. (2007).  Robust estimators for the fixed effects panel data model.  The econometrics journal, 10(3):521–540.
  * Cameron et al., (2011) Cameron, A. C., Gelbach, J. B., and Miller, D. L. (2011).  Robust inference with multiway clustering.  Journal of Business & Economic Statistics, 29(2):238–249.
  * Cameron and Trivedi, (2005) Cameron, A. C. and Trivedi, P. K. (2005).  Microeconometrics: methods and applications.  Cambridge university press.
  * Cattaneo et al., (2019) Cattaneo, M. D., Jansson, M., and Ma, X. (2019).  Two-step estimation and inference with possibly many included covariates.  The Review of Economic Studies, 86(3):1095–1122.
  * Cattaneo et al., (2018) Cattaneo, M. D., Jansson, M., and Newey, W. K. (2018).  Inference in linear regression models with many covariates and heteroscedasticity.  Journal of the American Statistical Association, 113(523):1350–1361.
  * Chesher and Jewitt, (1987) Chesher, A. and Jewitt, I. (1987).  The bias of a heteroskedasticity consistent covariance matrix estimator.  Econometrica: Journal of the Econometric Society, pages 1217–1222.
  * Cribari-Neto, (2004) Cribari-Neto, F. (2004).  Asymptotic inference under heteroskedasticity of unknown form.  Computational Statistics & Data Analysis, 45(2):215–233.
  * Cribari-Neto and da Silva, (2011) Cribari-Neto, F. and da Silva, W. B. (2011).  A new heteroskedasticity-consistent covariance matrix estimator for the linear regression model.  AStA Advances in Statistical Analysis, 95(2):129–146.
  * Cribari-Neto et al., (2007) Cribari-Neto, F., Souza, T. C., and Vasconcellos, K. L. (2007).  Inference under heteroskedasticity and leveraged data.  Communications in Statistics - Theory and Methods, 36(10):1877–1888.
  * Davidson et al., (1993) Davidson, R., MacKinnon, J. G., et al. (1993).  Estimation and inference in econometrics.  OUP Catalogue.
  * Efron, (1982) Efron, B. (1982).  The jackknife, the bootstrap, and other resampling plans, volume 38.  Siam.
  * Eicker, (1967) Eicker, F. (1967).  Limit theorems for regressions with unequal and dependent errors.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages 59–82.
  * Freedman and Peters, (1984) Freedman, D. A. and Peters, S. C. (1984).  Bootstrapping an econometric model: Some empirical results.  Journal of Business & Economic Statistics, 2(2):150–158.
  * Godfrey, (2006) Godfrey, L. (2006).  Tests for regression models with heteroskedasticity of unknown form.  Computational Statistics & Data Analysis, 50(10):2715–2733.
  * Hansen, (2019) Hansen, B. E. (2019).  Econometrics.  Unpublished manuscript. Latest version: February 2019.
  * Hansen, (2007) Hansen, C. B. (2007).  Asymptotic properties of a robust variance matrix estimator for panel data when t is large.  Journal of Econometrics, 141(2):597–620.
  * Hayes and Cai, (2007) Hayes, A. F. and Cai, L. (2007).  Using heteroskedasticity-consistent standard error estimators in ols regression: An introduction and software implementation.  Behavior research methods, 39(4):709–722.
  * Hendry, (1984) Hendry, D. F. (1984).  Monte carlo experimentation in econometrics.  Handbook of econometrics, 2:937–976.
  * Hinkley, (1977) Hinkley, D. V. (1977).  Jackknifing in unbalanced situations.  Technometrics, 19(3):285–292.
  * Horn et al., (1975) Horn, S. D., Horn, R. A., and Duncan, D. B. (1975).  Estimating heteroscedastic variances in linear models.  Journal of the American Statistical Association, 70(350):380–385.
  * Huber et al., (1967) Huber, P. J. et al. (1967).  The behavior of maximum likelihood estimates under nonstandard conditions.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1, pages 221–233. University of California Press.
  * Kezdi, (2003) Kezdi, G. (2003).  Robust standard error estimation in fixed-effects panel models.  Available at SSRN 596988.
  * Kiviet et al., (2012) Kiviet, J. F. et al. (2012).  Monte Carlo simulation for econometricians.  now publishers.
  * Long and Ervin, (2000) Long, J. S. and Ervin, L. H. (2000).  Using heteroscedasticity consistent standard errors in the linear regression model.  The American Statistician, 54(3):217–224.
  * MacKinnon, (2013) MacKinnon, J. G. (2013).  Thirty years of heteroskedasticity-robust inference.  In Recent advances and future directions in causality, prediction, and specification analysis, pages 437–461. Springer.
  * MacKinnon and White, (1985) MacKinnon, J. G. and White, H. (1985).  Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties.  Journal of econometrics, 29(3):305–325.
  * Silva, (2001) Silva, J. S. (2001).  Influence diagnostics and estimation algorithms for powell’s scls.  Journal of Business & Economic Statistics, 19(1):55–62.
  * Şimşek and Orhan, (2016) Şimşek, E. and Orhan, M. (2016).  Heteroskedasticity-consistent covariance matrix estimators in small samples with high leverage points.  Theoretical Economics Letters, 6(04):658.
  * Stock and Watson, (2008) Stock, J. H. and Watson, M. W. (2008).  Heteroskedasticity-robust standard errors for fixed effects panel data regression.  Econometrica, 76(1):155–174.
  * Verardi and Croux, (2009) Verardi, V. and Croux, C. (2009).  Robust regression in stata.  The Stata Journal, 9(3):439–453.
  * Verbeek, (2008) Verbeek, M. (2008).  A guide to modern econometrics.  John Wiley & Sons.
  * White, (1980) White, H. (1980).  A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  Econometrica: Journal of the Econometric Society, pages 817–838.
  * White, (1984) White, H. (1984).  Asymptotic theory for Econometricians.  Academic press.

##  Appendix A Tables and Figures

Table 1: Single hypothesis test, heteroskedasticity

| Heteroskedasticity ($\gamma=2$)
---|---
| PB | RP | RMSE | PB | RP | RMSE | PB | RP | RMSE | PB | RP | RMSE
$(N,T)$ | (25, 2) | (50, 2) | (150, 2) | (500, 2)
PHC0 | 0.713 | 0.516 | 0.305 | 0.625 | 0.407 | 0.263 | 0.450 | 0.223 | 0.147 | 0.316 | 0.098 | 0.081
PHC3 | -0.083 | 0.017 | 0.035 | -0.138 | 0.023 | 0.058 | -0.014 | 0.029 | 0.005 | 0.084 | 0.027 | 0.021
PHC6 | -0.042 | 0.020 | 0.018 | -0.123 | 0.025 | 0.052 | -0.010 | 0.030 | 0.003 | 0.085 | 0.027 | 0.022
PHCjk | -0.039 | 0.018 | 0.017 | -0.119 | 0.024 | 0.050 | -0.010 | 0.030 | 0.003 | 0.085 | 0.027 | 0.022
| (25, 5) | (50, 5) | (150, 5) | (500, 5)
PHC0 | 0.578 | 0.337 | 0.224 | 0.473 | 0.238 | 0.160 | 0.338 | 0.112 | 0.089 | 0.225 | 0.062 | 0.042
PHC3 | -0.121 | 0.022 | 0.047 | -0.024 | 0.027 | 0.008 | 0.066 | 0.026 | 0.017 | 0.099 | 0.033 | 0.018
PHC6 | -0.098 | 0.023 | 0.038 | -0.024 | 0.027 | 0.008 | 0.069 | 0.026 | 0.018 | 0.100 | 0.033 | 0.018
PHCjk | -0.086 | 0.023 | 0.033 | -0.010 | 0.029 | 0.003 | 0.069 | 0.026 | 0.018 | 0.100 | 0.033 | 0.018
| (25, 10) | (50, 10) | (150, 10) | (500, 10)
PHC0 | 0.472 | 0.219 | 0.024 | 0.395 | 0.154 | 0.117 | 0.277 | 0.079 | 0.062 | 0.184 | 0.052 | 0.027
PHC3 | -0.033 | 0.022 | 0.011 | 0.052 | 0.027 | 0.015 | 0.098 | 0.029 | 0.022 | 0.105 | 0.035 | 0.015
PHC6 | -0.012 | 0.025 | 0.004 | 0.061 | 0.029 | 0.018 | 0.101 | 0.030 | 0.023 | 0.106 | 0.035 | 0.015
PHCjk | -0.006 | 0.024 | 0.002 | 0.063 | 0.029 | 0.019 | 0.101 | 0.030 | 0.023 | 0.106 | 0.035 | 0.015
| (25, 20) | (50, 20) | (150, 20) | (500, 20)
PHC0 | 0.385 | 0.138 | 0.112 | 0.308 | 0.088 | 0.075 | 0.211 | 0.056 | 0.037 | 0.131 | 0.052 | 0.014
PHC3 | 0.029 | 0.021 | 0.008 | 0.076 | 0.024 | 0.019 | 0.096 | 0.031 | 0.017 | 0.084 | 0.042 | 0.009
PHC6 | 0.049 | 0.026 | 0.014 | 0.085 | 0.027 | 0.021 | 0.099 | 0.031 | 0.017 | 0.086 | 0.042 | 0.009
PHCjk | 0.052 | 0.025 | 0.015 | 0.086 | 0.026 | 0.021 | 0.099 | 0.031 | 0.017 | 0.085 | 0.042 | 0.009

  * _The number of replications is 10,000. The random variable associated with slope parameter $\beta_{1}$ is contaminated with leverage points and drives heteroskedasticity. PB: Proportional Bias. Positive values indicate by how much the standard error underestimates the “true” standard error. RP: Rejection Probability of 5%-level t-test on $\beta_{1}$ (i.e., size of test). RMSE: Root Mean Squared Error. _

Table 2: Single hypothesis test, homoskedasticity

| Homoskedasticity ($\gamma=0$)
---|---
| PB | RP | RMSE | PB | RP | RMSE | PB | RP | RMSE | PB | RP | RMSE
$(N,T)$ | (25, 2) | (50, 2) | (150, 2) | (500, 2)
PHC0 | 0.369 | 0.204 | 0.043 | 0.348 | 0.192 | 0.014 | 0.174 | 0.116 | 0.002 | 0.058 | 0.067 | 0.000
PHC3 | -0.411 | 0.028 | 0.048 | -0.374 | 0.040 | 0.015 | -0.140 | 0.049 | 0.002 | -0.049 | 0.046 | 0.000
PHC6 | -0.328 | 0.036 | 0.038 | -0.346 | 0.044 | 0.014 | -0.131 | 0.050 | 0.002 | -0.046 | 0.046 | 0.000
PHCjk | -0.361 | 0.030 | 0.042 | -0.353 | 0.041 | 0.014 | -0.135 | 0.049 | 0.002 | -0.048 | 0.046 | 0.000
| (25, 5) | (50, 5) | (150, 5) | (500, 5)
PHC0 | 0.324 | 0.160 | 0.008 | 0.205 | 0.118 | 0.002 | 0.074 | 0.075 | 0.000 | 0.017 | 0.055 | 0.000
PHC3 | -0.310 | 0.040 | 0.007 | -0.160 | 0.050 | 0.002 | -0.062 | 0.047 | 0.000 | -0.027 | 0.046 | 0.000
PHC6 | -0.284 | 0.045 | 0.007 | -0.148 | 0.052 | 0.002 | -0.059 | 0.048 | 0.000 | -0.026 | 0.047 | 0.000
PHCjk | -0.273 | 0.043 | 0.006 | -0.146 | 0.052 | 0.002 | -0.059 | 0.048 | 0.000 | -0.026 | 0.047 | 0.000
| (25, 10) | (50, 10) | (150, 10) | (500, 10)
PHC0 | 0.197 | 0.113 | 0.002 | 0.125 | 0.095 | 0.001 | 0.041 | 0.066 | 0.000 | 0.025 | 0.057 | 0.000
PHC3 | -0.179 | 0.043 | 0.002 | -0.071 | 0.050 | 0.000 | -0.032 | 0.050 | 0.000 | 0.002 | 0.051 | 0.000
PHC6 | -0.156 | 0.048 | 0.002 | -0.060 | 0.053 | 0.000 | -0.028 | 0.051 | 0.000 | 0.003 | 0.051 | 0.000
PHCjk | -0.151 | 0.047 | 0.002 | -0.059 | 0.052 | 0.000 | -0.028 | 0.051 | 0.000 | 0.003 | 0.051 | 0.000
| (25, 20) | (50, 20) | (150, 20) | (500, 20)
PHC0 | 0.114 | 0.084 | 0.001 | 0.065 | 0.068 | 0.001 | 0.019 | 0.053 | 0.000 | 0.008 | 0.052 | 0.000
PHC3 | -0.097 | 0.047 | 0.001 | -0.043 | 0.048 | 0.001 | -0.020 | 0.045 | 0.000 | -0.004 | 0.049 | 0.000
PHC6 | -0.076 | 0.053 | 0.001 | -0.033 | 0.050 | 0.000 | -0.016 | 0.046 | 0.000 | -0.003 | 0.049 | 0.000
PHCjk | -0.074 | 0.050 | 0.000 | -0.033 | 0.050 | 0.000 | -0.016 | 0.046 | 0.000 | -0.003 | 0.049 | 0.000

  * _The number of replications is 10,000. The random variable associated with slope parameter $\beta_{1}$ is contaminated with leverage points and drives heteroskedasticity. PB: Proportional Bias. Positive values indicate by how much the standard error underestimates the “true” standard error. RP: Rejection Probability of 5%-level t-test on $\beta_{1}$ (i.e., size of test). RMSE: Root Mean Squared Error. _

Table 3: Joint hypothesis test, heteroskedasticity

| Heteroskedasticity ($\gamma=2$)
---|---
| Wald Stats | RP | Wald Stats | RP | Wald Stats | RP | Wald Stats | RP
$(N,T)$ | (25,2) | (50, 2) | (150, 2) | (500, 2)
PHC0 | 3341260.250 | 0.923 | 29613.070 | 0.820 | 6.124 | 0.510 | 1.775 | 0.210
PHC3 | 117.439 | 0.112 | 9.203 | 0.150 | 1.513 | 0.129 | 1.105 | 0.076
PHC6 | 66.965 | 0.264 | 9.947 | 0.281 | 2.495 | 0.223 | 1.606 | 0.189
PHCjk | 124.333 | 0.115 | 9.505 | 0.154 | 1.525 | 0.130 | 1.107 | 0.077
| (25, 5) | (50, 5) | (150, 5) | (500, 5)
PHC0 | 1293.329 | 0.768 | 19.074 | 0.569 | 2.138 | 0.263 | 1.365 | 0.115
PHC3 | 9.351 | 0.169 | 1.964 | 0.140 | 1.175 | 0.093 | 1.081 | 0.060
PHC6 | 10.219 | 0.303 | 2.872 | 0.248 | 1.737 | 0.197 | 1.709 | 0.217
PHCjk | 10.377 | 0.178 | 2.019 | 0.145 | 1.184 | 0.094 | 1.084 | 0.060
| (25, 10) | (50, 10) | (150, 10) | (500, 10)
PHC0 | 16.202 | 0.568 | 3.348 | 0.355 | 1.540 | 0.153 | 1.210 | 0.076
PHC3 | 2.157 | 0.152 | 1.353 | 0.110 | 1.086 | 0.065 | 1.053 | 0.044
PHC6 | 3.565 | 0.256 | 2.107 | 0.213 | 1.639 | 0.191 | 1.864 | 0.264
PHCjk | 2.317 | 0.163 | 1.389 | 0.113 | 1.094 | 0.067 | 1.055 | 0.045
| (25, 20) | (50, 20) | (150, 20) | (500, 20)
PHC0 | 4.169 | 0.375 | 1.972 | 0.213 | 1.328 | 0.098 | 1.145 | 0.064
PHC3 | 1.568 | 0.119 | 1.201 | 0.082 | 1.077 | 0.052 | 1.055 | 0.046
PHC6 | 2.438 | 0.234 | 1.820 | 0.202 | 1.773 | 0.230 | 2.311 | 0.390
PHCjk | 1.658 | 0.130 | 1.228 | 0.087 | 1.085 | 0.054 | 1.057 | 0.046

  * _The number of replications is 10,000. Tested hypothesis $H_{0}:\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=1$. Random variables associated with slope parameters $\beta_{1}$ and $\beta_{3}$ are contaminated with leverage points. All random variables drive heteroskedasticity. RP: Rejection Probability of 5%-level t-test (i.e., size of test). _

Table 4: Joint hypothesis test, homoskedasticity

| Homoskedasticity ($\gamma=0$)
---|---
| Wald Stats | RP | Wald Stats | RP | Wald Stats | RP | Wald Stats | RP
$(N,T)$ | (25,2) | (50, 2) | (150, 2) | (500, 2)
PHC0 | 118.278 | 0.724 | 12.974 | 0.575 | 2.298 | 0.298 | 1.341 | 0.136
PHC3 | 1.444 | 0.103 | 1.478 | 0.132 | 1.235 | 0.112 | 1.098 | 0.082
PHC6 | 24.939 | 0.323 | 10.828 | 0.427 | 5.985 | 0.471 | 4.638 | 0.495
PHCjk | 1.521 | 0.109 | 1.514 | 0.136 | 1.244 | 0.113 | 1.100 | 0.083
| (25, 5) | (50, 5) | (150, 5) | (500, 5)
PHC0 | 8.066 | 0.522 | 2.845 | 0.340 | 1.463 | 0.162 | 1.144 | 0.087
PHC3 | 1.855 | 0.155 | 1.399 | 0.132 | 1.134 | 0.086 | 1.050 | 0.066
PHC6 | 10.577 | 0.470 | 7.094 | 0.484 | 5.191 | 0.491 | 5.164 | 0.661
PHCjk | 1.957 | 0.166 | 1.431 | 0.136 | 1.142 | 0.088 | 1.053 | 0.066
| (25, 10) | (50, 10) | (150, 10) | (500, 10)
PHC0 | 3.252 | 0.343 | 1.853 | 0.216 | 1.265 | 0.110 | 1.075 | 0.069
PHC3 | 1.588 | 0.139 | 1.286 | 0.110 | 1.102 | 0.078 | 1.026 | 0.058
PHC6 | 8.413 | 0.492 | 5.917 | 0.476 | 4.834 | 0.545 | 6.700 | 0.845
PHCjk | 1.665 | 0.147 | 1.313 | 0.115 | 1.109 | 0.079 | 1.028 | 0.059
| (25, 20) | (50, 20) | (150, 20) | (500, 20)
PHC0 | 2.104 | 0.227 | 1.465 | 0.143 | 1.152 | 0.079 | 1.055 | 0.066
PHC3 | 1.437 | 0.121 | 1.196 | 0.088 | 1.067 | 0.060 | 1.029 | 0.059
PHC6 | 6.842 | 0.234 | 5.376 | 0.508 | 5.584 | 0.705 | 10.257 | 0.986
PHCjk | 1.501 | 0.489 | 1.221 | 0.094 | 1.075 | 0.062 | 1.031 | 0.060

  * _The number of replications is 10,000. Tested hypothesis $H_{0}:\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=1$. Random variables associated with slope parameters $\beta_{1}$ and $\beta_{3}$ are contaminated with leverage points. All random variables drive heteroskedasticity. RP: Rejection Probability of 5%-level t-test (i.e., size of test). _

Figure 1: Power test for $\beta_{1}$, heteroskedasticity

(a)

(b)

(c)

(d)

Figure 2: Power test for $\beta_{1}$, homoskedasticity

(a)

(b)

(c)

(d)

Figure 3: Power test for joint coefficient test, heteroskedasticity

(a)

(b)

(c)

(d)

Figure 4: Power test for joint coefficient test, homoskedasticity

(a)

(b)

(c)

(d)

##  Appendix B Leave-One-Out (L1O) Estimator

Following the L1O estimator for RE models in Banerjee and Frees, (1997), we
derive $\widehat{\bm{\beta}}_{(i)}$ using Woodbury’s formula
$(\mathrm{A}+\mathrm{B}\mathrm{D}\mathrm{C})^{-1}=\mathrm{A}^{-1}-\mathrm{A}^{-1}\mathrm{B}\,(\mathrm{D}^{-1}+\mathrm{C}\mathrm{A}^{-1}\mathrm{B})^{-1}\mathrm{C}\mathrm{A}^{-1}$,
where $\mathrm{A}=\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}$,
$\mathrm{B}=-\widetilde{\mathbf{X}}_{i}^{\prime}$,
$\mathrm{C}=\widetilde{\mathbf{X}}_{i}$, and $\mathrm{D}=\mathbf{I}_{T}$.

$\displaystyle\widehat{\bm{\beta}}_{(i)}$
$\displaystyle=\biggl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\biggr{)}^{-1}\biggl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{Y}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{y}}_{i}\biggr{)}$
$\displaystyle=\Biggr{(}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}+\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\biggl{(}\underbrace{\mathbf{I}_{T}-\widetilde{\mathbf{X}}_{i}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}}_{=\,\mathbf{I}_{T}-\,\mathbf{H}_{i}}\biggr{)}^{-1}\widetilde{\mathbf{X}}_{i}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\Biggr{)}\times\biggl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{Y}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{y}}_{i}\biggr{)}$
$\displaystyle=\underbrace{\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{Y}}}_{=\,\widehat{\bm{\beta}}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{y}}_{i}+\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\widetilde{\mathbf{X}}_{i}\underbrace{\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{Y}}}_{=\,\widehat{\bm{\beta}}}$
$\displaystyle-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\underbrace{\widetilde{\mathbf{X}}_{i}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}}_{=\,\mathbf{H}_{i}}\widetilde{\mathbf{y}}_{i}$
$\displaystyle=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\bigl{[}(\mathbf{I}_{T}-\mathbf{H}_{i})\widetilde{\mathbf{y}}_{i}-\widetilde{\mathbf{X}}_{i}\widehat{\bm{\beta}}+\mathbf{H}_{i}\widetilde{\mathbf{y}}_{i}\bigr{]}$
$\displaystyle=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}(\widetilde{\mathbf{y}}_{i}-\widetilde{\mathbf{X}}_{i}\widehat{\bm{\beta}})$
$\displaystyle=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\widehat{\mathbf{u}}_{i}$
$\displaystyle=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i},$
(24)

where $\mathbf{M}_{i}^{-1}=(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}$. Result (24)
is the L1O estimator for FE in Belotti and Peracchi, (2020). The sample mean
of (24) is

$\displaystyle\bar{\bm{\beta}}\equiv$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\widehat{\bm{\beta}}_{(i)}=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\bm{\mu}^{*}\,,$
(25)

where
$\frac{1}{N}\sum_{i=1}^{N}\sum_{t=1}^{T}\,\widehat{\bm{\beta}}=N\,\widehat{\bm{\beta}}$,
$\bm{\mu}^{*}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}$
is a $k\times 1$ vector. Therefore, from (24) and (25) we get

$\displaystyle\widehat{\bm{\beta}}_{(i)}-\bar{\bm{\beta}}$
$\displaystyle=\widehat{\bm{\beta}}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}-\widehat{\bm{\beta}}+\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\bm{\mu}^{*}$
$\displaystyle=-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}-\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\bm{\mu}^{*}\bigr{)}$
(26)

##  Appendix C Derivation of the Jackknife Estimator

Following the procedure in Hansen, (2019, pp. 324–326), the jackknife
estimator of variance can be computed as

$\displaystyle\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{jk}$
$\displaystyle=\Biggl{(}\frac{N-1}{N}\biggl{)}\sum_{i=1}^{N}\Bigl{(}\widehat{\bm{\beta}}_{(i)}-\bar{\bm{\beta}}\Bigr{)}\Bigl{(}\widehat{\bm{\beta}}_{(i)}-\bar{\bm{\beta}}\Bigr{)}^{\prime}$
(27)
$\displaystyle=\Biggl{(}\frac{N-1}{N}\biggl{)}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\sum_{i=1}^{N}\,\Biggl{\\{}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}-\bm{\mu}^{*}\bigr{)}\bigl{(}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}-\bm{\mu}^{*}\bigr{)}^{\prime}\biggl{\\}}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}$
(28)
$\displaystyle=\Biggl{(}\frac{N-1}{N}\biggl{)}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\Biggl{\\{}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widetilde{\mathbf{X}}_{i}-N\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}\bm{\mu}^{*}{{}^{\prime}}$
$\displaystyle-N\bm{\mu}^{*}\frac{1}{N}\sum_{i=1}^{N}\widehat{\mathbf{u}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widetilde{\mathbf{X}}_{i}+N\bm{\mu}^{*}\bm{\mu}^{*}{{}^{\prime}}\Biggl{\\}}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}$
(29)
$\displaystyle=\Biggl{(}\frac{N-1}{N}\biggl{)}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}\Biggl{\\{}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widetilde{\mathbf{X}}_{i}-N\,\bm{\mu}^{*}\bm{\mu}^{*}{{}^{\prime}}\Biggl{\\}}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\bigr{)}^{-1}$
(30)
$\displaystyle=\Biggl{(}\frac{N-1}{N^{2}}\biggl{)}\Biggl{(}\frac{\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}}{N}\Biggr{)}^{-1}\Bigl{\\{}\widehat{\mathbf{V}}_{N}^{3}-\bm{\mu}^{*}\bm{\mu}^{*}{{}^{\prime}}\Bigl{\\}}\Biggl{(}\frac{\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}}{N}\Biggr{)}^{-1}\,,$
(31)

where
$\bm{\mu}^{*}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\mathbf{M}_{i}^{-1}\widehat{\mathbf{u}}_{i}$.

##  Appendix D Proof Consistency of Transformed Residuals

We show that
$\widehat{\mathbf{v}}_{i}=(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\widehat{\mathbf{u}}_{i}\overset{p}{\to}\widetilde{\mathbf{u}}_{i}$.
We start from

$\displaystyle\widehat{\mathbf{v}}_{i}-\widehat{\mathbf{u}}_{i}$
$\displaystyle=(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}\widehat{\mathbf{u}}_{i}-\widehat{\mathbf{u}}_{i}$
$\displaystyle=\bigl{(}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}-\mathbf{I}_{T}\bigr{)}\widehat{\mathbf{u}}_{i}$
$\displaystyle=\bigl{(}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}-\mathbf{I}_{T}\bigr{)}\bigl{(}\widetilde{\mathbf{u}}_{i}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\bigr{)}.$
(32)

Using Shwarz Inequality and Triangle Inequality for vectors and matrices (B.10
and B.13) in Hansen, (2019, p.795), Equation (32) can be rewritten as

$\displaystyle\|\widehat{\mathbf{v}}_{i}-\widehat{\mathbf{u}}_{i}\|$
$\displaystyle=\big{\|}\bigl{(}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}-\mathbf{I}_{T}\bigr{)}\bigl{(}\widetilde{\mathbf{u}}_{i}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\bigr{)}\big{\|}$
$\displaystyle\leq\big{\|}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}-\mathbf{I}_{T}\big{\|}\big{\|}\bigl{(}\widetilde{\mathbf{u}}_{i}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\bigr{)}\big{\|}$
(33)

Using Woodbury’s formula
$(\mathrm{A}+\mathrm{B}\mathrm{D}\mathrm{C})^{-1}=\mathrm{A}^{-1}-\mathrm{A}^{-1}\mathrm{B}(\mathrm{D}^{-1}+\mathrm{C}\mathrm{A}^{-1}\mathrm{B})^{-1}\mathrm{C}\mathrm{A}^{-1}$
with $\mathrm{A}=\mathbf{I}_{T}$, $\mathrm{B}=\widetilde{\mathbf{X}}_{i}$,
$\mathrm{C}=\widetilde{\mathbf{X}}_{i}^{\prime}$, and
$\mathrm{D}=\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}$,
$(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}$ can be rewritten as follows

$(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}=\big{(}\mathbf{I}_{T}-\widetilde{\mathbf{X}}_{i}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}})^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\big{)}^{-1}=\mathbf{I}_{T}+\widetilde{\mathbf{X}}_{i}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i})^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}$
(34)

and, hence, the first component in (33)

$\begin{split}\big{\|}(\mathbf{I}_{T}-\mathbf{H}_{i})^{-1}-\mathbf{I}_{T}\big{\|}&=\big{\|}\widetilde{\mathbf{X}}_{i}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i})^{-1}\widetilde{\mathbf{X}}_{i}^{\prime}\big{\|}\\\
&\leq\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\bigl{(}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\bigr{)}^{-1}\big{\|}\\\
&=\frac{1}{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\,\Bigg{\|}\biggl{(}\frac{1}{N}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\frac{1}{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\biggr{)}^{-1}\Bigg{\|}.\end{split}$
(35)

Then, using expression (35) in inequality (33)

$\displaystyle\|\widehat{\mathbf{v}}_{i}-\widehat{\mathbf{u}}_{i}\|$
$\displaystyle\leq\frac{1}{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\,\Bigg{\|}\biggl{(}\frac{1}{N}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}-\frac{1}{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\biggr{)}^{-1}\Bigg{\|}\bigl{(}\|\widetilde{\mathbf{u}}_{i}\|+\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}\big{\|}\widehat{\bm{\beta}}-\bm{\beta}\big{\|}\bigr{)}$
(36) $\displaystyle=o_{p}(1)$

where the first component of (36) is $O_{p}(N^{-1})$ under asm.4.asm.4i for
$r\geq 2$; the second component involves $\mathbf{S}_{XX}+o_{p}(1)$ as
$N^{-1}\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}}\overset{p}{\to}\mathbb{E}(\widetilde{\mathbf{X}}^{\prime}\widetilde{\mathbf{X}})=\mathbf{S}_{XX}>0$
by the Central Limit Theorem; the last component is $O_{p}(1)$ because the
random variables in parenthesis are $O_{p}(1)$ under asm.4.asm.4i-asm.4ii, and
$\big{\|}\widehat{\bm{\beta}}-\bm{\beta}\big{\|}\overset{p}{\to}0$ is
$o_{p}(1)$. Therefore, the overall expression is bounded above by an
$o_{p}(1)$ random variable. Note that

$\|\widehat{\mathbf{u}}_{i}-\widetilde{\mathbf{u}}_{i}\|\leq\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}\big{\|}\widehat{\bm{\beta}}-\bm{\beta}\big{\|}$
(37)

because $\widetilde{\mathbf{X}}_{i}$ is $O_{p}(1)$ by asm.4.asm.4i and
$\big{\|}\widehat{\bm{\beta}}-\bm{\beta}\big{\|}\overset{p}{\to}0$,
$\|\widehat{\mathbf{u}}_{i}-\widetilde{\mathbf{u}}_{i}\|$ is $o_{p}(1)$.
Therefore,
$\widehat{\mathbf{u}}_{i}\overset{p}{\to}\widetilde{\mathbf{u}}_{i}$ as
$N\to\infty$ and $T$ fixed. Using result in (36) and (37), we obtain the
desired result

$\widehat{\mathbf{v}}_{i}=\widehat{\mathbf{u}}_{i}+o_{p}(1)\overset{p}{\to}\widetilde{\mathbf{u}}_{i}\,\,\text{
as }N\to\infty\text{ and }T\text{ fixed}.$ (38)

Result (38) shows that the transformed standard errors
$\widehat{\mathbf{v}}_{i}$ are a uniformly consistent estimator for the error
term $\widetilde{\mathbf{u}}_{i}$. This result guarantees the consistency of
any other formula of alternative HC estimators, such as $PHC3$ and $PHC6$. In
fact, using Equation (38) we can show that

$\displaystyle\widehat{\widehat{\mathbf{V}}}_{N}-\widehat{\mathbf{V}}_{N}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}-\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\Big{(}\widehat{\mathbf{v}}_{i}(\widehat{\mathbf{v}}_{i}^{\prime}-\widehat{\mathbf{u}}_{i}^{\prime})(\widehat{\mathbf{v}}_{i}+\widehat{\mathbf{u}}_{i})\widehat{\mathbf{u}}_{i}^{\prime}\Big{)}\widetilde{\mathbf{X}}_{i}$
(39)

Then,

$\displaystyle\Big{\|}\widehat{\widehat{\mathbf{V}}}_{N}-\widehat{\mathbf{V}}_{N}\Big{\|}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}(\widehat{\mathbf{v}}_{i}^{\prime}-\widehat{\mathbf{u}}_{i}^{\prime})\widetilde{\mathbf{X}}_{i}\big{\|}+\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}^{\prime}(\widehat{\mathbf{v}}_{i}+\widehat{\mathbf{u}}_{i})\widehat{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}\big{\|}$
$\displaystyle\leq\max_{1\leq i\leq
N}\|\widehat{\mathbf{v}}_{i}-\widehat{\mathbf{u}}_{i}\|\Biggl{(}\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\big{\|}+\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{u}}_{i}^{\prime}\big{\|}\Biggr{)}$
$\displaystyle\leq\max_{1\leq i\leq
N}\|\widehat{\mathbf{v}}_{i}-\widehat{\mathbf{u}}_{i}\|\Biggl{(}\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\widehat{\mathbf{v}}_{i}\big{\|}+\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\widehat{\mathbf{u}}_{i}\big{\|}\Biggr{)}$
(40) $\displaystyle=o_{p}(1)$

where the first term of (D) is $o_{p}(1)$ from result (36); the two components
in parenthesis are the sums of random variables with finite means both
converging in probability to
$\mathbb{E}\Bigl{(}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}\Bigr{)}$
by asm.4 and results (37)-(38) and, therefore, $O_{p}(1)$. Their product is
$o_{p}(1)$.

The last step left to show is
$\widehat{\mathbf{V}}_{N}\overset{p}{\to}\mathbf{V}$ such that
$\widehat{\widehat{\mathbf{V}}}_{N}=\widehat{\mathbf{V}}_{N}+o_{p}(1)\overset{p}{\to}\mathbf{V}$
as $N\to\infty$ and $T$ fixed. Following Hansen, (2019, pp.230–232), we start
from the definition of conventional robust variance-covariance matrix

$\displaystyle\widehat{\mathbf{V}}_{N}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i}+\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}(\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime})\widetilde{\mathbf{X}}_{i}$
(41)

where the first component of Equation (D) converges in probability to
$\mathbb{E}(\widetilde{\mathbf{X}}_{i}^{\prime}\bm{\Sigma}_{i}\widetilde{\mathbf{X}}_{i})=\mathbf{V}_{i}$
by asm.2.asm.2ii by lie and asm.1 with finite limit $\mathbf{V}$ under thm
6.16 in Hansen, (2019, p.189) for sequences of _inid_ random variables,
provided that asm.5 and asm.4.asm.4i hold. The second component needs to
converge in probability to zero to claim consistency of
$\widehat{\mathbf{V}}_{N}$. Applying matrix norm to (D), we get

$\displaystyle\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}(\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime})\widetilde{\mathbf{X}}_{i}\Bigg{\|}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\Big{\|}\widetilde{\mathbf{X}}_{i}^{\prime}(\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime})\widetilde{\mathbf{X}}_{i}\Big{\|}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}\big{\|}$
(42)

Note that

$\displaystyle\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}$
$\displaystyle=\Bigl{(}\widetilde{\mathbf{u}}_{i}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\Bigr{)}\Bigl{(}\widetilde{\mathbf{u}}_{i}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\Bigr{)}^{\prime}$
$\displaystyle=\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})^{\prime}\widetilde{\mathbf{X}}_{i}^{\prime}-\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\widetilde{\mathbf{u}}_{i}^{\prime}+\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})(\widehat{\bm{\beta}}-\bm{\beta})^{\prime}\widetilde{\mathbf{X}}_{i}^{\prime}$
(43)

Rearranging last line of Equation (D) and using the Triangle Inequality (B.14)
and Schwarz Inequality (B.13), we obtain

$\displaystyle\big{\|}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}\big{\|}$
$\displaystyle\leq
2\Big{\|}\widetilde{\mathbf{X}}_{i}(\widehat{\bm{\beta}}-\bm{\beta})\widetilde{\mathbf{u}}_{i}^{\prime}\Big{\|}+\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}^{2}$
$\displaystyle\leq
2\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}+\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}^{2}$
(44)

Plugging (D) in (D)

$\displaystyle\Bigg{\|}\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}(\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime})\widetilde{\mathbf{X}}_{i}\Bigg{\|}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\big{\|}\widehat{\mathbf{u}}_{i}\widehat{\mathbf{u}}_{i}^{\prime}-\widetilde{\mathbf{u}}_{i}\widetilde{\mathbf{u}}_{i}^{\prime}\big{\|}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\Bigg{\\{}2\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}+\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{2}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}^{2}\Bigg{\\}}$
$\displaystyle\leq
2\Bigg{(}\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{3}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}\Bigg{)}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}+\Bigg{(}\frac{1}{N}\sum_{i=1}^{N}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{4}\Bigg{)}\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}^{2}$
$\displaystyle=o_{p}(1)$ (45)

where the average in the first parenthesis is $O_{p}(1)$ because it is the
mean of random variables bounded above by finite quantities, that is,
$\mathbb{E}\big{(}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{3}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}\big{)}\leq\big{(}\mathbb{E}\big{\|}\widetilde{\mathbf{X}}_{i}\big{\|}^{3}\big{)}^{\frac{3}{4}}\big{(}\mathbb{E}\big{\|}\widetilde{\mathbf{u}}_{i}\big{\|}^{4}\big{)}^{\frac{1}{4}}$
by evoking Hölder’s Inequality (B.28) in Hansen, (2019, p. 796) and under
asm.1 and asm.4.asm.4i-asm.4ii; the average in the second parenthesis is
$O_{p}(1)$ as the average of a random variable with finite mean by
asm.4.asm.4i;
$\Big{\|}\widehat{\bm{\beta}}-\bm{\beta}\Big{\|}\overset{p}{\to}0$ and, thus,
is $o_{p}(1)$. It follows that
$\widehat{\mathbf{V}}_{N}\overset{p}{\to}\mathbf{V}$ and, therefore, the
desired result

$\widehat{\widehat{\mathbf{V}}}_{N}=\widehat{\mathbf{V}}_{N}+o_{p}(1)\to\mathbf{V}.$
(46)

##  Appendix E Consistency of PHC6

The proposed estimator, PHC6, of the asymptotic variance-covariance matrix is

$\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{6}=c_{6}\,\mathbf{S}_{N}^{-1}\widehat{\mathbf{V}}_{N}^{6}\mathbf{S}_{N}^{-1}\,,$
(47)

where the variance-covariance matrix is
$\widehat{\mathbf{V}}_{N}^{6}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{\mathbf{X}}_{i}^{\prime}\widehat{\mathbf{v}}_{i}\widehat{\mathbf{v}}_{i}^{\prime}\widetilde{\mathbf{X}}_{i},$
and the matrix $\mathbf{M}_{i}$ has functional form

$\mathbf{M}_{i}=\begin{cases}\mathbf{I}_{T}&\text{if }h^{*}_{i}<2\\\
\mathbf{I}_{T}-\mathbf{H}_{i}&\text{otherwise}\end{cases}$ (48)

where
$h^{*}_{i}=max\big{\\{}h_{i11}/\overline{h}_{11},\dots,h_{iTT}/\overline{h}_{TT}\big{\\}}$
is the maximal individual leverage of unit $i$; and
$\overline{h}_{tt}=N^{-1}\sum_{i=i}^{N}h_{itt}$ is the average leverage at
time $t$, with $h_{itt}$ being the individual leverage of unit $i$ at time
$t$. The finite sample correction of PHC6 is

$c_{6}=\begin{cases}\frac{(NT-1)N}{(NT-k)(N-1)}&\text{if }h^{*}_{i}<2\\\
\frac{N-1}{N}&\text{otherwise}\end{cases}\\\ $

As $N\to\infty$ and $T$ is fixed, $\mathbf{M}=(\mathbf{I}_{T}-\mathbf{H}_{i})$
because the selection criterion is $2<\infty$. As argued above,
$\mathbf{H}_{i}\overset{p}{\to}\bm{0}$ as $N\to\infty$ and $T$ fixed because
leverage measures are asymptotically negligible (Hansen,, 2019, p.249).
Therefore, PHC6 collapses to PHC3 that converges to PHC0 which is a consistent
estimator of the asymptotic variance (Hansen,, 2019, Theorem 7.7, p.232). From
White,’s (1980) general result and under the above model assumptions and thm
7.7 in Hansen, (2019, p.232), it follows
$\widetilde{\bm{\Sigma}}=\widehat{\bm{\Sigma}}+o_{p}(1)\to\bm{\Sigma}$ as
$N\to\infty$ and $T$ fixed such that
$\widehat{\mathrm{AVar}(\widehat{\bm{\beta}})}_{6}\overset{p}{\to}\mathrm{AVar}(\widehat{\bm{\beta}})$,
making PHC6 consistent estimator of the sampling variance.

##  Appendix F Derivation of the Distribution of $\mathrm{W}$

The error term $u_{it}$ is intrinsically heteroskedastic but not on average
due to the presence of the scaling factor $z(\gamma)$. Let
$\mathrm{W}=\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}$ with
$\\{x_{it,j}\\}_{j=1}^{J}$. When $\gamma=1$, the mean and variance of a random
variable W with an unknown distribution are as follows

$\displaystyle\mathbb{E}(\mathrm{W})$
$\displaystyle=\mathbb{E}\Biggl{[}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\Biggr{]}=\beta_{0}+\sum_{j=1}^{J}\beta_{j}\,\mathbb{E}\bigl{(}x_{it,j}\bigr{)}=\beta_{0}+\sum_{j=1}^{J}\beta_{j}\mu_{x_{j}}$
(49) $\displaystyle\mathrm{Var}\bigl{(}\mathrm{W}\bigl{)}$
$\displaystyle=\mathrm{Var}\Biggl{[}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\Biggr{]}=\sum_{j=1}^{J}\beta^{2}_{j}\,\mathrm{Var}\bigl{(}x_{it,j}\bigr{)}+2\sum_{\begin{subarray}{c}j,k=1\\\
j\neq k\end{subarray}}^{J}\beta_{j}\beta_{k}\mathrm{Cov}(x_{it,j},x_{it,k})$
$\displaystyle=\sum_{j=1}^{J}\beta^{2}_{j}\sigma^{2}_{x_{j}},$ (50)

where $\mathbb{E}\bigl{(}x_{it,j}\bigr{)}=\mu_{x_{j}}$,
$\mathrm{Var}\bigl{(}x_{it,j}\bigr{)}=\sigma^{2}_{x_{j}}$, and
$\mathrm{Cov}(x_{it,j},x_{it,k})=0$ because the independence assumption
guarantees that
$\mathbb{E}\bigl{(}x_{it,j},x_{it,k}\bigr{)}=\mathbb{E}\bigl{(}x_{it,j}\bigr{)}\mathbb{E}\bigl{(}x_{it,k}\bigr{)}$.
The results are valid under independent and identically distributed (_iid_)
random variables. By the assumptions of _iid_ and normality of
$\mathbf{x}_{it}$, the random variable $\mathrm{W}$ is normally distributed
with mean (49) and variance (50). When the regressors are drawn from a
standard normal distribution, (49) and (50) reduce to $\beta_{0}$ and
$\sum_{j=1}^{J}\beta^{2}_{j}$, respectively. Thus,
$\mathrm{W}\mathbin{\overset{}{\kern
0.0pt\sim}}\mathcal{N}\Bigl{(}\beta_{0},\sum_{j=1}^{J}\beta^{2}_{j}\Bigr{)}$.
Standardising $\mathrm{W}$, we get
$\frac{\mathrm{W}-\mu_{w}}{\sigma^{2}_{w}}\mathbin{\overset{}{\kern
0.0pt\sim}}\mathcal{N}\bigl{(}0,1\bigr{)}$.

When $\gamma=2$, the mean and variance of $\mathrm{W}$ are as follows

$\displaystyle\mathbb{E}(\mathrm{W}^{2})$
$\displaystyle=\mathbb{E}\Biggl{[}\Biggl{(}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\Biggr{)}^{2}\Biggr{]}$
$\displaystyle=\beta^{2}_{0}+\sum_{j=1}^{J}\beta^{2}_{j}\,\mathbb{E}\bigl{(}x^{2}_{it,j}\bigr{)}+2\beta_{0}\sum_{j=1}^{J}\beta_{j}\,\mathbb{E}\bigl{(}x_{it,j}\bigr{)}+2\sum_{\begin{subarray}{c}j,k=1\\\
j\neq
k\end{subarray}}^{J}\beta_{j}\beta_{k}\mathbb{E}\bigl{(}x_{it,j}\bigr{)}\mathbb{E}\bigl{(}x_{it,j}\bigr{)}$
(51)
$\displaystyle=\beta^{2}_{0}+\sum_{j=1}^{J}\beta^{2}_{j}(\sigma^{2}_{x_{j}}+\mu^{2}_{x_{j}})+2\beta_{0}\sum_{j=1}^{J}\beta_{j}\mu_{x_{j}}+2\sum_{\begin{subarray}{c}j,k=1\\\
j\neq k\end{subarray}}^{J}\beta_{j}\beta_{k}\mu_{x_{j}}\mu_{x_{j}}$
$\displaystyle\mathrm{Var}\bigl{(}\mathrm{W}^{2}\bigl{)}$
$\displaystyle=\mathrm{Var}\Biggl{[}\biggl{(}\beta_{0}+\sum_{j=1}^{J}\beta_{j}x_{it,j}\biggr{)}^{2}\Biggr{]}$
$\displaystyle=\sum_{j=1}^{J}\beta^{4}_{k}\,\mathrm{Var}\bigl{(}x^{2}_{it,j}\bigr{)}+4\beta^{2}_{0}\sum_{j=1}^{J}\beta^{2}_{j}\,\mathrm{Var}\bigl{(}x_{it,j}\bigr{)}+4\sum_{\begin{subarray}{c}j,k=1\\\
j\neq
k\end{subarray}}^{J}\beta^{2}_{j}\beta^{2}_{k}\mathrm{Var}\bigl{(}x_{it,j},x_{it,j}\bigr{)}$
$\displaystyle=\sum_{j=1}^{J}\beta^{4}_{k}\sigma^{2}_{x_{j}}+4\beta^{2}_{0}\sum_{j=1}^{J}\beta^{2}_{j}\sigma^{2}_{x_{j}}+4\sum_{\begin{subarray}{c}j,k=1\\\
j\neq
k\end{subarray}}^{J}\beta^{2}_{j}\beta^{2}_{k}\bigl{(}\sigma^{2}_{x_{j}}\sigma^{2}_{x_{j}}+\sigma^{2}_{x_{j}}\mu^{2}_{x_{j}}+\sigma^{2}_{x_{j}}\mu^{2}_{x_{j}}\bigr{)}$
(52)

where
$\mathbb{E}\bigl{(}x^{2}_{it,j}\bigr{)}=\sigma^{2}_{x_{j}}+\mu^{2}_{x_{j}}$,
$\mathbb{E}\bigl{(}x_{it,j},x_{it,k}\bigr{)}=\mathbb{E}\bigl{(}x_{it,j}\bigr{)}\mathbb{E}\bigl{(}x_{it,k}\bigr{)}=\mu_{x_{j}}\mu_{x_{k}}$,
and
$\mathrm{Var}\bigl{(}x_{it,j},x_{it,k}\bigr{)}=\bigl{[}\mathbb{E}\bigl{(}x_{it,j},x_{it,ks}\bigr{)}\bigr{]}^{2}-\mathbb{E}\bigl{(}x_{it,j}\bigr{)}^{2}\mathbb{E}\bigl{(}x_{it,ks}\bigr{)}^{2}=\sigma^{2}_{x_{j}}\sigma^{2}_{x_{k}s}+\sigma^{2}_{x_{j}}\mu^{2}_{x_{k}}+\sigma^{2}_{x_{j}}\mu^{2}_{x_{k}}$
because of the _iid_ assumption. Any covariance among variables is null
because of the assumption of independence. With standardised $\mathrm{W}$,
$\Bigr{(}\frac{\mathrm{W}-\mu_{w}}{\sigma^{2}_{w}}\Bigr{)}^{2}\mathbin{\overset{}{\kern
0.0pt\sim}}\chi^{2}_{1}$.