The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
82,056 characters
\vspace{10mm}
\begin{center}
\long\def\symbolfootnote[1]{\begingroup
\footnote[1]{This draft: \today. We would like to thank the Co-editor Xiaohong Chen, an Associate Editor, and two referees for their detailed and constructive comments which have improved the paper considerably. We are also grateful to Antonio Galvao and Matt Harding for helpful comments and suggestions as well as seminar participants at the University of Kentucky, the 2019 CFE/CMStatistics conference, and the 2021 New York Camp Econometrics meeting.}
\endgroup}
{\bf WILD BOOTSTRAP INFERENCE FOR PENALIZED QUANTILE REGRESSION FOR LONGITUDINAL DATA}\symbolfootnote[1]
\long\def\symbolfootnote[2]{\begingroup
\footnote[2]{Carlos Lamarche: Department of Economics, University of Kentucky, 223G Gatton College of Business \& Economics, Lexington, KY 40506. Email: [email removed]. Thomas Parker: Department of Economics, University of Waterloo, 200 University Ave. West, Waterloo, ON, Canada N2L 3G1. Email: [email removed]}
\endgroup}
\vspace{2.5mm}
{\small CARLOS LAMARCHE AND THOMAS PARKER}\symbolfootnote[2]
\end{center}
\vspace{5mm}
{\small
\begin{quote}
{\bf Abstract}: The existing theory of penalized quantile regression for longitudinal data has focused primarily on point estimation. In this work, we investigate statistical inference. We propose a wild residual bootstrap procedure and show that it is asymptotically valid for approximating the distribution of the penalized estimator. The model puts no restrictions on individual effects, and the estimator achieves consistency by letting the shrinkage decay in importance asymptotically. The new method is easy to implement and simulation studies show that it has accurate small sample behavior in comparison with existing procedures. Finally, we illustrate the new approach using U.S. Census data to estimate a model that includes more than eighty thousand parameters.
\bigskip
\noindent \emph{Keywords}: Quantile regression; panel data; penalized estimator; bootstrap inference.
\noindent \emph{JEL classification}: C15; C21; C23.
\end{quote}}
\onehalfspacing
\section{Introduction}
We consider a longitudinal data model of conditional quantiles with individual
intercepts. Variations of this model have been extensively studied in the literature since at least
\citet{NS48}. Recent contributions to the literature using this model for
quantile regression have emphasized the drawbacks of estimating a large number
of individual intercepts ($N$) when the number of time periods ($T$) is small (see Galvao and Kato, 2018,\nocite{Galvao2018} for an excellent survey). Koenker (2004)\nocite{rK04}
proposed an estimator where $N$ individual parameters are regularized by a
Lasso-type penalty, shrinking them towards a common value. As in the case of the Gaussian random effect estimator, shrinkage can reduce the variability of the estimator of the slope parameter in the quantile regression model (Koenker, 2004). In models with short $T$, shrinkage can reduce the bias of the fixed effects estimator of the slope parameter as well (Harding and Lamarche, 2019\nocite{harding2019}).
Although the regularization procedure has advantages, the asymptotic
distribution of the estimator is difficult to approximate. It is known that
Lasso-type estimators have non-standard limiting distributions (Knight and Fu,
2000\nocite{kK00}), but in the case of quantile regression, there are new
challenges. Because individual intercepts are treated as parameters, the
increasing dimension of the parameter vector as the number of units increases
can be an issue. In the case of estimators without regularization, Kato,
Galvao, and Montes-Rojas (2012) and Galvao, Gu, and Volgushev
(2020)\nocite{AGalvao2020} found that $T$ must grow faster than $N$ for
consistency and asymptotic normality at rates that are, at best, similar to
standard non-linear panel data models \citep{Newey2004}. Second, the
covariance matrix of quantile regression estimators typically depends on
conditional densities and the penalized estimator of Koenker (2004) is no
exception. Inference based on the asymptotic distribution requires
non-parametric estimation of nuisance parameters, which can lead to important
size distortions (He, 2018\nocite{He2018}).
Motivated by these limitations, cross-sectional pairs (or block)
bootstrap, which samples sets of covariate and response vectors over
individuals with replacement, appears to be a natural
alternative method for inference. However, we demonstrate that the cross-sectional pairs
bootstrap does not approximate well the limiting distribution of the penalized
estimator. We consider instead a wild residual bootstrap
procedure, which was previously employed by Feng, He, and Hu (2011)\nocite{xumingHe2011}, and Wang, Van Keilegom,
and Maidman (2018)\nocite{lanWang2018} in cross-sectional settings. We investigate the application of the procedure to
longitudinal data and show that the proposed wild bootstrap procedure is a
consistent estimator of the distribution of the penalized estimator.
We begin by deriving consistency and asymptotic normality results for $\ell_1$ penalized
estimators of a longitudinal model in which individual effects can be
correlated with the regressors. Although our model might be considered to
be high-dimensional, the number of parameters is smaller than the number of
observations, as in the pioneering work by Koenker (2004), and thus our results are obtained without
assuming sparsity in terms of the individual intercepts. Consistency and asymptotic normality with $T$
growing faster than $N$ are achieved by letting the penalty parameter that
controls shrinkage diminish in importance asymptotically. Thus, relative to \citet{rK04},
the asymptotic bias of the estimator is zero in our case. The consistency and asymptotic normality results are new — they extend the heuristic
results in Koenker (2004) obtained for a model with individual effects as location shifts and are not
included in \citet{kK12} and \citet{AGalvao2020} because they did not consider
penalized estimation.
The main theoretical contribution is to show that the distribution of the wild
bootstrap estimator consistently estimates the asymptotic distribution and
covariance of the penalized estimator. The results include the special case of no penalization, and thus, these results also show the consistency of the wild bootstrap for the quantile regression estimator with fixed effects. The consistency
of the wild bootstrap is established using developments that are critically
different to those used in Wang, Van Keilegom, and Maidman
(2018)\nocite{lanWang2018}. We also consider bootstrap estimation of the asymptotic
covariance matrix of the slope parameter estimator, which is novel in the panel quantile
literature. As emphasized in \citet{GoncalvesWhite5}, \citet{Andreas2017}, and \citet{HahnLiao21}, the weak convergence of the bootstrap estimator does not necessarily imply convergence of the bootstrap second moment estimator. Therefore, we provide conditions and establish a result that supports using the second moment of the bootstrap distribution to estimate the asymptotic variance of the estimator.
Several penalized estimators for quantile regression models have been proposed
in the literature since \citet{rK04}. \citet{belloni2011} propose quantile regression estimators for
high-dimensional sparse models using cross-sectional data. \citet{WANG2013} considers a penalized least absolute deviation
estimator, and \citet{Wang2019} derives error bounds for the
penalized estimator under weak conditions. \citet{cL06}
investigates the selection of a regularization parameter, and \citet{SLee2018}
study estimation of a high-dimensional quantile regression model with a change
point, or threshold. \citet{harding2016,harding2019} investigate estimation of models with attrition and
correlated random effects. \citet{GU2019} propose a
method for estimation of models with unknown group membership. \citet{CHEN2009} establish the validity of a related weighted bootstrap procedure for the limiting distribution of a penalized sieve estimator and consider applications using quantile regression \citep[see also][]{CHEN2015}. The literature
on penalized estimation methods for linear panel data models has also grown in
the last decade \citep[see, e.g.,][among others.]{kock2013,KOCK2016,aBelloni2016, lSu2016, lSu2018, CANER2018143, kock_tang_2019}
This paper is organized as follows. The next section provides background and discusses the motivation of our study. It also introduces the proposed wild residual bootstrap approach. Section 3 presents theoretical results. Section 4 investigates the small sample performance of the method, showing that the estimator has satisfactory performance under different specifications and it performs better than the cross-sectional pairs bootstrap procedure. Section 5 presents extensions to the basic model. Section 6 illustrates the theory and provides practical guidelines from an application of the method. Considering data from the U.S. Census, we estimate a quantile function with more than eighty thousand parameters to study how wages of U.S. workers have been affected by the North American Free Trade Agreement. Finally, Section 7 concludes. One appendix contains proof of the main results, while a supplementary appendix contains additional technical results and proofs.
\section{Inference for penalized quantile regression}
In this section, we first introduce the model and the estimator, and then we discuss the validity of a cross-sectional pairs bootstrap method. Motivated by the limitations of existing procedures, we propose a new approach to estimate the asymptotic distribution of the estimator.
\subsection{Background and Motivation}
We observe repeated measures $\{ ( y_{it},\bm{x}_{it}') \}_{t=1}^{T}$ for each
subject $1 \leq i \leq N$. The variable $y_{it} \in \mathbb{R}$ denotes the response
for $i$ at time $t$ and $\bm{x}_{it}$ denotes a $p$-dimensional vector of
covariates. Although the number of repeated observations does not vary with $i$,
the analysis can be trivially extended to consider $T_i$ as long as $\max T_i /
\min T_i$ is bounded for $1 \leq i \leq N$ (Gu and Volgushev, 2019).
The model considered in this paper is
\begin{equation}
Q_y(\tau | \bm{x}_{it}) = \bm{x}_{it}' \bm{\beta}_0(\tau) + \alpha_{i0}(\tau),
\end{equation}
where $\tau \in (0,1)$ and $ Q_y(\tau | \bm{x}_{it}) $ is the $\tau$-th quantile of the
conditional distribution of $y_{it}$ given $\bm{x}_{it}$.
It is assumed that the vector $\bm{x}_{it}$ does not contain an intercept.
The parameter of interest is $\bm{\beta}_0(\tau) \in \mathbb{R}^p$ and
$\alpha_{i0}(\tau)$ is treated as a nuisance parameter. Because we consider
just one value of $\tau$, we supress the dependence of the
parameters on $\tau$ in the sequel.
Let $\bm{\theta} = (\bm{\beta}',\bm{\alpha}')' \in \bm{\Theta} \subseteq
\mathbb{R}^{p+N}$, where $\bm{\alpha} = (\alpha_{1},...,\alpha_{N})'$, and let
$\bm{\theta}_0 = (\bm{\beta}_0',\bm{\alpha}_0')'$. To estimate $\bm{\theta}_0$,
we consider the following estimator:
\begin{equation} \label{pqr}
\hat{\bm{\theta}} = (\hat{\bm{\beta}}', \hat{\bm{\alpha}}')'
= \operatorname*{argmin}_{\bm{\theta} \in \bm{\Theta}} \sum_{i=1}^N \sum_{t=1}^T \rho_{\tau}
(y_{it} - \bm{x}_{it}' \bm{\beta} - \alpha_i) + \lambda_T
\sum_{i=1}^{N} | \alpha_i |,
\end{equation}
where $\rho _{\tau }(u)=$ $u(\tau -I(u<0))$ is the quantile regression loss
function. The tuning parameter $\lambda_T \geq 0$ depends on $T$ and it can
also depend on data, as discussed below.
The penalty term in (2.2) helps improve the finite sample performance of the fixed effects estimator, which is defined for $\lambda_T=0$. Shrinkage of the individual effects can lead to reductions of the variance of the estimator. In models with incidental parameters, the penalty term reduces the noise in the estimation of individual intercepts, and consequently, it can also reduce the bias of the fixed effects estimator of $\bm{\beta}_0$. The online appendix presents simulation evidence to illustrate finite sample improvements when the time dimension is short, complementing the evidence presented in Koenker (2004) and Harding and Lamarche (2019). See Bester and Hansen (2009)\nocite{cH2009} for a related penalty approach to bias reduction in nonlinear models with fixed effects.
We establish conditions that result in
a tractable asymptotic distribution for the estimator defined in \eqref{pqr}. However, we expect that
resampling methods offer a more accurate description of the
distribution of the estimator in finite samples. In practice, the cross-sectional pairs bootstrap, which samples over $i$ with replacement keeping the entire
block of time series observations for each $i$, has been used as a
method for inference, primarily in the fixed effects case when $\lambda_T = 0$.
However, the cross-sectional pairs bootstrap does not provide a good
approximation to the sampling distribution of the penalized estimator
\eqref{pqr}, as in the case of the pairs bootstrap procedure for the Lasso
estimator \citep{camponovo15}.
\subsection{A cross-sectional pairs bootstrap procedure}\label{subsection:cs}
We now offer a heuristic illustration of some problems with using a cross-sectional pairs bootstrap and the penalized quantile regression estimator. The cross-sectional pairs bootstrap can be used successfully to estimate the distribution of the quantile regression model with unpenalized fixed effects, but it will be shown below that the penalty causes problems for this approach to resampling. We fix $N$ in this section to avoid the effect of a diverging number of parameters as the sample size increases (later, asymptotic approximations will be found assuming that $T$ grows faster than $N$). This allows us to see problems with the cross-sectional pairs without the additional incidental parameters problem. Define $\bm{\gamma} = (\bm{\delta}', \bm{\eta}')' \in \mathbb{R}^{p + N}$, where $\bm{\delta} = \sqrt{NT} (\bm{\beta} - \bm{\beta}_0)$ and for $i = 1, \ldots N$, $\eta_i = \sqrt{T} (\alpha_i - \alpha_{i0})$. Then let
\begin{equation}
\mathbb{V}_T(\bm{\gamma}) = \sum_{i=1}^N \sum_{t=1}^T \left\{ \rho_{\tau} \left( u_{it} - \frac{\bm{x}_{it}' \bm{\delta}}{\sqrt{NT}} - \frac{\eta_i}{\sqrt{T}} \right) - \rho_{\tau} (u_{it}) \right\} + \lambda_T \sum_{i=1}^{N} \left\{ \left| \alpha_{i0} + \frac{\eta_i}{\sqrt{T}} \right| - | \alpha_{i0} | \right\}, \label{pqr2}
\end{equation}
where $u_{it} = y_{it} - \bm{x}_{it}' \bm{\beta}_0 - \alpha_{i0}$. This objective function is equivalent to~\eqref{pqr}. \citet{kK00} developed a method for dealing with the asymptotic behavior of this objective function, stated here as a lemma.
\begin{lemma}[\citet{kK00}] \label{L1}
Under Assumptions \ref{assume:data}-\ref{assume:Avar} below, if $N$ is fixed, $T
\rightarrow \infty$ and $\lambda_T/\sqrt{T} \to \lambda_0 \geq 0$, the
minimizer of \eqref{pqr2}, $\hat{\bm{\gamma}}$, converges weakly to the
minimizer of $\mathbb{V}: \mathbb{R}^{p + N} \rightarrow \mathbb{R}$ defined by
\begin{equation*}
\mathbb{V}(\bm{\gamma}) = - \bm{\gamma}' \bm{B} + \frac{1}{2} \bm{\gamma}' \bm{D}_1 \bm{\gamma} + \lambda_0 \sum_{i=1}^{N} \left( \eta_i \operatorname*{sgn}(\alpha_{i0}) I(\alpha_{i0} \neq 0) + |\eta_i| I(\alpha_{i0} = 0) \right),
\end{equation*}
where $\bm{D}_1$ is positive definite and $\bm{B} \sim
\mathcal{N}(\bm{\mathbf{0}},\bm{D}_0)$.
\end{lemma}
To examine the validity of the cross-sectional pairs bootstrap, consider an analog loss function for resampled data. Letting $\bm{y}_i$ and $\bm{X}_i$ denote the vector and matrix of response and covariate observations corresponding to unit $i$, a cross-sectional pairs bootstrap procedure resamples $N$ pairs $(\bm{y}_i, \bm{X}_i)$ for $1 \leq i \leq N$ with replacement. Let $n_i^*$ denote the number of times unit $i$ is redrawn from the original sample. Thus, the asymptotic distribution of $\hat{\bm{\gamma}}$ is approximated with $\tilde{\bm{\gamma}} = ( \sqrt{NT} (\tilde{\bm{\beta}} - \hat{\bm{\beta}})', \sqrt{T} ( \tilde{\bm{\alpha}} - \hat{\bm{\alpha}})' )'$ where
\begin{equation} \label{pairblock}
\tilde{\bm{\theta}} = \left( \tilde{\bm{\beta}}', \tilde{\bm{\alpha}}' \right)' = \operatorname*{argmin}_{ \bm{\theta} \in \bm{\Theta} } \sum_{i=1}^N n_i^* \sum_{t=1}^T \rho_{\tau} \left( y_{it} - \bm{x}'_{it} \bm{\beta} - \alpha_i \right) + \lambda_T \sum_{i=1}^N n_i^* | \alpha_i |.
\end{equation}
Since $n_i^*$ is a multinomial weight with probability $1/N$, it is straightforward to calculate that the expected value of the objective function with respect to the bootstrap weights (i.e., conditional on the observations) is minimized at $\hat{\bm{\theta}} = (\hat{\bm{\beta}}, \hat{\bm{\alpha}})$. However, a finite sample problem is associated with the presence of the penalty in the objective function. To see this, let $\alpha_i^\ast = n_i^* |\alpha_i|$ and $\mathcal{A}=\{i : \alpha_i^* \neq 0\}$ denote the ``active" set corresponding to the penalty term in \eqref{pairblock}. In each bootstrap repetition, the cardinality of $\mathcal{A} < N$, leading to solutions $\tilde{\bm{\theta}}$ that can be potentially very different than the minimizer $\hat{\bm{\theta}}$. This may be especially so when $\alpha_i$ is correlated with $\bm{x}_{it}$.
To see other problems with the cross-sectional bootstrap, we can find the weak limit of the bootstrap objective function \eqref{pairblock} similarly to Lemma~\ref{L1}. When we recenter~\eqref{pairblock} employing $\hat{\bm{\theta}}$, using the $i$ chosen by resampling, we find a naive bootstrap analog of the original objective function~\eqref{pqr2}, denoting $\hat{u}_{it} = y_{it} - \hat{\bm{\beta}}' \bm{x}_{it} - \hat{\alpha}_i$:
\begin{equation}
\tilde{\mathbb{V}}_{T}(\bm{\gamma}) = \sum_{i=1}^N n_i^* \sum_{t=1}^T \left\{ \rho_{\tau} \left( \hat{u}_{it} - \frac{\bm{\delta}'\bm{x}_{it}}{\sqrt{NT}} - \frac{\eta_i}{\sqrt{T}} \right) - \rho_{\tau}( \hat{u}_{it} ) \right\} + \lambda_T \sum_{i=1}^N n_i^* \left\{ \left| \hat{\alpha}_i + \frac{\eta_i}{\sqrt{T}} \right| - | \hat{\alpha}_i | \right\}.
\end{equation}
As $T \rightarrow \infty$, assuming $\hat{\eta}_i = \sqrt{T}(\hat{\alpha}_i - \alpha_{i0}) \overset{d}{\longrightarrow} A_i$ for $i = 1, \ldots N$ as $T \rightarrow \infty$, $\tilde{\mathbb{V}}_T$ converges weakly to
\begin{equation*}
\tilde{\mathbb{V}}(\bm{\gamma}) = - \bm{\gamma}' \tilde{\bm{B}} + \frac{1}{2} \bm{\gamma}' \tilde{\bm{D}}_1 \bm{\gamma} + \lambda_0 \sum_{i=1}^N n_i^* \big( \eta_i \operatorname*{sgn}(\alpha_{i0}) I(\alpha_{i0} \neq 0) + \left( |\eta_i + A_i| - |A_i| \right) I(\alpha_{i0} = 0) \big).
\end{equation*}
However, there are two key differences with the resulting expression. The first problem with this limiting objective function is that $\tilde{\bm{B}} \neq \bm{B}$ and $\tilde{\bm{D}_1} \neq \bm{D}_1$ from Lemma~\ref{L1}, due to the fact that recentering uses $\hat{\bm{\theta}}$, which is asymptotically biased if $\lambda_0>0$. Second, there is additional randomness arising from variable selection and resampling. (In the online appendix, we illustrate these issues with fixed $N$ and $T$). In the next section, we propose a wild residual bootstrap that does not suffer from these shortcomings. Then we expect that the distribution of the wild bootstrap estimator $\bm{\gamma}^\ast$ provides a better approximation to the distribution of $\hat{\bm{\gamma}}$ in Lemma \ref{L1}.
\subsection{Wild bootstrap procedures} \label{subsec:wild}
Let $\hat{u}_{it} = y_{it} - \bm{x}_{it}' \hat{\bm{\beta}} - \hat{\alpha}_i$ be the $\tau$-th quantile residual. Let $u_{it}^\ast = w_{it} | \hat{u}_{it} |$ denote bootstrap residuals, where $w_{it}$ is drawn randomly from a pre-determined distribution $G_W$ that satisfies the following conditions:
\begin{assumptionA} \label{C1}
The $\tau$-th quantile of $G_W$ is equal to zero, i.e. $G_W(0) = \tau$.
\end{assumptionA}
\begin{assumptionA} \label{C2}
The support of $G_W$ is bounded and contained in the interval $(-\infty,-c_1] \cup [c_2, \infty)$, where $c_1 > 0$ and $c_2 > 0$.
\end{assumptionA}
\begin{assumptionA} \label{C3}
The weight distribution $G_W$ satisfies $- \int_{-\infty}^0 w^{-1} dG_W(w) = \int_0^{+\infty} w^{-1} dG_W(w) = \frac{1}{2}$.
\end{assumptionA}
Several weight distributions have been proposed in the
quantile regression literature that satisfy these conditions. \citet{xumingHe2011} propose, for $1/8 \leq \tau \leq 7/8$, the continuous weight density
$g_W(w) = - w I(-2 \tau - 1/4 \leq w \leq -2 \tau + 1/4) + w I(2 (1-\tau) - 1/4
\leq w \leq 2 (1-\tau) + 1/4)$. Another distribution that satisfies
\ref{C1}-\ref{C3} is the two-point distribution at $w = 2 (1-\tau)$ with
probability $\tau$ and at $w = - 2 \tau$ with probability $(1-\tau)$. We adopt
this distribution in the numerical examples. See Appendix 3 in
\citet{lanWang2018} for additional examples of the weight
distribution.
Using the bootstrap sample of residuals and the penalized quantile estimator as defined in equation \eqref{pqr}, we can form $y_{it}^\ast = \bm{x}_{it}' \hat{\bm{\beta}} + \hat{\alpha}_i + u_{it}^\ast$ to obtain the bootstrap estimator:
\begin{equation} \label{penboot}
\bm{\theta}^\ast = (\bm{\beta}^{\ast'}, \bm{\alpha}^{\ast'})' = \operatorname*{argmin}_{\bm{\theta} \in \bm{\Theta}} \sum_{i=1}^N \sum_{t=1}^T \rho_{\tau} (y_{it}^\ast - \bm{x}_{it}' \bm{\beta} - \alpha_i) + \lambda_T \sum_{i=1}^{N} | \alpha_{i} | .
\end{equation}
Given a bootstrap sample $\{\bm{\beta}_b^\ast\}_{b=1}^B$, we can obtain confidence intervals that are asymptotically valid, as demonstrated in Theorem \ref{thm:boot} below. Let $G_{j}^*(\alpha/2)$ and $G_{j}^*(1-\alpha/2)$ be the $(\alpha/2)$-th quantile and $(1-\alpha/2)$-th quantile of the bootstrap distribution of $\sqrt{NT} ( \beta_{j}^\ast - \hat{\beta}_{j})$ for $j = 1,2,\hdots,p$. We obtain asymptotically valid $100 (1-\alpha)\%$ confidence intervals for $\beta_{j}$ by $[ \hat{\beta}_{j} - (NT)^{-1/2} G_{j}^*(1 - \alpha/2), \hat{\beta}_{j} - (NT)^{-1/2} G_{j}^*(\alpha/2)]$. Alternatively, Theorem~\ref{thm:var} shows that we may also estimate the covariance matrix of $\sqrt{NT}(\hat{\bm{\beta}} - \bm{\beta}_0)$ using the estimated covariance matrix from the bootstrap sample, which can be used to estimate the variance without requiring density estimation and to construct bootstrap-$t$ statistics for inference.
We may also consider a threshold estimator for $1 \leq i \leq N$, $\alpha_i^{**} = \hat{\alpha}_i I( | \hat{\alpha}_i | \geq a_{T})$, where $a_T$ is a constant that satisfies $a_T \to 0$ as $T \to \infty$. Define $v_{it}^\ast = w_{it} | \hat{v}_{it} |$, where $\hat{v}_{it} = y_{it} - \bm{x}_{it}' \hat{\bm{\beta}} - \alpha_i^{**}$. The response variable is generated as $y_{it}^{**} = \bm{x}_{it}' \hat{\bm{\beta}} + \alpha_i^{**} + v_{it}^\ast$, and the threshold estimator is defined as
\begin{equation} \label{penbootthr}
\bm{\theta}^{**} = \operatorname*{argmin}_{\bm{\theta} \in \bm{\Theta}} \sum_{i=1}^N \sum_{t=1}^T \rho_{\tau} (y_{it}^{**} - \bm{x}_{it}' \bm{\beta} - \alpha_i) + \lambda_T \sum_{i=1}^{N} | \alpha_{i} | .
\end{equation}
As in the case of the estimator defined in \eqref{penboot}, we estimate the
distribution of $\hat{\bm{\theta}}$ based on the estimator $\bm{\theta}^{**}$. Given the similarities between estimators \eqref{penboot} and \eqref{penbootthr}, we derive below consistency and asymptotic normality results for \eqref{penboot} only. The performance of the bootstrap with this estimator is examined in the online appendix.
\subsection{Tuning parameter selection}
The tuning parameter $\lambda_T$ controls the degree of shrinkage of the
individual effect $\alpha_i$ towards zero and the penalty helps to control the bias
and variance of $\hat{\bm{\beta}}$. We restrict the tuning parameter to
$\lambda_T \in \mathcal{L} \subset [0, \lambda_U]$, where $\lambda_U$ is an
upper bound. As shown in Lemma~\ref{lem:upperbnd} in the supplementary
appendix, $\lambda_U = \max\{\tau, 1 - \tau\} T$ is a natural choice because if
$\lambda_T$ is set larger than this value, all the individual effects will be
set equal to zero. If the number of observed time periods $T_i$ vary over $i$,
then one would need to replace the $T$ in these bounds with $\max_i T_i$. This
estimator accommodates the choice of $\lambda_T = 0$, which means that the
results below continue to hold for the corresponding unpenalized estimator.
The selection $\lambda_T$ in related settings has been investigated in several papers \citep[see, e.g.,][]{cL06,ERLee2014,lanWang2018}. We follow \citet{lanWang2018} and employ cross-validation for tuning parameter selection. To the best of our knowledge, theory has not yet been developed for the stochastic order of $\lambda_T$ when chosen using cross-validation, but in extensive simulations we have found that it tends to grow much more slowly than $T$, as required in Theorems~\ref{thm:consistent} and~\ref{thm:AN} below.
\section{Asymptotic theory}
This section investigates the large sample properties of the proposed estimator. We consider the following assumptions:
\begin{assumptionB} \label{assume:data}
Suppose that $\{ (y_{it}, \bm{x}_{it}): t \geq 1 \}$ are independent across $i$ and independent and identically distributed (i.i.d.) within each unit $i$.
\end{assumptionB}
\begin{assumptionB}\label{assume:ID}
For each $\phi > 0$,
\begin{equation*}
\inf_{i\geq1} \inf_{\| \bm{\theta}_i \|_1 = \phi} \textnormal{E} \left[ \int_0^{ (\alpha_i - \alpha_{i0}) + \bm{x}_{it}' (\bm{\beta} - \bm{\beta}_0)} \left( F_i(s | \bm{x}_{it}) - \tau \right) \textnormal{d} s \right] = \epsilon_\phi > 0,
\end{equation*}
where $F_i := F_{u_{it}|\bm{x}_{it}}$ is the distribution function of $u_{it} = y_{it} - \alpha_{i0} - \bm{x}_{it}' \bm{\beta}_0$ conditional on $\bm{x}_{it}$.
\end{assumptionB}
\begin{assumptionB} \label{assume:xsupport}
The covariate vector $\bm{x}_{it}$ satisfies $\sup_{i,t} \| \bm{x}_{it} \| < M < \infty$ a.s.
\end{assumptionB}
These conditions are standard in the literature on quantile regression with individual effects. Conditions \ref{assume:data} and \ref{assume:ID} are the same as Assumptions (A1) and (A3) in Kato, Galvao, and Montes-Rojas (2012). Condition \ref{assume:data} is relaxed in Kato et al. (2012) and in Section 5 below to allow for time dependence. Condition \ref{assume:ID} is an identification condition and it is sufficient for consistency. Slightly weaker than the assumption that $F_i$ has a continuous density given $\bm{x}_{it}$, it allows an expansion that guarantees the convexity of the limiting objective function, and therefore, the uniqueness of $(\bm{\beta}_0',\alpha_{i0})$ for all $1 \leq i \leq N$. Assumption \ref{assume:xsupport} is a simple way to assume appropriate moment conditions on the covariates and it is similar to (B1) in Kato, Galvao, and Montes-Rojas (2012) and (A1) in \citet{GU2019}. The condition can be relaxed as in Kato, Galvao and Montes-Rojas (2012). Condition B3 can be replaced with the moment condition $\sup_{i \geq 1} \textnormal{E} \left[ \| \bm{x}_{i1} \|^{2s} \right] < \infty$ for some $s \geq 1$. The implication of this weaker condition is that $N / T^s \to 0$ instead of $\log(N)/T \to \infty$ to achieve consistency, as demonstrated in Theorem \ref{thm:consistent}.
The consistency of the estimator $\hat{\bm{\theta}}$ is needed to establish the main result stated in Theorem \ref{thm:boot}.
\begin{theorem} \label{thm:consistent}
Under Assumptions \ref{assume:data}-\ref{assume:xsupport}, if $\log(N)/T \to 0$ and $\lambda_T = o_p(T)$ as $N, T \to \infty$, then the estimator $\hat{\bm{\theta}}$ defined in equation \eqref{pqr} is a consistent estimator of $\bm{\theta}_0$.
\end{theorem}
\begin{remark}
Theorem \ref{thm:consistent} is of independent interest as it has not been established the consistency of the penalized estimator under arbitrary dependence between regressors and individual effects. The result depends on the condition that $\lambda_T$, the parameter governing penalization of the individual effects, grows slowly as $T$ increases.
\end{remark}
We now focus our attention on weak convergence and we present a series of results to facilitate the estimation of standard errors and confidence intervals. To show asymptotic normality of the estimator, it is necessary to strengthen the conditions required for consistency slightly with the following conditions routinely adopted in the panel quantile regression literature (see, e.g., assumptions (B2) and (B3) in Kato, Galvao, and Montes-Rojas, 2012, and assumption (A2) in Gu and Volgushev, 2019).
\begin{assumptionB}\label{assume:ID_weakconv}
The conditional density function $f_i := f_{u_{it}|\bm{x}_{it}}$ corresponding to $F_i$ is uniformly bounded and has a bounded first derivative:
\begin{equation*}
\overline{f} := \sup_i \sup_{u \in \mathbb{R}, \bm{x} \in \mathbb{R}^p} | f_i(u | \bm{x}) | < \infty
\end{equation*}
and
\begin{equation*}
\overline{f'} := \sup_i \sup_{u \in \mathbb{R}, \bm{x} \in \mathbb{R}^p} | f'_i(u | \bm{x}) | < \infty.
\end{equation*}
Assume that in an open neighborhood $\mathcal{U}$ of $0$, $f_i$ is bounded away from zero for all realizations of $\bm{x}_{it}$:
\begin{equation*}
\underline{f} := \inf_i \inf_{u \in \mathcal{U}, \bm{x} \in \mathbb{R}^p} | f_i(u | \bm{x}) | < \infty.
\end{equation*}
\end{assumptionB}
\begin{assumptionB} \label{assume:Avar}
Let $\varphi_i := \textnormal{E} \left[ f_i(0 | \bm{x}_{i1}) \right]$, $\bm{E}_i := \textnormal{E} \left[ f_i(0 | \bm{x}_{i1}) \bm{x}_{i1} \right]$ and $\bm{J}_i := \textnormal{E} \left[ f_i(0 | \bm{x}_{i1}) \bm{x}_{i1} \bm{x}_{i1}' \right]$. Let
\begin{equation*}
\bm{D}_N = \frac{1}{N} \sum_{i=1}^N \left( \bm{J}_i - \varphi_i^{-1} \bm{E}_i \bm{E}_i' \right).
\end{equation*}
Suppose that $\bm{D}_N$ is positive definite for all $N$ and there is a positive definite matrix $\bm{D}$ such that $\bm{D} = \lim_{N \rightarrow \infty} \bm{D}_N$. Also assume that
\begin{equation*}
\bm{V} = \tau(1-\tau) \times \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N \textnormal{E} \left[ \left( \bm{x}_{i1} - \varphi_i^{-1} \bm{E}_i \right) \left( \bm{x}_{i1} - \varphi_i^{-1} \bm{E}_i \right)' \right]
\end{equation*}
is positive definite.
\end{assumptionB}
Then we have the following result:
\begin{theorem} \label{thm:AN}
Under Assumptions \ref{assume:data}-\ref{assume:Avar}, if $N^2 (\log N)^3 / T \to 0$ and $\lambda_T = o_p(T^{1/2} (\log N)^{1/2})$ as $N, T \rightarrow \infty$, then
\begin{equation*}
\sqrt{NT} (\hat{\bm{\beta}} - \bm{\beta}_0) \overset{d}{\longrightarrow} \mathcal{N}(\bm{0}, \bm{\Omega}),
\end{equation*}
where $\bm{\Omega} = \bm{D}^{-1} \bm{V} \bm{D}^{-1}$.
\end{theorem}
\begin{remark}
As with Condition G in Theorem 3.2 in \citet{GU2019}, Theorem \ref{thm:AN} provides a selection rule for candidate values of the tuning parameters that are justified by theory. The limiting distribution for this estimator matches that of the conventional fixed effects estimator derived in \citet{kK12} because the tuning parameter $\lambda_T$ diverges at a slow rate.
\end{remark}
\begin{remark}
Because the goal of the shrinkage estimator here is not variable selection but regularization of the estimated $\hat{\alpha}_i$, the rate of growth of $\lambda_T$ is different than what would usually be used in high-dimensional models (\citet[p. 86]{belloni2011}, \citet[eq. 2.3]{SLee2018}, and \citet[Theorem 3.2]{Wang2019}). This difference in stochastic order is because the individual effects $\{\alpha_{i0}\}_i$ are not assumed sparse and this condition on $\lambda_T$ is needed for consistency in models with regressors correlated with individual latent effects.
Moreover, perhaps not surprisingly, the rates derived for linear models \citep[see, e.g.,][]{kock2013,KOCK2016} are also different to the rate required for establishing the asymptotic normality of the quantile estimator.
\end{remark}
The wild residual bootstrap procedure is consistent as an estimator of the asymptotic distribution of $\hat{\bm{\beta}}$, as the next theorem shows.
\begin{theorem} \label{thm:boot}
Under Assumptions \ref{C1}-\ref{C3} and the conditions of Theorem~\ref{thm:AN},
\begin{equation*}
\sup_{b \in \mathbb{R}^p} \left| \textnormal{P} \left\{ \sqrt{NT}( \bm{\beta}^* - \hat{\bm{\beta}}) \leq b | \bm{S} \right\} - \textnormal{P} \left\{ \sqrt{NT}( \hat{\bm{\beta}} - \bm{\beta}_0) \leq b \right\} \right| \overset{p}{\longrightarrow} 0
\end{equation*}
where $\bm{S}$ denotes the observed sample and $\bm{\beta}^\ast$ denotes the slope estimator defined by~\eqref{penboot}.
\end{theorem}
\begin{remark}
By setting $\lambda_T = 0$, Theorem~\ref{thm:boot} also implies consistency of the wild residual bootstrap for the unpenalized estimator with individual effects and i.i.d. errors.
\end{remark}
\begin{remark} \label{rk:bootlmda}
The results allow for a data-dependent $\lambda_T$ but they do not allow selecting the tuning parameter at each bootstrap repetition. While theoretical developments are out of the scope of this paper, we investigated if this idea leads to improvements in the finite sample performance of the estimator. We did not find significant changes relative to the results presented in Section 4, although the computational cost of the procedure is higher.
\end{remark}
Theorem~\ref{thm:boot} only shows consistency of the bootstrap distribution estimator. Theorem~\ref{thm:var} ahead shows that the bootstrap covariance matrix, defined as
\begin{equation*}
\bm{\Omega}^* = \textnormal{E}^* \left[ NT \left( \bm{\beta}^* - \hat{\bm{\beta}} \right) \left( \bm{\beta}^* - \hat{\bm{\beta}} \right)' \right],
\end{equation*}
may be used to estimate the covariance of $\sqrt{NT}(\hat{\bm{\beta}} - \bm{\beta}_0)$. In practice, one simply uses the sample covariance of all the bootstrap repetitions, increasing the number of repetitions to bring the sample average as close as desired to the bootstrap expectation. Variance estimation using the bootstrap was formally investigated for quantile regression with clustered data in \citet{Andreas2017}, but the model in this paper is complicated by the diverging number of individual effects as $N \rightarrow \infty$ and the penalty term in~\eqref{penboot}.
\begin{theorem} \label{thm:var}
Under Assumptions \ref{C1}-\ref{C3} and the conditions of Theorem~\ref{thm:AN}, if $\bm{\theta}_i$ for $1 \leq i \leq N$ lie in a compact set and
\begin{equation*}
\sup_{N,T} \textnormal{E} \left[ | \sqrt{N} \lambda_T / \sqrt{T} |^q \right] < \infty
\end{equation*}
for $q > 2$, then
\begin{equation*}
\| \bm{\Omega}^* - \bm{\Omega} \| \stackrel{p^\ast}{\longrightarrow} 0.
\end{equation*}
\end{theorem}
The assumptions that are required for Theorem~\ref{thm:var} are slightly stronger than those used in Theorem~\ref{thm:boot}. The requirement on $\lambda_T$ is due to its presence in asymptotic expansions leading to the Bahadur representation of $\bm{\beta}^\ast$ and is similar to the moment requirement made on the covariates in \citet{Andreas2017}. The compactness assumption must be made to ensure that expansions used in the asymptotic approximation are uniformly bounded.
\section{Simulation Study}
In this section, we report the results of several simulation experiments designed to evaluate the performance of the method in finite samples. We consider a data generating process similar to the ones considered in Koenker (2004) and
Kato, Galvao and Montes-Rojas (2012). The dependent variable is $y_{it} = \alpha_i + x_{it} + (1 + \zeta x_{it}) u_{it}$, where $x_{it} = 0.5 \alpha_i + z_i + \epsilon_{it}$, and $z_i$
and $\epsilon_{it}$ are i.i.d. random variables distributed as $\chi^2$ with 3 degrees of freedom ($\chi_3^2$). The corresponding quantile regression function is
$Q_y (\tau | x_{it}) = \alpha_{0i} + \beta_0 x_{it}$, where $\alpha_{0i} = \alpha_i + F_u(\tau)^{-1}$, $\beta_0 = 1 + \zeta F_u(\tau)^{-1}$, and $F_u(\cdot)$ denotes the distribution of the error term, $u_{it}$.
We generate data from several variations of the basic model. In one variant of the model, $\alpha_i$ is an i.i.d. Gaussian random variable.
In another, we generate $\alpha_i = i/N$ for $1 \leq i \leq N$ as in Galvao, Gu, and Volgushev (2020). We use $\zeta \in \{0,0.5\}$, and thus, $\beta_0 = 1$ in the location shift version of the model and $\beta_0 = 1 + 0.5 F_u(\tau)^{-1}$ in the location-scale shift case. Lastly, we consider three
different distributions for the error term. We assume that $u_{it}$ is
distributed as $\mathcal{N}(0,1)$, a $t$ distribution with 3 degrees of
freedom ($t_3$), or $\chi_3^2$.
\begin{singlespace}
\begin{table}
\begin{center}\footnotesize
\begin{tabular}{c c c c c c c c c c c c c c} \hline
& & \multicolumn{6}{c}{Quantile 0.5} & \multicolumn{6}{c}{Quantile 0.75} \\
& & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} \\ \cmidrule(lr){3-5} \cmidrule(lr){6-8} \cmidrule(lr){9-11} \cmidrule(lr){12-14}
$N$ & $T$ & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} \\ \cmidrule(lr){4-5} \cmidrule(lr){7-8} \cmidrule(lr){10-11} \cmidrule(lr){13-14}
& & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim \mathcal{N}(0,1)$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.697 & 0.902 & 0.905 & 0.675 & 0.909 & 0.913 & 0.702 & 0.854 & 0.854 & 0.640 & 0.848 & 0.859 \\
100 & 10 & 0.720 & 0.908 & 0.868 & 0.683 & 0.915 & 0.876 & 0.737 & 0.887 & 0.900 & 0.683 & 0.885 & 0.903 \\
200 & 5 & 0.677 & 0.923 & 0.925 & 0.670 & 0.916 & 0.920 & 0.668 & 0.862 & 0.861 & 0.641 & 0.873 & 0.877 \\
200 & 10 & 0.650 & 0.925 & 0.857 & 0.700 & 0.927 & 0.870 & 0.662 & 0.898 & 0.909 & 0.691 & 0.897 & 0.915 \\ \hline
25 & 50 & 0.886 & 0.908 & 0.904 & 0.779 & 0.887 & 0.882 & 0.857 & 0.881 & 0.880 & 0.758 & 0.885 & 0.888 \\
25 & 100 & 0.903 & 0.910 & 0.911 & 0.811 & 0.902 & 0.905 & 0.908 & 0.898 & 0.901 & 0.816 & 0.892 & 0.898 \\
50 & 50 & 0.833 & 0.905 & 0.903 & 0.769 & 0.884 & 0.880 & 0.831 & 0.893 & 0.888 & 0.768 & 0.889 & 0.892 \\
50 & 100 & 0.847 & 0.900 & 0.895 & 0.831 & 0.903 & 0.900 & 0.854 & 0.902 & 0.897 & 0.827 & 0.898 & 0.902 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim t_3$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.710 & 0.906 & 0.911 & 0.674 & 0.902 & 0.912 & 0.701 & 0.828 & 0.833 & 0.623 & 0.819 & 0.840 \\
100 & 10 & 0.713 & 0.923 & 0.880 & 0.660 & 0.919 & 0.881 & 0.734 & 0.881 & 0.893 & 0.680 & 0.864 & 0.880 \\
200 & 5 & 0.681 & 0.932 & 0.936 & 0.645 & 0.907 & 0.927 & 0.669 & 0.841 & 0.852 & 0.573 & 0.816 & 0.852 \\
200 & 10 & 0.650 & 0.922 & 0.834 & 0.661 & 0.931 & 0.845 & 0.641 & 0.881 & 0.892 & 0.678 & 0.859 & 0.889 \\ \hline
25 & 50 & 0.887 & 0.921 & 0.916 & 0.784 & 0.906 & 0.906 & 0.852 & 0.887 & 0.886 & 0.738 & 0.881 & 0.881 \\
25 & 100 & 0.905 & 0.901 & 0.901 & 0.819 & 0.892 & 0.891 & 0.870 & 0.883 & 0.891 & 0.801 & 0.891 & 0.900 \\
50 & 50 & 0.850 & 0.898 & 0.895 & 0.759 & 0.884 & 0.884 & 0.839 & 0.886 & 0.889 & 0.761 & 0.900 & 0.895 \\
50 & 100 & 0.848 & 0.886 & 0.887 & 0.816 & 0.899 & 0.892 & 0.837 & 0.885 & 0.890 & 0.770 & 0.863 & 0.875 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim \chi_3^2$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.730 & 0.906 & 0.912 & 0.633 & 0.897 & 0.907 & 0.690 & 0.728 & 0.745 & 0.589 & 0.716 & 0.716 \\
100 & 10 & 0.673 & 0.894 & 0.871 & 0.647 & 0.870 & 0.855 & 0.718 & 0.828 & 0.844 & 0.707 & 0.827 & 0.838 \\
200 & 5 & 0.708 & 0.920 & 0.927 & 0.657 & 0.915 & 0.910 & 0.652 & 0.745 & 0.753 & 0.602 & 0.763 & 0.763 \\
200 & 10 & 0.686 & 0.901 & 0.860 & 0.642 & 0.880 & 0.838 & 0.703 & 0.833 & 0.853 & 0.703 & 0.818 & 0.831 \\ \hline
25 & 50 & 0.791 & 0.882 & 0.883 & 0.716 & 0.891 & 0.896 & 0.749 & 0.839 & 0.839 & 0.707 & 0.861 & 0.862 \\
25 & 100 & 0.845 & 0.887 & 0.890 & 0.749 & 0.873 & 0.874 & 0.781 & 0.871 & 0.879 & 0.683 & 0.856 & 0.860 \\
50 & 50 & 0.768 & 0.869 & 0.872 & 0.726 & 0.898 & 0.897 & 0.758 & 0.862 & 0.862 & 0.727 & 0.876 & 0.877 \\
50 & 100 & 0.798 & 0.879 & 0.880 & 0.735 & 0.891 & 0.892 & 0.770 & 0.861 & 0.867 & 0.710 & 0.864 & 0.877 \\ \hline
\end{tabular}
\vspace{3mm}
\end{center}
\caption{\emph{Empirical coverage probabilities of the bootstrap confidence interval for a nominal 90\% level. A location shift model is considered. CS denotes cross-sectional pairs bootstrap, WB denotes wild bootstrap, PQR denotes the penalized estimator, and FE is the unpenalized fixed effects estimator.}}
\label{mc.table1}
\end{table}
\end{singlespace}
\begin{singlespace}
\begin{table}
\begin{center}\footnotesize
\begin{tabular}{c c c c c c c c c c c c c c} \hline
& & \multicolumn{6}{c}{Quantile 0.5} & \multicolumn{6}{c}{Quantile 0.75} \\
& & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} \\ \cmidrule(lr){3-5} \cmidrule(lr){6-8} \cmidrule(lr){9-11} \cmidrule(lr){12-14}
$N$ & $T$ & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} \\ \cmidrule(lr){4-5} \cmidrule(lr){7-8} \cmidrule(lr){10-11} \cmidrule(lr){13-14}
& & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale shift model ($\zeta=0.5$) and $u \sim \mathcal{N}(0,1)$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.688 & 0.861 & 0.881 & 0.699 & 0.884 & 0.893 & 0.592 & 0.798 & 0.824 & 0.690 & 0.819 & 0.827 \\
100 & 10 & 0.619 & 0.899 & 0.865 & 0.657 & 0.896 & 0.868 & 0.611 & 0.864 & 0.867 & 0.652 & 0.877 & 0.887 \\
200 & 5 & 0.660 & 0.871 & 0.892 & 0.674 & 0.868 & 0.892 & 0.540 & 0.825 & 0.823 & 0.605 & 0.850 & 0.848 \\
200 & 10 & 0.650 & 0.913 & 0.852 & 0.651 & 0.904 & 0.860 & 0.567 & 0.876 & 0.880 & 0.627 & 0.868 & 0.884 \\ \hline
25 & 50 & 0.698 & 0.898 & 0.893 & 0.676 & 0.900 & 0.899 & 0.696 & 0.881 & 0.882 & 0.672 & 0.866 & 0.865 \\
25 & 100 & 0.762 & 0.906 & 0.906 & 0.678 & 0.894 & 0.896 & 0.749 & 0.894 & 0.897 & 0.678 & 0.892 & 0.893 \\
50 & 50 & 0.685 & 0.896 & 0.892 & 0.678 & 0.888 & 0.883 & 0.679 & 0.888 & 0.889 & 0.669 & 0.886 & 0.888 \\
50 & 100 & 0.720 & 0.900 & 0.902 & 0.674 & 0.905 & 0.903 & 0.697 & 0.879 & 0.881 & 0.677 & 0.901 & 0.902 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale shift model ($\zeta=0.5$) and $u \sim t_3$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.723 & 0.857 & 0.886 & 0.758 & 0.866 & 0.890 & 0.620 & 0.751 & 0.769 & 0.721 & 0.769 & 0.781 \\
100 & 10 & 0.664 & 0.902 & 0.882 & 0.676 & 0.895 & 0.873 & 0.638 & 0.862 & 0.874 & 0.672 & 0.855 & 0.857 \\
200 & 5 & 0.754 & 0.875 & 0.910 & 0.722 & 0.873 & 0.909 & 0.513 & 0.767 & 0.775 & 0.707 & 0.780 & 0.795 \\
200 & 10 & 0.650 & 0.899 & 0.857 & 0.649 & 0.910 & 0.843 & 0.554 & 0.816 & 0.835 & 0.634 & 0.822 & 0.824 \\ \hline
25 & 50 & 0.742 & 0.910 & 0.910 & 0.682 & 0.914 & 0.916 & 0.694 & 0.874 & 0.875 & 0.672 & 0.869 & 0.873 \\
25 & 100 & 0.750 & 0.885 & 0.884 & 0.678 & 0.896 & 0.897 & 0.711 & 0.874 & 0.876 & 0.688 & 0.885 & 0.888 \\
50 & 50 & 0.683 & 0.891 & 0.888 & 0.675 & 0.880 & 0.877 & 0.690 & 0.881 & 0.879 & 0.678 & 0.885 & 0.891 \\
50 & 100 & 0.714 & 0.886 & 0.885 & 0.661 & 0.883 & 0.883 & 0.691 & 0.879 & 0.889 & 0.668 & 0.864 & 0.872 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale model ($\zeta=0.5$) and $u \sim \chi_3^2$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.735 & 0.840 & 0.862 & 0.693 & 0.818 & 0.842 & 0.746 & 0.719 & 0.676 & 0.784 & 0.727 & 0.701 \\
100 & 10 & 0.666 & 0.868 & 0.859 & 0.645 & 0.849 & 0.837 & 0.641 & 0.783 & 0.790 & 0.679 & 0.774 & 0.762 \\
200 & 5 & 0.714 & 0.831 & 0.856 & 0.685 & 0.839 & 0.842 & 0.722 & 0.750 & 0.669 & 0.770 & 0.731 & 0.612 \\
200 & 10 & 0.661 & 0.886 & 0.850 & 0.658 & 0.865 & 0.835 & 0.637 & 0.777 & 0.771 & 0.642 & 0.751 & 0.758 \\ \hline
25 & 50 & 0.662 & 0.850 & 0.859 & 0.663 & 0.880 & 0.884 & 0.613 & 0.847 & 0.843 & 0.653 & 0.848 & 0.846 \\
25 & 100 & 0.674 & 0.888 & 0.886 & 0.676 & 0.883 & 0.886 & 0.662 & 0.877 & 0.879 & 0.651 & 0.864 & 0.872 \\
50 & 50 & 0.633 & 0.863 & 0.869 & 0.685 & 0.896 & 0.899 & 0.661 & 0.852 & 0.856 & 0.672 & 0.866 & 0.868 \\
50 & 100 & 0.646 & 0.861 & 0.874 & 0.685 & 0.887 & 0.887 & 0.670 & 0.863 & 0.873 & 0.668 & 0.856 & 0.860 \\ \hline
\end{tabular}
\vspace{3mm}
\end{center}
\caption{\emph{Empirical coverage probabilities of the bootstrap confidence interval for a nominal 90\% level. A location-scale shift model is considered. CS denotes cross-sectional pairs bootstrap, WB denotes wild bootstrap, PQR denotes the penalized estimator, and FE is the unpenalized fixed effects estimator.}}
\label{mc.table2}
\end{table}
\end{singlespace}
\begin{singlespace}
\begin{table}
\begin{center}\footnotesize
\begin{tabular}{c c c c c c c c c c c c c c} \hline
& & \multicolumn{6}{c}{Quantile 0.5} & \multicolumn{6}{c}{Quantile 0.75} \\
& & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} \\ \cmidrule(lr){3-5} \cmidrule(lr){6-8} \cmidrule(lr){9-11} \cmidrule(lr){12-14}
$N$ & $T$ & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} \\ \cmidrule(lr){4-5} \cmidrule(lr){7-8} \cmidrule(lr){10-11} \cmidrule(lr){13-14}
& & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim \mathcal{N}(0,1)$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.951 & 0.908 & 0.910 & 0.862 & 0.922 & 0.919 & 0.930 & 0.896 & 0.894 & 0.845 & 0.899 & 0.899 \\
100 & 10 & 0.971 & 0.917 & 0.862 & 0.864 & 0.913 & 0.870 & 0.969 & 0.909 & 0.913 & 0.861 & 0.910 & 0.922 \\
200 & 5 & 0.953 & 0.925 & 0.927 & 0.849 & 0.919 & 0.922 & 0.918 & 0.902 & 0.902 & 0.831 & 0.912 & 0.913 \\
200 & 10 & 0.974 & 0.933 & 0.853 & 0.868 & 0.925 & 0.868 & 0.975 & 0.910 & 0.920 & 0.876 & 0.904 & 0.924 \\ \hline
25 & 50 & 0.998 & 0.911 & 0.909 & 0.916 & 0.894 & 0.894 & 1.000 & 0.889 & 0.892 & 0.912 & 0.886 & 0.886 \\
25 & 100 & 1.000 & 0.910 & 0.910 & 0.972 & 0.908 & 0.907 & 1.000 & 0.901 & 0.904 & 0.944 & 0.902 & 0.904 \\
50 & 50 & 1.000 & 0.907 & 0.905 & 0.922 & 0.888 & 0.885 & 1.000 & 0.895 & 0.897 & 0.918 & 0.907 & 0.908 \\
50 & 100 & 1.000 & 0.908 & 0.904 & 0.973 & 0.904 & 0.903 & 1.000 & 0.900 & 0.900 & 0.965 & 0.913 & 0.907 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim t_3$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.945 & 0.922 & 0.918 & 0.860 & 0.920 & 0.916 & 0.917 & 0.887 & 0.885 & 0.814 & 0.865 & 0.901 \\
100 & 10 & 0.969 & 0.923 & 0.873 & 0.856 & 0.923 & 0.878 & 0.964 & 0.908 & 0.916 & 0.857 & 0.894 & 0.907 \\
200 & 5 & 0.946 & 0.940 & 0.939 & 0.826 & 0.913 & 0.931 & 0.911 & 0.898 & 0.901 & 0.774 & 0.841 & 0.903 \\
200 & 10 & 0.967 & 0.920 & 0.822 & 0.850 & 0.933 & 0.838 & 0.960 & 0.901 & 0.915 & 0.848 & 0.880 & 0.897 \\ \hline
25 & 50 & 0.997 & 0.926 & 0.921 & 0.924 & 0.912 & 0.906 & 0.989 & 0.890 & 0.889 & 0.901 & 0.890 & 0.889 \\
25 & 100 & 1.000 & 0.895 & 0.893 & 0.950 & 0.893 & 0.892 & 0.999 & 0.879 & 0.880 & 0.935 & 0.902 & 0.902 \\
50 & 50 & 0.998 & 0.904 & 0.899 & 0.900 & 0.887 & 0.881 & 0.999 & 0.895 & 0.897 & 0.910 & 0.906 & 0.909 \\
50 & 100 & 1.000 & 0.896 & 0.891 & 0.955 & 0.902 & 0.898 & 0.998 & 0.890 & 0.890 & 0.917 & 0.875 & 0.874 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location shift model ($\zeta=0$) and $u \sim \chi_3^2$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.897 & 0.917 & 0.918 & 0.832 & 0.892 & 0.916 & 0.857 & 0.805 & 0.824 & 0.791 & 0.750 & 0.800 \\
100 & 10 & 0.905 & 0.898 & 0.873 & 0.830 & 0.879 & 0.844 & 0.886 & 0.835 & 0.852 & 0.852 & 0.836 & 0.848 \\
200 & 5 & 0.897 & 0.914 & 0.922 & 0.824 & 0.893 & 0.917 & 0.828 & 0.784 & 0.820 & 0.798 & 0.784 & 0.820 \\
200 & 10 & 0.912 & 0.904 & 0.845 & 0.823 & 0.877 & 0.831 & 0.884 & 0.834 & 0.858 & 0.849 & 0.831 & 0.846 \\ \hline
25 & 50 & 0.968 & 0.881 & 0.883 & 0.864 & 0.894 & 0.896 & 0.904 & 0.851 & 0.852 & 0.846 & 0.868 & 0.870 \\
25 & 100 & 0.994 & 0.889 & 0.890 & 0.870 & 0.878 & 0.879 & 0.956 & 0.878 & 0.885 & 0.831 & 0.856 & 0.860 \\
50 & 50 & 0.965 & 0.878 & 0.874 & 0.871 & 0.895 & 0.894 & 0.916 & 0.865 & 0.864 & 0.859 & 0.876 & 0.878 \\
50 & 100 & 0.992 & 0.882 & 0.882 & 0.890 & 0.894 & 0.894 & 0.965 & 0.870 & 0.870 & 0.853 & 0.868 & 0.868 \\ \hline
\end{tabular}
\vspace{3mm}
\end{center}
\caption{\emph{Empirical coverage probabilities of the asymptotic Gaussian confidence interval for a nominal 90\% level. A location shift model is considered. CS denotes cross-sectional pairs bootstrap, WB denotes wild bootstrap, PQR denotes the penalized estimator, and FE is the unpenalized fixed effects estimator.}}
\label{mc.table3}
\end{table}
\end{singlespace}
\begin{singlespace}
\begin{table}
\begin{center}\footnotesize
\begin{tabular}{c c c c c c c c c c c c c c} \hline
& & \multicolumn{6}{c}{Quantile 0.5} & \multicolumn{6}{c}{Quantile 0.75} \\
& & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} & \multicolumn{3}{c}{Method:} \\ \cmidrule(lr){3-5} \cmidrule(lr){6-8} \cmidrule(lr){9-11} \cmidrule(lr){12-14}
$N$ & $T$ & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} & CS & \multicolumn{2}{c}{WB} \\ \cmidrule(lr){4-5} \cmidrule(lr){7-8} \cmidrule(lr){10-11} \cmidrule(lr){13-14}
& & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE & PQR & PQR & FE \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale shift model ($\zeta=0.5$) and $u \sim \mathcal{N}(0,1)$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.851 & 0.893 & 0.886 & 0.863 & 0.909 & 0.905 & 0.812 & 0.841 & 0.837 & 0.830 & 0.876 & 0.843 \\
100 & 10 & 0.833 & 0.912 & 0.865 & 0.821 & 0.910 & 0.860 & 0.819 & 0.873 & 0.871 & 0.825 & 0.887 & 0.877 \\
200 & 5 & 0.838 & 0.899 & 0.896 & 0.834 & 0.906 & 0.895 & 0.784 & 0.838 & 0.831 & 0.810 & 0.884 & 0.837 \\
200 & 10 & 0.833 & 0.917 & 0.845 & 0.817 & 0.913 & 0.847 & 0.812 & 0.873 & 0.870 & 0.801 & 0.884 & 0.879 \\ \hline
25 & 50 & 0.897 & 0.903 & 0.903 & 0.837 & 0.909 & 0.905 & 0.883 & 0.894 & 0.893 & 0.813 & 0.873 & 0.874 \\
25 & 100 & 0.938 & 0.912 & 0.909 & 0.827 & 0.902 & 0.901 & 0.922 & 0.901 & 0.900 & 0.824 & 0.898 & 0.898 \\
50 & 50 & 0.887 & 0.901 & 0.897 & 0.801 & 0.896 & 0.886 & 0.882 & 0.890 & 0.891 & 0.821 & 0.888 & 0.890 \\
50 & 100 & 0.943 & 0.909 & 0.906 & 0.824 & 0.904 & 0.900 & 0.919 & 0.891 & 0.888 & 0.838 & 0.907 & 0.907 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale shift model ($\zeta=0.5$) and $u \sim t_3$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.864 & 0.915 & 0.898 & 0.882 & 0.918 & 0.894 & 0.816 & 0.816 & 0.789 & 0.857 & 0.871 & 0.793 \\
100 & 10 & 0.845 & 0.912 & 0.864 & 0.845 & 0.910 & 0.870 & 0.847 & 0.897 & 0.873 & 0.832 & 0.874 & 0.854 \\
200 & 5 & 0.884 & 0.919 & 0.921 & 0.855 & 0.910 & 0.916 & 0.778 & 0.803 & 0.756 & 0.836 & 0.860 & 0.792 \\
200 & 10 & 0.830 & 0.913 & 0.842 & 0.821 & 0.919 & 0.827 & 0.786 & 0.842 & 0.817 & 0.793 & 0.849 & 0.817 \\ \hline
25 & 50 & 0.894 & 0.919 & 0.912 & 0.851 & 0.916 & 0.914 & 0.856 & 0.877 & 0.878 & 0.814 & 0.883 & 0.881 \\
25 & 100 & 0.902 & 0.887 & 0.883 & 0.826 & 0.903 & 0.903 & 0.874 & 0.886 & 0.886 & 0.831 & 0.890 & 0.893 \\
50 & 50 & 0.882 & 0.894 & 0.889 & 0.791 & 0.885 & 0.873 & 0.856 & 0.882 & 0.884 & 0.820 & 0.899 & 0.897 \\
50 & 100 & 0.908 & 0.886 & 0.882 & 0.815 & 0.889 & 0.885 & 0.893 & 0.892 & 0.892 & 0.805 & 0.879 & 0.878 \\ \hline
\multicolumn{2}{c}{} & \multicolumn{12}{c}{Location-scale model ($\zeta=0.5$) and $u \sim \chi_3^2$} \\
& & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} & \multicolumn{3}{c}{$\alpha_i \sim \mathcal{N}(0,1)$} & \multicolumn{3}{c}{$\alpha_i = i/N$} \\ \hline
100 & 5 & 0.857 & 0.887 & 0.868 & 0.829 & 0.857 & 0.847 & 0.881 & 0.818 & 0.672 & 0.896 & 0.843 & 0.714 \\
100 & 10 & 0.837 & 0.886 & 0.847 & 0.822 & 0.871 & 0.832 & 0.821 & 0.830 & 0.786 & 0.822 & 0.830 & 0.779 \\
200 & 5 & 0.833 & 0.881 & 0.858 & 0.823 & 0.860 & 0.844 & 0.875 & 0.850 & 0.641 & 0.884 & 0.842 & 0.584 \\
200 & 10 & 0.845 & 0.903 & 0.835 & 0.814 & 0.883 & 0.821 & 0.815 & 0.823 & 0.778 & 0.799 & 0.818 & 0.758 \\ \hline
25 & 50 & 0.803 & 0.858 & 0.858 & 0.815 & 0.890 & 0.886 & 0.804 & 0.851 & 0.850 & 0.806 & 0.859 & 0.857 \\
25 & 100 & 0.847 & 0.893 & 0.900 & 0.819 & 0.885 & 0.885 & 0.829 & 0.881 & 0.884 & 0.808 & 0.863 & 0.868 \\
50 & 50 & 0.819 & 0.874 & 0.872 & 0.836 & 0.907 & 0.907 & 0.808 & 0.859 & 0.861 & 0.826 & 0.878 & 0.876 \\
50 & 100 & 0.827 & 0.876 & 0.877 & 0.829 & 0.895 & 0.892 & 0.825 & 0.875 & 0.878 & 0.811 & 0.866 & 0.868 \\ \hline
\end{tabular}
\vspace{3mm}
\end{center}
\caption{\emph{Empirical coverage probabilities of the asymptotic Gaussian confidence interval for a nominal 90\% level. A location-scale shift model is considered. CS denotes cross-sectional pairs bootstrap, WB denotes wild bootstrap, PQR denotes the penalized estimator, and FE is the unpenalized fixed effects estimator.}}
\label{mc.table4}
\end{table}
\end{singlespace}
Tables \ref{mc.table1}, \ref{mc.table2}, \ref{mc.table3}, and \ref{mc.table4} present coverage probabilities for a nominal 90\% confidence interval for the slope parameter $\beta_{0}$. We present coverage probabilities using the empirical distribution of the bootstrap estimator (Tables \ref{mc.table1} and \ref{mc.table2}), as well as coverage probabilities of the asymptotic Gaussian confidence interval (Tables \ref{mc.table3} and \ref{mc.table4}). In the latter case, the coverage is constructed using the standard error of the corresponding bootstrap procedure. Tables \ref{mc.table1} and \ref{mc.table3} present results for the location shift model ($\zeta=0$), while Tables \ref{mc.table2} and \ref{mc.table4} present results for the location-scale shift model ($\zeta=0.5$). The tables present results for $\tau \in
\{0.50,0.75\}$, based on different combinations of $N\in \{25,50,100,200\}$ and $T\in \{5,10,50,100\}$. The number of
bootstrap repetitions is set to 400, and the results are obtained by using 1000 random samples.
The tables show results for two bootstrap methods. The cross-sectional pairs bootstrap (CS) samples over $i$ with replacement, keeping the entire block of
time series observations. The
wild bootstrap (WB) is implemented as discussed in Section \ref{subsec:wild}. We first obtain residuals $\hat{u}_{it}$ using the penalized quantile regression \eqref{penboot}, which is labeled `PQR' in the tables. The tuning parameter is obtained as $\hat{\lambda}_T = b_T \tilde{\lambda}$ where $\tilde{\lambda}$ is obtained by cross-validation and $b_T=0.5 T^{-\nu}$ controls the bias. The selection of $\nu = 1$ performed well in the simulations and it is consistent with Theorem \ref{thm:consistent}. As in the case of the wild bootstrap estimator proposed by Feng, He, and Hu (2011),
a finite sample correction is recommended. We adopt an adjustment following closely the {\tt R} package {\tt quantreg} by Koenker
(2021)\nocite{rK21}. In our case, we adjust the residuals with the influence function and sign function following the Bahadur representation of the estimator derived in Theorem \ref{thm:AN}. Then, we generate $u^\ast_{it} = w_{it} | \hat{u}_{it}
|$, where $w_{it}$ is an i.i.d. random variable distributed as a two-point
distribution with probabilities $\tau$ and $1-\tau$ at $w_{it} = -2 \tau$ and
$w_{it} = 2 (1-\tau)$. Lastly, we generate the dependent variable as
$y_{it}^\ast = \hat{\alpha}_i + \hat{\beta} x_{it} + u^\ast_{it}$. The performance of the estimator \eqref{penbootthr} was similar and the results are not presented here to save space. Finally, we include the estimator \eqref{penboot} defined for $\lambda_T = 0$ and it is labeled `FE'.
Following the result presented in Theorem \ref{thm:boot}, the coverage probabilities in Table \ref{mc.table1} are obtained considering the quantiles of the empirical distribution of $\sqrt{NT} (\beta^\ast - \hat{\beta})$. As can be seen in the upper block of Table \ref{mc.table1}, the performance of the WB bootstrap estimators are excellent, and they are in general around the specified coverage probability. Furthermore, performance improves with $T$, and tends to be similar for both $0.5$ and $0.75$ quantiles. On the other hand, the performance of the CS estimator is poor, with estimates not approaching to specified nominal values. In the lower parts of the table, we present the performance of the estimators for different distributions $F_u$. The WB method continues to perform better than CS, and, as expected, the estimation of the higher quantile is more challenging in the $\chi_3^2$ case. In all the variations of the model considered in the table, the WB estimator performs much better than the CS estimator.
The results for the location-scale shift model presented in Table \ref{mc.table2} are similar. We continue to see that the WB bootstrap performs better than the CS method. This conclusion holds when we consider asymptotic Gaussian confidence intervals obtained using bootstrap standard errors $\mbox{se}(\beta^\ast)$ (see Tables \ref{mc.table3} and \ref{mc.table4}). Moreover, the tables confirm two results that were expected. First, as $T$ increases relative to $N$, the coverage of the WB improves. Second, the performance of WB in the case of $\lambda_T=0$ reveals that, in general, the procedure proposed in this paper is valid for approximating the distribution of the fixed effects estimator.
\begin{figure}
\begin{center}
\centerline{\includegraphics[width=.75\textwidth]{figure41.new.pdf}}
\caption{\emph{The performance of the bootstrap estimators as $\lambda_T$ increases. SD denotes standard deviation of the penalized estimator, CS denotes cross-sectional pair bootstrap, and WB denotes wild bootstrap estimator \eqref{penboot}.}}
\label{mc.figure0}
\end{center}
\end{figure}
We finish the section by briefly documenting the relative performance of the estimators of the standard errors. We generate data from a location-scale shift model ($\zeta=0.5$) when the error term $u_{it} \sim \mathcal{N}(0,1)$ and $\alpha_i \sim \mathcal{N}(0,1)$, by setting $N=100$, $T=10$, and $\tau = 0.5$.
The left panel of Figure \ref{mc.figure0} shows CS and WB bootstrap estimates of the standard error, $\mbox{se}(\beta^\ast)$, and the standard deviation of the penalized estimator, $\mbox{sd}(\hat{\beta})$. The figure shows the advantage of the penalized estimator relative to the fixed effects estimator, as the standard deviation of the estimator is decreasing as $\lambda_T$ increases. We also see that the WB procedure performs better than CS when $\lambda_T$ is relatively small, and the performance of the WB estimator does not seem to change over the degree of shrinkage of the individual effects, as the bias appears to be roughly constant over $\lambda_T$. Using the right panel in Figure \ref{mc.figure0}, we explore further the difference in performance between approaches. The empirical distribution obtained by the CS procedure is not centered at the true value, and the distribution of the standard error of the WB is centered at $\mbox{se}(\hat{\beta}) = 0.081$ (with $\lambda_T = 0.05$).
\section{Extensions}
In this section, we investigate the consistency of the wild bootstrap under different conditions. First, we extend the results of Theorems \ref{thm:consistent} and \ref{thm:AN} to allow for dependent data, and then we focus on the consistency of the wild bootstrap. In such case, we use the following assumptions:
\begin{assumptionC} \label{assume:stationary}
The processes $\{(y_{it},\bm{x}_{it}), t \in 1, 2, \ldots\}$ are strictly stationary for each $i$ and $\beta$-mixing, and independent across $i$. Letting $\{\beta_i(j)\}_j$ denote the $\beta$-mixing coefficients, assume that there are constants $0 < a < 1$ and $B > 0$ such that $\sup_i \beta_i(j) \leq B a^j$ for all $j \geq 1$.
\end{assumptionC}
\begin{assumptionC} \label{assume:density}
The random vector $(u_{it}, u_{it+j})$ has a density conditional on $(\bm{x}_{it}, \bm{x}_{it+j})$ that is bounded uniformly over $i$ and $j \geq 1$.
\end{assumptionC}
\begin{assumptionC} \label{assume:Avar_dependent}
Assume that the matrix $\bm{D}_N$ as defined in Assumption~\ref{assume:Avar} exists and is positive definite for all $N$ under Assumptions \ref{assume:stationary} and~\ref{assume:density} and that $\bm{D} = \lim_{N \rightarrow \infty} \bm{D}_N$ exists and is positive definite. Also assume that
\begin{equation*}
\tilde{\bm{V}} = \lim_{N,T \to \infty} \frac{1}{NT} \sum_{i=1}^N \operatorname{Var} \left( \sum_{t=1}^T (\tau - I(y_{it} < \bm{x}_{it}'\bm{\beta}_0 + \alpha_{i0})) \left( \bm{x}_{it} - \varphi_i^{-1} \bm{E}_i \right) \right)
\end{equation*}
is positive definite.
\end{assumptionC}
Theorem~\ref{thm:AN_dependent} presents both consistency and asymptotic normality results for the estimator with dependent error terms.
\begin{theorem} \label{thm:AN_dependent}
Under Assumptions \ref{assume:stationary}-\ref{assume:Avar_dependent}, \ref{assume:xsupport} and \ref{assume:ID_weakconv}, if $\log(N)^2/T \to 0$ and $\lambda_T = o_p(T)$ as $N, T \to \infty$, the estimator $\hat{\bm{\beta}}$ is consistent. Moreover, if $N^2 (\log N)^3 / T \to 0$ and $\lambda_T = o_p( T^{1/2} (\log N)^{1/2} )$ as $N, T \rightarrow \infty$, then
\begin{equation*}
\sqrt{NT} (\hat{\bm{\beta}} - \bm{\beta}_0) \overset{d}{\longrightarrow} \mathcal{N}(\bm{0}, \tilde{\bm{\Omega}}),
\end{equation*}
where $\tilde{\bm{\Omega}} = \bm{D}^{-1} \tilde{\bm{V}} \bm{D}^{-1}$.
\end{theorem}
Theorem~\ref{thm:boot_dependent} shows consistency of the bootstrap distribution estimator in the case of dependent errors. This more complex situation requires another assumption:
\begin{assumptionA} \label{assume:weight_dependent}
Suppose that
\begin{equation*} \label{joint_G_condition}
\lim_{N, T \rightarrow \infty} \frac{1}{N} \sum_{i=1}^N \sum_{j=1}^{T-1} \left( 1 - \frac{j}{T} \right) \left( \text{P}^*\{w_{it} < 0, w_{it+j} < 0\} - \textnormal{P} \left\{ u_{it} \leq 0, u_{it+j} \leq 0 | \bm{x}_{it}, \bm{x}_{it+j} \right\} \right) = 0.
\end{equation*}
\end{assumptionA}
Assumption~\ref{assume:weight_dependent} is a high-level assumption on the distribution of bootstrap weights. The assumption guarantees that the variance of the bootstrap estimator is bounded and sufficiently close to the true variance, because the weights mimic the within-unit dependence structure of the errors. A feasible version could use a plug-in estimate of the average of the joint conditional CDFs of $(u_{it}, u_{it+j})$ to generate weights that satisfy the average probability.
\begin{theorem} \label{thm:boot_dependent}
Suppose that the bootstrap weights satisfies assumptions~\ref{C1}-\ref{assume:weight_dependent} and the data satisfy assumptions \ref{assume:stationary}-\ref{assume:Avar_dependent}, \ref{assume:xsupport} and \ref{assume:ID_weakconv}. If $N^2 (\log N)^3 / T \to 0$ and $\lambda_T = o_p( T^{1/2} (\log N)^{1/2})$ as $N, T \rightarrow \infty$, then
\begin{equation*}
\sup_{b \in \mathbb{R}^p} \left| \textnormal{P} \left\{ \sqrt{NT}( \bm{\beta}^* - \hat{\bm{\beta}}) \leq b | \bm{S} \right\} - \textnormal{P} \left\{ \sqrt{NT}( \hat{\bm{\beta}} - \bm{\beta}_0) \leq b \right\} \right| \overset{p}{\longrightarrow} 0,
\end{equation*}
where $\bm{S}$ denotes the observed sample and $\bm{\beta}^\ast$ denotes the slope estimator defined by~\eqref{penboot}.
\end{theorem}
Finally, we investigate if the conditions on the size of $T$ relative to $N$ needed for the asymptotic normality in Theorem~\ref{thm:AN} can be improved, especially in the light of recent work by~\citet{AGalvao2020}. If instead of focusing on the stochastic order of the terms of the Bahadur representation of the penalized estimator, we focus on the expected values of the remainder terms, it is possible to show that the rates can be improved substantially. In order to show asymptotic normality, we employ the following assumption about the behavior of the penalty parameter.
\begin{assumptionB} \label{assume:lambda_tail}
For some $\kappa \geq 2$, there exists a constant $K > 0$ such that $\textnormal{P} \left\{ \lambda_T > K T^{1/2} (\log T)^{1/2} \right\} = O(T^{-\kappa})$.
\end{assumptionB}
Assumption~\ref{assume:lambda_tail} dictates the rate at which the probability of observing large a $\lambda_T$ becomes small asymptotically. As illustrated in remark \ref{lambdaT_bound}, it is needed to provide a tail bound for the distribution of individual effects, which figure in the remainder terms of the Bahadur representation used to find the asymptotic distribution of $\bm{\hat{\beta}}$ (such a bound holds naturally for terms related to minimizing the quantile regression objective function with bounded regressors, a fact used extensively in~\citet{AGalvao2020}). In the theorem below, we require $\lambda_T = O_p(\log T) = o_p(T^{1/2} (\log T)^{1/2})$, so this assumption only mildly strengthens the other regularity conditions.
\begin{theorem} \label{thm:rates}
Under Assumptions \ref{assume:data} and \ref{assume:xsupport}-\ref{assume:lambda_tail}, if $N (\log T)^2 / T \to 0$ and $\lambda_T = O_p(\log T)$ as $N, T \rightarrow \infty$, then
\begin{equation*}
\sqrt{NT} (\hat{\bm{\beta}} - \bm{\beta}_0) \overset{d}{\longrightarrow} \mathcal{N}(\bm{0}, \bm{\Omega}),
\end{equation*}
where $\bm{\Omega} = \bm{D}^{-1} \bm{V} \bm{D}^{-1}$.
\end{theorem}
The proof in Theorem \ref{thm:rates} uses an infeasible estimator $\tilde{\alpha}_i$ that is obtained considering $T$ observations $y_{it} - \bm{x}_{it}'\bm{\beta}_0$. The difference between $\tilde{\alpha}_i$ and $\hat{\alpha}_i$ converges to zero as the slope coefficient $\hat{\bm{\beta}}$ converges in probability towards $\bm{\beta}_0$, under the condition on $\lambda_T$. Therefore, the remainder terms of the corresponding Bahadur representations are sufficiently close, leading to the improvements in the rates first obtained in \citet{AGalvao2020} for the fixed effects estimator. We now show the consistency of the bootstrap distribution estimator under these relatively closer orders of $N$ and $T$.
\begin{theorem} \label{thm:boot_rates}
Under Assumptions \ref{C1}-\ref{C3} and the conditions of Theorem \ref{thm:rates},
\begin{equation*}
\sup_{b \in \mathbb{R}^p} \left| \textnormal{P} \left\{ \sqrt{NT}( \bm{\beta}^* - \hat{\bm{\beta}}) \leq b | \bm{S} \right\} - \textnormal{P} \left\{ \sqrt{NT}( \hat{\bm{\beta}} - \bm{\beta}_0) \leq b \right\} \right| \overset{p}{\longrightarrow} 0.
\end{equation*}
where $\bm{S}$ denotes the observed sample and $\bm{\beta}^\ast$ denotes the slope estimator defined by~\eqref{penboot}.
\end{theorem}
\section{An Empirical Illustration}
In recent years, policy makers and the general public have been debating and re-evaluating several aspects of trade, including the benefits of trade agreements \citep[][among others]{mB2001, sH2016}. An important question is whether workers have been negatively affected by the North American Free Trade Agreement (NAFTA), which was signed by the governments of the United States of America, Canada, and Mexico in 1993. Hakobyan and McLaren (2016) find that the effect of NAFTA on \textit{average} wage growth in the period 1990-2000 was negative. In this section, we use similar data and apply our approach to study the distributional impact of NAFTA. Our findings suggest that the agreement increased wage inequality. Low-wage workers experienced significant negative wage growth, while high-wage workers experienced, in general, significant positive wage growth. Our results are similar to evidence on the effect of Chinese imports on low-wage American workers \citep{denisChet2016}.
\subsection{Data}
Following Hakobyan and McLaren (2016), we use a 5\% sample from the U.S. Census. We employ two cross-sectional samples in the year 1990 and 2000, and therefore, workers in the sample are observed once. The longitudinal nature of the analysis comes from exploiting the fact that we observe multiple individuals in a given industry and location. The sample includes workers between 25 and 64 years of age who reported positive income. We have demographic information including age, gender, marital status, race, and educational attainment of the worker classified in four categories: high school dropout, high school graduate, some college, and college graduate.
The data on U.S. tariffs and Mexico's revealed comparative advantage (RCA) are obtained from Hakobyan and McLaren (2016). Using their data, we have access to average U.S. tariffs by industry of employment of the worker and location (or Consistent Public-Use Microdata Area, abbreviated \emph{conspuma}) of residence of the worker. In 1990, the average tariff by industry in 1990 was 2.1\% percent (with a standard deviation of 3.9\%), while the average local tariff by conspuma level was 1.03\% (with a standard deviation of 0.67\%). In the period 1990-2000, the tariffs decreased 1.7\% at the industry level and 0.9\% at the conspuma level. These descriptive statistics are used in the next section to estimate the percentage change in wages associated with the reduction in tariffs. We consider all industries with the exception of agriculture.
\subsection{Model}
To investigate the effect of NAFTA on the wages of American workers, we consider a specification that allows for the impact of the trade agreement to vary by industry, location, and educational attainment of the worker. To that end, we consider the following model as in Hakobyan and McLaren (2016):
\begin{equation}
y_{ijc} = \bm{\beta}_{1L}' \bm{L}_{ic} + \bm{\beta}_{2L}' \Delta \bm{L}_{ic} + \bm{\beta}_{1I}' \bm{I}_{ij} + \bm{\beta}_{2I}' \Delta \bm{I}_{ij} + \bm{X}_{ijc}' \bm{\Pi} + \alpha_{jc} + u_{ijc}, \label{main}
\end{equation}
where the response variable $y_{ijc}$ is the logarithm of wages for worker $i$, who is employed in industry $j$ and resides in conspuma $c$, $\bm{L}_{ic}$ and $\Delta \bm{L}_{ic}$ are location variables to be described below, $\bm{I}_{ij}$ and $\Delta \bm{I}_{ij}$ are industry variables, $\bm{X}_{ijc}$ is the vector of control variables considered in Hakobyan and McLaren (2016), and $\alpha_{jc}$ is a industry-conspuma effect. The error term is denoted by $u_{ijc}$.
The location variables are defined as $\bm{L}_{ic} = (L_{ic,1},L_{ic,2},L_{ic,3},L_{ic,4})'$, where $L_{ic,k}$ is the product of an indicator for educational category $k$ of worker $i$, an indicator variable for whether $i$ is in the 2000 sample, and the average tariff in the conspuma of residence of worker $i$. Similarly, we can define $\Delta \bm{L}_{ic} = (\Delta L_{ic,1},\Delta L_{ic,2},\Delta L_{ic,3},\Delta L_{ic,4})'$, as the change in $\bm{L}_{ic}$ due to the change in tariffs between 1990 and 2000 in the conspuma of residence of worker $i$. In terms of the industry variables, $\bm{I}_{ij} = (I_{ij,1},I_{ij,2},I_{ij,3},I_{ij,4})'$, where $I_{ij,k}$ is the product of an indicator for educational category $k$, the RCA in industry $j$, an indicator variable for whether $i$ is in the 2000 sample, and the tariff of the industry that employs worker $i$. Similarly, we define $\Delta \bm{I}_{ij} = (\Delta I_{ij,1},\Delta I_{ij,2},\Delta I_{ij,3},\Delta I_{ij,4})'$, as the change in $\bm{I}_{ij,k}$ due to the tariff change between 1990 and 2000 in the industry that employs worker $i$.
Because industry latent factors and trends in some areas can affect wages and also the changes in tariffs, we employ the penalized estimator \eqref{pqr} to estimate a high-dimensional model with more than 84,000 parameters $\alpha_{jc}$. The parameters of interest in equation \eqref{main} are $\bm{\beta}_{1L}$, $\bm{\beta}_{2L}$, $\bm{\beta}_{1I}$, and $\bm{\beta}_{2I}$, which measure the initial effect of tariffs by location and industry ($\bm{\beta}_{1L}$ and $\bm{\beta}_{1I}$), and the impact effect of a reduction of tariffs by location and industry ($\bm{\beta}_{2L}$ and $\bm{\beta}_{2I}$). Using these parameters, it is possible to obtain the effect of the trade agreement on wages. For instance, for locations that lost all of their protection after the introduction of NAFTA, the effect of the local average tariff is measured by $\bm{\beta}_{1L} - \bm{\beta}_{2L}$. Similarly, for industries that lost all of their protection, the effect of the industry tariff is $\bm{\beta}_{1I} - \bm{\beta}_{2I}$.
\begin{singlespace}
\begin{table}
\begin{center}\small
\begin{tabular}{l c c c c c c } \hline
& Mean & \multicolumn{5}{c}{Quantiles} \\
& Effect & 0.1 & 0.25 & 0.5 & 0.75 & 0.9 \\ \hline
& \multicolumn{6}{c}{High school dropouts} \\ \hline
Initial tariff effect, $\beta_{1I,1}$ & 2.018 & 1.156 & 2.603 & 1.880 & 0.991 & 0.434 \\
&(1.274)&(0.820)&(1.120)&(1.047)&(0.847)&(1.105)\\
Impact effect, $\beta_{2I,1}$ &3.569&3.082&4.625&3.245&1.666&0.600\\
&(1.544)&(0.945)&(1.314)&(1.191)&(1.024)&(1.290)\\
Industry effect: $\beta_{1I,1}-\beta_{2I,1}$&-1.551&-1.925&-2.022&-1.365&-0.675&-0.166\\
&[0.000]&[0.000]&[0.000]&[0.000]&[0.000]&[0.556]\\ \hline
& \multicolumn{6}{c}{High school graduates} \\ \hline
Initial tariff effect, $\beta_{1I,2}$ &1.081&5.015&2.224&0.426&-2.216&-2.933\\
&(0.870)&(0.523)&(0.626)&(0.747)&(0.515)&(0.436)\\
Impact effect, $\beta_{2I,2}$ &2.315&9.259&4.318&1.337&-2.469&-3.855\\
&(1.086)&(0.595)&(0.736)&(0.873)&(0.618)&(0.543)\\
Industry effect: $\beta_{1I,2}-\beta_{2I,2}$&-1.234&-4.245&-2.094&-0.911&0.253&0.922\\
&[0.000]&[0.000]&[0.000]&[0.000]&[0.022]&[0.000]\\ \hline
& \multicolumn{6}{c}{Some college} \\ \hline
Initial tariff effect, $\beta_{1I,3}$ &-0.181&3.187&2.631&-0.921&-2.963&-3.765\\
&(1.146)&(0.820)&(1.172)&(1.151)&(0.779)&(0.879)\\
Impact effect, $\beta_{2I,3}$&1.070&7.360&4.889&-0.263&-3.452&-4.662\\
&(1.396)&(0.972)&(1.468)&(1.359)&(0.954)&(1.026)\\
Industry effect: $\beta_{1I,3}-\beta_{2I,3}$&-1.234&-4.245&-2.094&-0.911&0.253&0.922\\
&[0.000]&[0.000]&[0.000]&[0.000]&[0.022]&[0.000]\\ \hline
& \multicolumn{6}{c}{College graduate} \\ \hline
Initial tariff effect, $\beta_{1I,4}$&-2.438&7.623&-1.363&-6.538&-7.681&-8.688\\
&(1.839)&(1.826)&(1.362)&(1.856)&(1.041)&(1.181)\\
Impact effect, $\beta_{2I,4}$&-2.095&12.840&-0.024&-8.066&-9.828&-11.490\\
&(2.175)&(2.215)&(1.630)&(2.291)&(1.178)&(1.301)\\
Industry effect: $\beta_{1I,4}-\beta_{2I,4}$&-0.343&-5.217&-1.339&1.528&2.147&2.801\\
&[0.439]&[0.000]&[0.000]&[0.000]&[0.000]&[0.000]\\ \hline
Location variables & Yes & Yes & Yes & Yes & Yes & Yes \\
Control variables & Yes & Yes & Yes & Yes & Yes & Yes \\
Number of $\alpha_{jc}$ effects & 84,266 & 84,266 & 84,266 & 84,266 & 84,266 & 84,266 \\
Observations & 9,580,568 & 9,580,568 & 9,580,568 & 9,580,568 & 9,580,568 & 9,580,568 \\ \hline
\end{tabular}
\vspace{3mm}
\end{center}
\caption{\emph{Regression results for the industry effects by educational category of the worker. We present standard errors in parenthesis, and p-values of a test for the equality of initial and impact effects in brackets.}}
\label{table1.results}
\end{table}
\end{singlespace}
\subsection{Main empirical results}
Table \ref{table1.results} reports results for the coefficients $\bm{\beta}_{1I}$,
and $\bm{\beta}_{2I}$ for the four educational categories. The table also shows results for $\beta_{1I,k} - \beta_{2I,k}$ for each educational category $k$ and p-values (in brackets) of Wald-type tests for the null hypothesis $\mbox{H}_0: \beta_{1I,k} = \beta_{2I,k}$. The variance of the test is obtained using the proposed wild residual bootstrap procedure. The first column presents mean fixed effects regression results, that is, estimation of model \eqref{main} by least squares methods. The last five columns show penalized quantile regression (PQR) results with $\lambda_T$ selected by cross-validation. The standard errors are obtained by the proposed wild residual bootstrap procedure. To save space, we do not present results on the control variables included in the vector $\bm{L}_{ic}$, $\Delta \bm{L}_{ic}$, and $\bm{X}_{ijc}$, but the fixed effects results shown in the first column are similar to the results in Table 4 (column (2)) in Hakobyan and McLaren (2016).
Looking at the first set of estimates in the first rows, we see that an initial tariff estimate equal to 2.02 and an impact effect of 3.57. Based on the standard deviation of tariffs at the industry level, a 1\% standard deviation increase in the initial industry tariff has an effect of reducing wages by $3.9\% \times -1.55$, or $-6.05\%$ in the period 1990-2000. This implies that, among industries with tariff declining after the introduction of NAFTA, average wage growth is negative for high school dropouts. The results, however, show that the average response does not summarize well the distributional impact of NAFTA. While the industry effect, which is measured as the difference between the initial effect and the impact effect, is negative ($-1.93$, or $-7.50\%$) and significant for high school dropouts at the 0.1 quantile, it is small ($-0.17$, or $-0.65\%$) and insignificant at the 0.9 quantile. Moreover, we find that the largest differences between the 0.1 and 0.9 effects are among college graduates in industries that lost all of their protection, suggesting that wage growth has been also unequal by educational attainment.
\begin{singlespace}
\begin{figure}
\begin{center}
\centerline{\includegraphics[width=.6\textwidth]{nafta-figure1r.pdf}}
\caption{\emph{Conditional wage growth impacts. PQR denotes penalized quantile regression and the dashed areas are 95\% confidence intervals.}\label{nafta.figure}}
\end{center}
\end{figure}
\end{singlespace}
Lastly, using Figure \ref{nafta.figure}, we report point estimates and confidence intervals for the location and industry effects for high school dropouts and college graduates.
The evidence reveals that inequality increased in the period after the implementation of the trade agreement.
\section{Conclusion}
In this article, we address the problem of estimating the distribution of the penalized quantile regression estimator for longitudinal data using a wild residual bootstrap procedure. Originally introduced by Koenker (2004) as a convenient alternative to the quantile regression estimator with fixed effects, the practical use of the penalized estimator has been limited by challenges involving inference. We show that the wild bootstrap procedure is asymptotically valid for approximating the distribution of the penalized estimator. We derive a series of new asymptotic results and carry out a simulation study that indicates that the wild residual bootstrap performs better than an alternative bootstrap approach commonly used in practice for similar estimators that do not include a penalty term.
Although the paper makes an important contribution by providing a valid method for statistical inference, there are several questions that remain to be answered. We believe that the procedure leads to valid inference in the case of $J$ quantiles estimated simultaneously, but we leave this to future research. Moreover, under an assumption of sparsity as in other high-dimensional models, we expect changes in the consistency and asymptotic normality results. In terms of theoretical developments, we did not consider the case where $\alpha_i$ is a random effect. Lastly, the practical implementation of the wild bootstrap in the case of dependent data involves a few challenges. We hope to investigate these directions in future work.