EconBase
← Back to paper

Low-rank Panel Quantile Regression: Estimation and Inference

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.

111,590 characters

Low-rank Panel Quantile Regression: Estimation and Inference



\title{Low-rank Panel Quantile Regression: Estimation and Inference\thanks{
Su gratefully acknowledges the support from the National Natural Science
Foundation of China under Grant No. 72133002. Zhang acknowledges the financial support from a Lee Kong Chian fellowship. Any and all errors are our own. }}
\author{Yiren Wang$^{a}$, Liangjun Su$^{b}$ and Yichong Zhang$^{a}$ \\
$^{a}$School of Economics, Singapore Management University, Singapore\\
$^{b}$School of Economics and Management, Tsinghua University, China}
\maketitle

\begin{abstract}
In this paper, we propose a class of low-rank panel quantile regression
models which allow for unobserved slope heterogeneity over both individuals
and time. We estimate the heterogeneous intercept and slope matrices via
nuclear norm regularization followed by sample splitting, row- and
column-wise quantile regressions and debiasing. We show that the estimators
of the factors and factor loadings associated with the intercept and slope
matrices are asymptotically normally distributed. In addition, we develop
two specification tests: one for the null hypothesis that the slope
coefficient is a constant over time and/or individuals under the case that
true rank of slope matrix equals one, and the other for the null hypothesis
that the slope coefficient exhibits an additive structure under the case
that the true rank of slope matrix equals two. We illustrate the finite
sample performance of estimation and inference via Monte Carlo simulations
and real datasets.\medskip

\noindent \textbf{Key words:} Debiasing, heterogeneity, nuclear norm
regularization, panel quantile regression, sample splitting, specification
test. \medskip

\noindent \textbf{JEL Classification:} C23, C31, C32, C52\medskip\

\ \ \ \ \ \ \ \ \ \ \

\ \ \ \ \ \ \ \ \ \ \
\end{abstract}

\ \ \ \ \pagebreak

\section{Introduction}

Panel quantile regressions are widely used to estimate the conditional
quantiles, which can capture the heterogeneous effects that may vary across
the distribution of the outcomes. Such effects are usually assumed to be
homogeneous across individuals and over time periods. However, in empirical
analyses, it is usually unknown whether the slope coefficients are
homogeneous across individuals and/or time. Mistakenly forcing slopes to be
homogeneous across time and individuals may lead to inconsistent estimation
and misleading inferences. This prompts two questions to be answered: how
can we estimate the true model at different quantiles when we allow for
heterogeneous slopes across individuals and time at the same time? How to
conduct specification tests for homogeneous effects over individuals or time
and tests for the additive structure of the slope coefficients?

To answer the first question, we propose an estimation procedure for
heterogeneous panel quantile regression models where we allow the fixed
effects to be either additive or interactive, and the slope coefficients to
be heterogeneous over both individuals and time. We impose a low-rank
structure for both the intercept and slope coefficient matrices and estimate
them via nuclear norm regularization (NNR) followed by the sample splitting,
row- and column-wise quantile regressions and debiasing steps. The
estimation algorithm is inspired by \cite{chernozhukov2019inference}, where
the main difference is that we split the full sample into three subsamples
rather than two because we need certain uniform results which require
independence of regressors and regressand used in the debiasing step, and we
do not have the closed form for the quantile regression estimates. At last,
we derive the asymptotic distributions for the estimators of the factors and
factor loadings associated with slope coefficient matrices.

To answer the second question, under the case when the rank of slope
coefficient matrix equals one, we conduct sup-type specification tests for
homogeneous effects over individuals or time following the lead of \cite
{castagnetti2015inference} and \cite{lu2021uniform}. We show that our
sup-test statistics follow the Gumbel distribution under the null, and the
tests have non-trivial power against certain classes of local alternatives.
Under the case when the rank of slope matrix equals two, our sup-type test
statistic is also shown to follow the Gumbel distribution under the null
that the slope coefficient exhibits an additive structure.

This paper relates to three bunches of literature. First, we contribute to
the large literature on panel quantile regressions (PQRs). Since \cite
{koenker2004quantile} studied the PQRs with individual fixed effects, there
has been an increasing number of papers on PQRs. \cite{galvao2010penalized},
\cite{kato2012asymptotics}, \cite{galvao2015efficient}, \cite
{galvao2016smoothed}, \cite{machado2019quantiles}, and \cite
{galvao2020unbiased} study the asymptotics for PQRs with individual fixed
effects. \cite{chen2021quantile} study quantile factor models and \cite
{chen2019two} considers PQRs with interactive fixed effects (IFEs). We
complement the literature by allowing for unobserved heterogeneity in the
slope coefficients of PQRs.

Second, our paper also pertains to slope heterogeneity in panel data models.
Latent group structures across individuals and structural changes over time
are two common types of slope heterogeneity that have received vast
attention in the literature. To recover the unobserved group structures,
various methods have been proposed. For example, \cite{lin2012estimation},
\cite{bonhomme2015grouped} and \cite{ando2016panel} use the K-means
algorithm; \cite{su2016identifying} propose the C-lasso algorithm which is
further studied and extended by \cite{su2018identifying}, \cite{su2019sieve}
and \cite{wang2019heterogeneous}; \cite{wang2018homogeneity} propose an
clustering algorithm in regression via data-driven segmentation called
CARDS; \cite{wang2021identifying} propose a sequential binary segmentation
algorithm to identify the latent group structures in nonlinear panels.
Recent literature on the estimation with structural changes in panel data
models includes, but is not limited to, \cite{chen2015estimating}, \cite
{cheng2016shrinkage}, \cite{ma2018estimation}, \cite{baltagi2021estimating}.
In addition, \cite{galvao2018testing} and \cite{zhang2019quantile} consider
individual heterogeneity in PQRs while they assume homogeneity across time.
To allow for both latent groups and structural breaks, \cite
{okui2021heterogeneous} study a linear panel data model with individual
fixed effects where each latent group has common breaks and the breaking
points can be different across different groups, and they propose a grouped
adaptive group fused lasso (GAGFL) approach to estimate slope coefficients.
\cite{lumsdaine2021estimation} consider a linear panel data model with a
grouped pattern of heterogeneity where the latent group membership structure
and/or the values of slope coefficients can change at a breaking point, and
they propose a K-means-type estimation algorithm and establish the
asymptotic properties of the resulting estimators. Compared with the models
studied above, our model combines both individual and time heterogeneity and
only requires certain low-rank structure in the slope coefficient matrix. So
the unobserved heterogeneity takes a more flexible form in our model than
those in the literature such as \cite{okui2021heterogeneous} and \cite
{lumsdaine2021estimation}.


Last, our paper also connects with the burgeoning literature on nuclear norm
regularization. Such a method has been widely adopted to study panel and
network models. See, \cite{alidaee2020recovering}, \cite{athey2021matrix},
\cite{bai2019rank}, \cite{belloni2019high}, \cite{chen2020noisy}, \cite
{chernozhukov2019inference}, \cite{feng2019regularized}, \cite
{Hong_Su_Jiang2022}, \cite{Miao_Phillips_Su2022}, among others. In the least
squares panel framework, \cite{moon2018nuclear} consider a homogeneous panel
with IFEs by using NNR-based estimator as an initial estimator to construct
iterative estimators that are asymptotically equivalent to the least squares
estimators; \cite{chernozhukov2019inference} study a heterogenous panel
where both the intercept and slope coefficient matrices exhibit a low-rank
structure and establish the asymptotic distribution theory based on NNR. In
the presence of endogeneity, \cite{Hong_Su_Jiang2022} proposes a profile GMM
method to estimate panel data models with IFEs. In the panel quantile
regression setting, \cite{feng2019regularized} develops error bounds for the
low-rank estimates in terms of Frobenius norms under independence
assumption; \cite{belloni2019high} relaxes the independence assumption to
the $\beta $-mixing condition along the time dimension. Our paper extends
\cite{chernozhukov2019inference} from the least squares framework to the PQR
framework, derives the asymptotic distribution theory and develops various
specification tests under some strong mixing conditions along the time
dimension that is weaker than the $\beta $-mixing condition. We also rely on
the sequential symmetrization technique developed by \cite
{rakhlin2015sequential} to obtain the convergence rates of the nuclear norm
regularized estimators.


The rest of the paper is organized as follows. We first introduce the
low-rank structure PQR model and the estimation algorithm in Section 2. We
study the asymptotic properties of our estimators in Section 3. In Section
4, we propose two specification tests: one for the no-factor structure and
one for the additive structure, and study the asymptotic properties of the
test statistics. In Section 5, we show the finite sample performance of our
method via Monte Carlo simulations. In Section 6, we apply our method to two
datasets: one is to study how Tobin's q and cash flows affect corporate
investment and whether firm's external investment to its internal financing
exhibits heterogeneity structure, and the other is to study the relationship
between economics growth, foreign direct investment and unemployment.
Section 7 concludes. All proofs are related to the online supplement.

\textit{Notation.} $\left\Vert \cdot \right\Vert _{1}$, $\left\Vert \cdot
\right\Vert _{op}$, $\left\Vert \cdot \right\Vert _{\infty }$, $\left\Vert
\cdot \right\Vert _{\max }$ $\left\Vert \cdot \right\Vert _{2}$, $\left\Vert
\cdot \right\Vert _{F}$, $\left\Vert \cdot \right\Vert _{\ast }$ denote the
matrix norm induced by 1-norms, the matrix norm induced by 2-norms, the
matrix norm induced by $\infty $-norms, the maximum norm, the Euclidean
norm, the Frobenius norm and the nuclear norm. $\odot $ is the element-wise
product. $\lfloor \cdot \rfloor $ and $\lceil \cdot \rceil $ denote the
floor and ceiling functions, respectively. $a\vee b$ and $a\wedge b$ return
the max and the min of $a$ and $b,$ respectively. The symbol $\lesssim $
means \textquotedblleft the left is bounded by a positive constant times the
right\textquotedblright . Let $A=\{A_{it}\}_{i \in [n], t\in [T]}$ be a matrix with its $(i,t)$-th
entry denoted as $A_{it}$, where $[n]$ to denote the set $
\{1,\cdots ,n\}$ for any positive integer $n$. Let $\{A_{j}\}_{j=0}^{p}$ denote the collection
of matrices $A_{j}$ for all $j\in \{0,\cdots ,p\}$. When $A$ is symmetric, $
\lambda _{\max }(A)$ and $\lambda _{\min }(A)$ denote its largest and
smallest eigenvalues, respectively. The operators $\rightsquigarrow $ and $
\operatornamewithlimits{\to}\limits^{p}$ denote convergence in distribution
and in probability, respectively. Besides, we use w.p.a.1 and
a.s. to abbreviate \textquotedblleft with probability approaching
1\textquotedblright\ and \textquotedblleft almost surely\textquotedblright ,
respectively.

\section{Model and Estimation}

In this section, we introduce the PQR model and estimation algorithm.

\subsection{Model}

Consider the PQR model
\begin{equation}
\mathscr{Q}_{\tau }\left( Y_{it}\bigg|\left\{ X_{j,it}\right\} _{j\in
\lbrack p],t\in \lbrack T]},\left\{ \Theta _{j,it}^{0}\left( \tau \right)
\right\} _{j\in \lbrack p]\cup \{0\},t\in \lbrack T]}\right) =\Theta
_{0,it}^{0}(\tau )+\sum_{j=1}^{p}X_{j,it}\Theta _{j,it}^{0}(\tau ),
\label{eq:model}
\end{equation}
where $i\in \left[ N\right] ,$ $t\in \left[ T\right] ,$ $\tau \in (0,1)$ is
the quantile index, $Y_{it}$ is the dependent variable, $X_{j,it}$ is the $j$
-th regressor for individual $i$ at time $t$, $\{\Theta _{j,it}^{0}\}_{j\in
\lbrack p]}$ is the corresponding slope coefficient, $\Theta _{0,it}^{0}$ is
the intercept, and $\mathscr{Q}_{\tau }\left( Y_{it}\bigg|\left\{
X_{j,it}\right\} _{j\in \lbrack p],t\in \lbrack T]},\left\{ \Theta
_{j,it}^{0}\left( \tau \right) \right\} _{j\in \lbrack p]\cup \{0\},t\in
\lbrack T]}\right) $ denotes the conditional $\tau $-quantile of $Y_{it}$
given the regressors $\left\{ X_{j,it}\right\} _{j\in \lbrack p],t\in
\lbrack T]}$ and\ the parameters $\left\{ \Theta _{j,it}^{0}\left( \tau
\right) \right\} _{j\in \lbrack p]\cup {0},t\in \lbrack T]}$.\footnote{
We will assume that both the intercept term $\Theta _{0,it}^{0}$ and the
slope coefficients $\{\Theta _{j,it}^{0}\}_{j\in \lbrack p]}$ have low-rank
structures, and follow the convention in the panel data literature by
treating the factors to be random. Therefore, $\{\Theta _{j,it}^{0}\}_{j\in
\lbrack p]\cup \{0\}}$ are random as well.} Alternatively, we can rewrite
the above model as
\begin{align}
& Y=\Theta _{0}^{0}(\tau )+\sum_{j=1}^{p}X_{j}\odot \Theta _{j}^{0}(\tau
)+\epsilon (\tau)\quad \text{and}  \notag \\
& \mathscr{Q}_{\tau }\left( \epsilon _{it} (\tau)\bigg|\left\{ X_{j,it}\right\}
_{j\in \lbrack p],t\in \lbrack T]},\left\{ \Theta _{j,it}^{0}\left( \tau
\right) \right\} _{j\in \lbrack p]\cup \{0\},t\in \lbrack T]}\right) =0,
\end{align}
where $\epsilon (\tau)$ is the idiosyncratic error matrix with the $(i,t)$
-th entry being $\epsilon _{it} (\tau)$. Similarly, $X_{j}$, $\Theta
_{j}\left( \tau \right) $, and $Y$ are matrices with the $(i,t)$-th entry
being $X_{j,it}$, $\Theta _{j,it}\left( \tau \right) $, and $Y_{it}$,
respectively. In this model, we assume $p$, the number of regressors, is
fixed and both $N$ and $T$ pass to infinity. In Assumption \ref{ass:1}
below, we characterize the dependence of the data, under which
\eqref{eq:model} holds.

In the paper, we focus on the panel quantile regression for a fixed $\tau $
and thus suppress the dependence of $\Theta _{j}^{0}(\tau )$ and $\epsilon
(\tau )$ on $\tau $ for notation simplicity. In addition, we impose low-rank
structures for the intercept and slope matrices, i.e., $\text{rank}(\Theta
_{j}^{0})=K_{j}$ for some positive constant $K_{j}$ and for each $j\in
\{0,\cdots ,p\}$. By the singular value decomposition (SVD), we have
\begin{equation*}
\Theta _{j}^{0}=\sqrt{NT}\mathcal{U}_{j}^{0}\Sigma _{j}^{0}\mathcal{V}
_{j}^{0\prime }=U_{j}^{0}V_{j}^{0\prime }\text{ }\,\forall \,j=0,\cdots ,p,
\end{equation*}
where $\mathcal{U}_{j}^{0}\in \mathbb{R}^{N\times K_{j}}$, $\mathcal{V}
_{j}^{0}\in \mathbb{R}^{T\times K_{j}}$, $\Sigma _{j}^{0}=\text{diag}(\sigma
_{1,j},\cdots ,\sigma _{K_{j},j})$, $U_{j}^{0}=\sqrt{N}\mathcal{U}
_{j}^{0}\Sigma _{j}^{0}$ with each row being $u_{i,j}^{0\prime }$, and $
V_{j}^{0}=\sqrt{T}\mathcal{V}_{j}^{0}$ with each row being $v_{t,j}^{0\prime
}$.

The low-rank structure assumption includes several popular cases. For the
intercept term, one commonly assumes that $\Theta _{0,it}^{0}$ to take the
forms $\alpha _{i}^{0},$ $\mu _{t}^{0},$ or $\alpha _{i}^{0}+\mu _{t}^{0}$
in classical PQRs. Then the matrix $\Theta _{0}^{0}$ has rank 1, 1, and 2,
respectively. It is also possible to assume $\Theta _{0,it}^{0}$ to take an
interactive form, say, $\Theta _{0,it}^{0}=\lambda _{0,i}^{0\prime
}f_{0,t}^{0},$ where both $\lambda _{0,i}^{0}$ and $f_{0,t}^{0}$ are $K_{0}$
-vectors. For the slope matrix $\Theta _{j}^{0}$, $j\in \left[ p\right] ,$
the early PQR models frequently assume that $\Theta _{j,it}^{0}$ is a
constant across $\left( i,t\right) \ $ to yield a homogenous PQR model.
Obviously, such a model is very restrictive by assuming homogenous slope
coefficients. It is possible to allow the slope coefficients to change over
either $i$, or $t$, or both. See the following examples for different
low-rank structures.\bigskip

\noindent \textbf{Example 1.} When $\Theta _{j,it}^{0}=\Theta _{j,i}^{0}$ $
\forall t\in \left[ T\right] ,$ or $\Theta _{j,it}^{0}=\Theta _{j,t}^{0}$ $
\forall i\in \left[ N\right] ,$ or $\Theta _{j,it}^{0}=\Theta _{j}^{0}$ $
\forall \left( i\text{,}t\right) \in \left[ N\right] \times \left[ T\right] $
, and this holds for all $j\in \left[ p\right] ,$ we have the PQR models
with only individual heterogeneity, with only time heterogeneity, and with
homogeneity, respectively. We observe that $K_{j}=1$ for these three
cases.\bigskip

\noindent \textbf{Example 2.} When $\Theta _{j,it}^{0}=\lambda
_{j,i}^{0}+f_{j,t}^{0}$, we notice that
\begin{equation*}
\frac{\Theta _{j}^{0}}{\sqrt{NT}}=
\begin{bmatrix}
\frac{1}{\sqrt{N}} & \frac{\lambda _{j,1}^{0}}{\sqrt{N}} \\
\vdots  & \vdots  \\
\frac{1}{\sqrt{N}} & \frac{\lambda _{j,N}^{0}}{\sqrt{N}}
\end{bmatrix}
\begin{bmatrix}
\frac{f_{j,1}^{0}}{\sqrt{T}} & \cdots  & \frac{f_{j,T}^{0}}{\sqrt{T}} \\
\frac{1}{\sqrt{T}} & \cdots  & \frac{1}{\sqrt{T}}
\end{bmatrix}
:=A_{j}B_{j}^{\prime }.
\end{equation*}
Let $\Sigma _{A,j}:=A_{j}^{\prime }A_{j}$ and $\Sigma _{B,j}:=B_{j}^{\prime
}B_{j}$. Let $\Sigma _{A,j}^{\frac{1}{2}}$ (resp. $\Sigma _{B,j}^{\frac{1}{2}
}$) be the symmetric square root of $\Sigma _{A,j}$ (resp. $\Sigma _{B,j}$).
By eigendecomposition, we have $\Sigma _{A,j}^{\frac{1}{2}
}=P_{j,1}S_{j,1}P_{j,1}^{\prime }$ and $\Sigma _{B,j}^{\frac{1}{2}
}=P_{j,2}S_{j,2}P_{j,2}^{\prime }$. Besides, we apply singular value
decomposition to matrix $S_{j,1}P_{j,1}^{\prime }P_{j,2}S_{j,2}$: $
S_{j,1}P_{j,1}^{\prime }P_{j,2}S_{j,2}=Q_{j,1}R_{j}Q_{j,2}^{\prime }$. Then
it follows that
\begin{align*}
\frac{\Theta _{j}^{0}}{\sqrt{NT}}& =A_{j}B_{j}^{\prime }=A_{j}\Sigma
_{A,j}^{-\frac{1}{2}}P_{j,1}S_{j,1}P_{j,1}^{\prime
}P_{j,2}S_{j,2}P_{j,2}^{\prime }\Sigma _{B,j}^{-\frac{1}{2}}B_{j}^{\prime }
\\
& =A_{j}\Sigma _{A,j}^{-\frac{1}{2}}P_{j,1}Q_{j,1}R_{j}Q_{j,2}^{\prime
}P_{j,2}^{\prime }\Sigma _{B,j}^{-\frac{1}{2}}B_{j}^{\prime }:=\mathcal{U}
_{j}^{0}\Sigma _{j}^{0}\mathcal{V}_{j}^{0\prime },
\end{align*}
where $\mathcal{U}_{j}^{0}=A_{j}\Sigma _{A,j}^{-\frac{1}{2}}P_{j,1}Q_{j,1}$,
$\Sigma _{j}^{0}=R_{j}$ and $\mathcal{V}_{j}^{0}=B_{j}\Sigma _{B,j}^{-\frac{1
}{2}}P_{j,2}Q_{j,2}$. Given $P_{j,1}$, $P_{j,2}$, $Q_{j,1}$ and $Q_{j,2}$
are orthonormal matrices, it's easy to  that $\mathcal{U}_{j}^{0}$ and $
\mathcal{V}_{j}^{0}$ are also orthonormal so that $\mathcal{U}_{j}^{0\prime }
\mathcal{U}_{j}^{0}=\mathcal{V}_{j}^{0\prime }\mathcal{V}_{j}^{0}=I_{2}$.
When $j=0$, $\left\{ \lambda _{0,i}^{0}\right\} _{i=1}^{N}$ and $\left\{
f_{0,t}^{0}\right\} _{t=1}^{T}$ are usually referred to as the individual
and time fixed effects, respectively, so that the intercept term exhibits an
additive fixed effects structure.\bigskip

\noindent \textbf{Example 3.} Let $\Theta _{j,it}^{0}=\sum_{k\in \lbrack
K_{j,t}]}\alpha _{j,kt}\mathbf{1}\{i\in G_{j,kt}\}$, where $\left\{
G_{j,kt}\right\} $ forms a partition of $[N]$ for each specific time $t$ and
$K_{j,t}$ is the number of groups at time $t$. Moreover, let
\begin{equation*}
\alpha _{j,kt}=\left\{ \begin{aligned} &\alpha_{j,k}^{(1)},
\quad\text{for}\quad t=1,\dots,T_{b},\\ &
\alpha_{j,k}^{(2)},\quad\text{for}\quad t=T_{b}+1,\dots,T,\\ \end{aligned}
\right.
\end{equation*}
\begin{equation*}
G_{j,kt}=\left\{ \begin{aligned} &G_{j,k}^{(1)},\quad \text{for}\quad
t=1,\dots,T_{b}, k=1,\dots,K_{j}^{(1)},\\ & G_{j,k}^{(2)},\quad
\text{for}\quad t=T_{b}+1,...,T, k=1,...,K_{j}^{(2)},\\ \end{aligned}\right.
\end{equation*}
where $K_{j}^{(1)}$ and $K_{j}^{(2)}$ are the number of groups before and
after the break point $T_{b}$. If $K_{j}^{(1)}=K_{j}^{(2)}$, it is clear
that $rank(\Theta _{j}^{0})=1$. If the group structure does not change after
the break but $\alpha _{j,k}^{(1)}=c\alpha _{j,k}^{(2)}$ for some constant $
c $, we also have $rank(\Theta _{j}^{0})=1$. Except for these two cases, we
can show that
\begin{equation*}
\Theta _{j}^{0}=
\begin{bmatrix}
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(1)}]}\alpha
_{j,k}^{(1)}\mathbf{1}\left\{ 1\in G_{j,k}^{(1)}\right\} , &
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(2)}]}\alpha
_{j,k}^{(2)}\mathbf{1}\left\{ 1\in G_{j,k}^{(2)}\right\} \\
\vdots & \vdots \\
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(1)}]}\alpha
_{j,k}^{(1)}\mathbf{1}\left\{ i\in G_{j,k}^{(1)}\right\} , &
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(2)}]}\alpha
_{j,k}^{(2)}\mathbf{1}\left\{ i\in G_{j,k}^{(2)}\right\} \\
\vdots & \vdots \\
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(1)}]}\alpha
_{j,k}^{(1)}\mathbf{1}\left\{ N\in G_{j,k}^{(1)}\right\} , &
\operatornamewithlimits{\sum}\limits_{k\in \lbrack K^{(2)}]}\alpha
_{j,k}^{(2)}\mathbf{1}\left\{ N\in G_{j,k}^{(2)}\right\}
\end{bmatrix}
\begin{bmatrix}
\iota _{T_{b}} &  & \mathbf{0}_{T_{b}} \\
\mathbf{0}_{T-T_{b}} &  & \iota _{T-T_{b}}
\end{bmatrix}
^{\prime }
\end{equation*}
where $\iota _{T_{b}}$ is a $T_{b}\times 1$ vector of ones and $\mathbf{0}
_{T_{b}}$ is a $T_{b}\times 1$ vector of zeros. In this case, we notice that
$rank(\Theta _{j}^{0})=2$.\bigskip

\noindent \textbf{Example 4.} When $\Theta _{j,it}^{0}=\lambda
_{j,i}^{0\prime }f_{j,t}^{0}$ with $\lambda _{j,i}^{0}$ and $f_{j,t}^{0}$
being two $K_{j}$-vectors, we have the IFEs structure. This is the most
general example without further restrictions.

Like \cite{chernozhukov2019inference}, we assume that for each $j\in \left[ p
\right] ,$ $X_{j,it}$ exhibits a factor structure: $X_{j,it}=\mu
_{j,it}+e_{j,it}=l_{j,i}^{0\prime }w_{j,t}^{0}+e_{j,it},$ where $w_{j,t}^{0}$
and $l_{j,i}^{0}$ are the factors and factor loadings of dimension $r_{j}.$

\subsection{Estimation Algorithm}

In this subsection we provide the estimation algorithm by assuming that $
K_{j}$ are all known for all $j$. In the next subsection, we will introduce
a rank estimation method to estimate $K_{j}$ consistently.

Define the check function $\rho _{\tau }(u)=u\left( \tau -\mathbf{1}\{u\leq
0\}\right) $. The estimation procedure goes as follows:

\begin{itemize}[leftmargin=40pt]
\item[Step 1:] \textbf{Sample Splitting and Nuclear Norm Regularization.} Along the cross-section span, randomly split the sample into three
subsets denoted as $I_{1}$, $I_{2}$ and $I_{3}$, where $I_{\ell }$ has $
N_{\ell }$ individuals such that $N_{1}\approx N_{2}\approx N_{3}\approx N/3$
. Using the data with $(i,t)\in I_{1}\times \lbrack T]$, we run the nuclear
norm regularized quantile regression (QR) and obtain $\{\tilde{\Theta}
_{j}^{(1)}\}_{j\in \{0,\cdots ,p\}}$, i.e.,
\begin{equation}
\{\tilde{\Theta}_{j}^{(1)}\}_{j=0}^{p}=\operatorname*{arg\,min}\limits_{\left\{ \Theta
_{j}\right\} _{j=0}^{p}}\frac{1}{N_{1}T}\sum_{i\in I_{1}}\sum_{t=1}^{T}\rho
_{\tau }(Y_{it}-\sum_{j=1}^{p}X_{j,it}\Theta _{j,it}-\Theta
_{0,it})+\sum_{j=0}^{p}\nu _{j}\left\Vert \Theta _{j}\right\Vert _{\ast },
\label{obj1}
\end{equation}
where $\nu _{j}$ is a tuning parameter. For each $j$, conduct the SVD: $
\frac{1}{\sqrt{N_{1}T}}\tilde{\Theta}_{j}^{(1)}=\hat{\tilde{\mathcal{U}}}
_{j}^{(1)}\hat{\tilde{\Sigma}}_{j}^{(1)}\hat{\tilde{\mathcal{V}}}
_{j}^{(1)\prime }$, where $\hat{\tilde{\Sigma}}_{j}^{(1)}$ is the diagonal
matrix with the diagonal elements being the descending singular values of $
\tilde{\Theta}_{j}^{(1)}$. Let $\tilde{\mathcal{V}}_{j}^{(1)}$ consists the
first $K_{j}$ columns of $\hat{\tilde{\mathcal{V}}}_{j}^{(1)}.$ Let $\tilde{V
}_{j}^{(1)}=\sqrt{T}\tilde{\mathcal{V}}_{j}^{(1)}$ and $\tilde{v}
_{t,j}^{(1)\prime }$ be the $t$-th row of $\tilde{V}_{j}^{(1)}$ $\forall
t\in \lbrack T]$.

\item[Step 2:] \textbf{Row- and Column-Wise Quantile Regression.} Using the
data with $(i,t)\in I_{2}\times \lbrack T]$, we first run the row-wise QR
of\ $Y_{it}$ on $\left( \tilde{v}_{t,0}^{(1)},\left\{ \tilde{v}
_{t,j}^{(1)}X_{j,it}\right\} _{j\in \lbrack p]}\right) $ to obtain $\{\dot{u}
_{i,j}^{(1)}\}_{j=0}^{p}$ for $i\in I_{2}$, and then run the column-wise QR
of $Y_{it}$ on $(\dot{u}_{i,0}^{(1)},\{\dot{u}_{i,j}^{(1)}X_{j,it}\}_{j\in
\lbrack p]})$ to obtain $\{\dot{v}_{t,j}^{(1)}\}_{j=0}^{p}$ for $t\in
\lbrack T]$. That is,
\begin{align}
\{\dot{u}_{i,j}^{(1)}\}_{j=0}^{p}& =\operatorname*{arg\,min}\limits_{\{u_{i,j}\}_{j\in
\lbrack p]\cup \{0\}}}\frac{1}{T}\sum_{t\in \lbrack T]}\rho _{\tau }\left(
Y_{it}-u_{i,0}^{\prime }\tilde{v}_{t,0}^{(1)}-\sum_{j=1}^{p}u_{i,j}^{\prime }
\tilde{v}_{t,j}^{(1)}X_{j,it}\right) ,\forall i\in I_{2},  \label{obj2} \\
\{\dot{v}_{t,j}^{(1)}\}_{j=0}^{p}& =\operatorname*{arg\,min}\limits_{\{v_{t,j}\}_{j\in
\lbrack p]\cup \{0\}}}\frac{1}{N_{2}}\sum_{i\in I_{2}}\rho _{\tau }\left(
Y_{it}-v_{t,0}^{\prime }\dot{u}_{i,0}^{(1)}-\sum_{j=1}^{p}v_{t,j}^{\prime }
\dot{u}_{i,j}^{(1)}X_{j,it}\right) ,\forall t\in \lbrack T].
\end{align}
Similarly, we run the row-wise QR of $Y_{it}$ on $(\dot{v}_{t,0}^{(1)},\{
\dot{v}_{t,j}^{(1)}X_{j,it}\}_{j\in \lbrack p]})$ to obtain $\{\dot{u}
_{i,j}^{(1)}\}_{j=0}^{p}$ for $i\in I_{3}$, i.e.,
\begin{equation*}
\{\dot{u}_{i,j}^{(1)}\}_{j=0}^{p}=\operatorname*{arg\,min}\limits_{\{u_{i,j}\}_{j\in \lbrack
p]\cup \{0\}}}\frac{1}{T}\sum_{t\in \lbrack T]}\rho _{\tau }\left(
Y_{it}-u_{i,0}^{\prime }\dot{v}_{t,0}^{(1)}-\sum_{j=1}^{p}u_{i,j}^{\prime }
\dot{v}_{t,j}^{(1)}X_{j,it}\right) ,\forall i\in I_{3}.
\end{equation*}

\item[Step 3:] \textbf{Debiasing.}

\begin{itemize}
\item[Step 3.1:] For each $j\in \lbrack p]$, we conduct the principle
component analysis (PCA) for $X_{j,it}$ with $\left( i,t\right) \in \lbrack
N]\times \lbrack T]$ to obtain the factor and factor loading estimates as
\begin{equation}
\left\{ \hat{l}_{j,i},\hat{w}_{j,t}\right\} _{i\in \lbrack N],t\in \lbrack
T]}=\operatornamewithlimits{\operatorname*{arg\,min}}\limits_{\left\{ l_{j,i},w_{j,t}\right\}
_{i\in \lbrack N],t\in \lbrack T]}}\frac{1}{NT}\sum_{i\in \lbrack
N]}\sum_{t\in \lbrack T]}\left( X_{j,it}-l_{j,i}^{\prime }w_{j,t}\right)
^{2},  \label{debias_0}
\end{equation}
subject to the normalizations: $\frac{1}{N}\sum_{i=1}^{N}l_{i,j}l_{i,j}^{
\prime }=I_{r_{j}}$ and $\frac{1}{T}\sum_{t=1}^{T}w_{j,t}w_{j,t}^{\prime }$
is a diagonal matrix with descending diagonal elements. Then we define $\hat{
\mu}_{j,it}=\hat{l}_{j,i}^{\prime }\hat{w}_{j,t}$ and $\hat{e}
_{j,it}=X_{j,it}-\hat{\mu}_{j,it}$.

\item[Step 3.2:] For $(i,t)\in I_{3}\times [T]$, let $\tilde{Y}_{it}=Y_{it}-
\operatornamewithlimits{\sum}\limits_{j=1}^{p}\hat{\mu}_{j,it}\dot{u}
_{i,j}^{(1)\prime}\dot{v}_{t,j}^{(1)}.$ We run the row-wise QR\ of $\tilde{Y}
_{it}$ on $(\dot{v}_{t,0}^{(1)},\{\dot{v}_{t,j}^{(1)}\hat{e}_{j,it}\}_{j\in[p
]})$ to obtain the final estimates $\hat{u}_{i,j}^{(3,1)}$, i.e.,
\begin{equation}
\{\hat{u}_{i,j}^{(3,1)}\}_{j=0}^{p}=\operatorname*{arg\,min}\limits_{\{u_{i,j}\}_{j=0}^{p}}
\frac{1}{T}\sum_{t\in[T]}\rho_{\tau }\left(\tilde{Y}_{it}-u_{i,0}^{\prime}
\dot{v}_{t,0}^{(1)}-\sum_{j=1}^{p}u_{i,j}^{\prime}\dot{v}_{t,j}^{(1)}\hat{e}
_{j,it}\right) ,\forall i\in I_{3}.  \label{debias_1}
\end{equation}
Updating $\hat{Y}_{it}=Y_{it}-\sum_{j=1}^{p}\hat{\mu}_{j,it}\hat{u}
_{i,j}^{(3,1)\prime}\dot{v}_{t,j}^{(1)}$, we run the column-wise QR\ of $
\hat{Y}_{it}$ on $(\hat{u}_{i,0}^{(3,1)},\{\hat{u}_{i,j}^{(3,1)}\hat{e}
_{j,it}\}_{j\in[p]})$ to obtain $\hat{v}_{t,j}^{(3,1)}$, i.e.,
\begin{equation}
\{\hat{v}_{t,j}^{(3,1)}\}_{j=0}^{p}=\operatorname*{arg\,min}\limits_{\{v_{t,j}\}_{j=0}^{p}}
\frac{1}{N_{3}}\sum_{i\in I_{3}}\rho_{\tau }\left(\hat{Y}_{it}-v_{t,0}^{
\prime}\hat{u}_{i,0}^{(3,1)}-\sum_{j=1}^{p}v_{t,j}^{\prime}\hat{u}
_{i,j}^{(3,1)}\hat{e}_{j,it}\right) ,\forall t\in[T].  \label{debias_2}
\end{equation}
\end{itemize}
\end{itemize}

In order to obtain the final estimators for the full sample, we propose to
switch the role of each subsample for the low-rank estimation, row- and
column-wise QR and debiasing, then repeat Steps 1-3 to obtain $\left\{ \hat{u
}_{i,j}^{(a,b)}\right\} _{j=0}^{p}$ and $\left\{ \hat{v}_{t,j}^{(a,b)}\right
\} _{j=0}^{p}$ for $a\in \lbrack 3]$ and $b\in \lbrack 3]\setminus \{a\}$.
Here $(a,b)$ denotes the final estimates for subsample $I_{a}$ obtained from
the first step NNR estimates with subsample $I_{b}$. Table \ref
{tab:estimator symbol} shows the final estimators we obtain by using
different combination of subsamples.
\begin{table}[th]
\caption{Estimators using different subsamples at different steps in
algorithm.}
\label{tab:estimator symbol}\centering
\begin{tabular}{cccc}
\toprule \toprule Step 1 $(b)$ & Step 2 & Step 3 $(a)$ & estimators $(a,b)$ \\
\midrule $I_{1}$ & $I_{2}$ & $I_{3}$ & $\hat{u}_{i,j}^{(3,1)}$, $\hat{v}
_{t,j}^{(3,1)}$ \\
$I_{2}$ & $I_{1}$ & $I_{3}$ & $\hat{u}_{i,j}^{(3,2)}$, $\hat{v}
_{t,j}^{(3,2)} $ \\
$I_{1}$ & $I_{3}$ & $I_{2}$ & $\hat{u}_{i,j}^{(2,1)}$, $\hat{v}
_{t,j}^{(2,1)} $ \\
$I_{3}$ & $I_{1}$ & $I_{2}$ & $\hat{u}_{i,j}^{(2,3)}$, $\hat{v}
_{t,j}^{(2,3)} $ \\
$I_{2}$ & $I_{3}$ & $I_{1}$ & $\hat{u}_{i,j}^{(1,2)}$, $\hat{v}
_{t,j}^{(1,2)} $ \\
$I_{3}$ & $I_{2}$ & $I_{1}$ & $\hat{u}_{i,j}^{(1,3)}$, $\hat{v}
_{t,j}^{(1,3)} $ \\
\bottomrule &  &  &
\end{tabular}
\end{table}

Several remarks are in order. First, we randomly split the full sample into
three subsamples, each playing a significant role in the algorithm. We use
the first subsample for the low-rank estimation to obtain the preliminary
NNR estimators of the submatrices of the intercept and slope matrices. But
these estimators are only consistent in terms of Frobenius norm, and one
cannot derive the pointwise or uniform convergence rates for them. With the
low-rank estimates, we use the second subsample to do the row- and
column-wise QRs and can now establish the uniform convergence rates for each
row of factor and factor loading estimators. Then we use the remaining
subsample to debias the second-stage estimator and to obtain the final
estimators that have the desirable asymptotic properties.

Second, to reduce the randomness of sample splitting, one can run the
estimation algorithm several times with different splittings in practice.
Once one obtains factor and factor loading estimates, one can construct
estimators for $\Theta _{j}^{0}$ under different splittings and then choose
the one specific splitting which yields the minimum quantile objective function.

Third, the bias in the second-stage estimator is inherent from the
first-stage NNR estimator. We follow the lead of \cite
{chernozhukov2019inference} to assume that $X_{j,it}$ has a factor structure
with an additive idiosyncratic term, and remove the bias by a QR with the
demeaned $X_{j,it}$ as regressors. In the least squares panel regression
framework, the objective function is smooth and one has closed-form
solutions in the last stage so that \cite{chernozhukov2019inference} only
need to split the sample into two subsamples. In contrast, in the PQR
framework, the objective function is non-smooth, we do not have closed-form
solutions in any stage. In order to remove the bias from the early stage
estimation and to derive the distributional results, we need to split the
sample into three subsamples.





To save space, we relegate the detailed algorithm for the nuclear norm
regularization to the online supplement.

\subsection{Rank Estimation}

In this subsection we discuss how to estimate the ranks $K_{j}$
consistently. To estimate the ranks, we consider the full sample NNR QR
estimation:
\begin{equation}
\{\tilde{\Theta}_{j}\}_{j=0}^{p}=\operatorname*{arg\,min}\limits_{\left\{ \Theta _{j}\right\}
_{j=0}^{p}}\,\,\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\rho _{\tau
}(Y_{it}-\sum_{j=1}^{p}X_{j,it}\Theta _{j,it}-\Theta
_{0,it})+\sum_{j=0}^{p}\nu _{j}\left\Vert \Theta _{j}\right\Vert _{\ast }.
\label{pre rank}
\end{equation}
For $j\in \left\{ 0,\cdots ,p\right\} $, we estimate $K_{j}$ by the popular
singular value thresholding (SVT) as follows
\begin{equation*}
\hat{K}_{j}=\sum_{m}\mathbf{1}\left\{ \lambda _{m}\left( \tilde{\Theta}
_{j}\right) \geq 0.6\left( NT\nu _{j}\left\Vert \tilde{\Theta}
_{j}\right\Vert _{op}\right) ^{1/2}\right\} .
\end{equation*}
It is standard to show that $\mathbb{P}(\hat{K}_{j}=K_{j})\rightarrow 1$ as $
\left( N,T\right) \rightarrow \infty $ under some regularity conditions
given in the next section; see also Proposition D.1 in \cite
{chernozhukov2019inference} and Theorem 2 in \cite{Hong_Su_Jiang2022}. Since
the ranks can be estimated consistently, we assume that they are known in
the asymptotic theory below.


\section{Asymptotic Theory}

In this section, we study the asymptotic properties of the estimators
introduced in the last section.

\subsection{First Stage Estimator}

Recall that $X_{j,it}=\mu _{j,it}+e_{j,it}=l_{j,i}^{0\prime
}w_{j,t}^{0}+e_{j,it}$ for each $j\in \lbrack p]$. Let $
X_{it}=(X_{1,it},...,X_{p,it})^{\prime }$ and $
e_{it}=(e_{1,it},...,e_{p,it})^{\prime }.$ Define $\epsilon _{i}=\left(
\epsilon _{i1},\cdots ,\epsilon _{it}\right) ^{\prime }$, $e_{j,i}=\left(
e_{j,i1},\cdots ,e_{j,iT}\right) ^{\prime }$, $W_{j}^{0}$ as the $T\times
r_{j}$ matrix with each row being $w_{j,t}^{0\prime }$, and $V_{j}^{0}$ as
the $T\times K_{j}$ matrix with each row being $v_{t,j}^{0\prime }$. Further
define $a_{it}=\tau -\mathbf{1}\left\{ \epsilon _{it}\leq 0\right\} $ with $
a_{i}=\left( a_{i1},\cdots ,a_{iT}\right) ^{\prime }$ and $a=\left(
a_{1},\cdots ,a_{N}\right) ^{\prime }$. Throughout the paper, we treat the
factors $\{v_{t,j}^{0}\}_{t\in \lbrack T],j\in \lbrack p]\cup \{0\}}$ and $
\{w_{j,t}^{0}\}_{t\in \lbrack T],j\in \lbrack p]}$ as random and their
loadings $\{u_{i,j}^{0}\}_{i\in \lbrack N],j\in \lbrack p]\cup \{0\}}$ and $
\{l_{j,i}^{0}\}_{i\in \lbrack N],j\in \lbrack p]}$ as deterministic.

Table \ref{sigma fields} defines several $\sigma$-fields. We use $\mathscr{D}
$ to denote the minimal $\sigma$-field generated by $\left\{V_{j}^{0}\right
\}_{j\in [p]\cup \{0\}}\bigcup \left\{W_{j}^{0}\right\}_{j\in [p]};$ the
superscripts $I_{1}$ and $I_{1}\cup I_{2}$ are associated with the first
subsample and the first two subsamples, respectively.  For example, $
\mathscr{D}_{e_{i}}^{I_{1}}$ denotes the minimal $\sigma$-field generated by
$\mathscr{D},$ $\left\{e_{it}\right\}_{t\in \left[T\right] }$ and $\left\{
\epsilon_{it},e_{it}\right\}_{i\in I_{1},t\in \left[T\right] }.$
\begin{table}[h]
\caption{Definition of various $\protect\sigma $-fields}
\label{sigma fields}\centering
\begin{tabular}{cc}
\toprule\toprule Notation & $\sigma$-fields generated by \\
\midrule$\mathscr{D}$ & $\left\{ V_{j}^{0}\right\}_{j\in [p]\cup
\{0\}}\bigcup \left\{ W_{j}^{0}\right\}_{j\in [p]}$ \\
\midrule$\mathscr{D}_{e_{it}}$ & $\mathscr{D}\bigcup e_{it}$ \\
\midrule$\mathscr{D}_{e_{i}}$ & $\mathscr{D}\bigcup \left\{ e_{it}\right\}
_{t\in[T]}$ \\
\midrule$\mathscr{D}_{e}$ & $\mathscr{D}\bigcup \left\{ e_{it}\right\} _{i\in
[N], t\in[T]}$ \\
\midrule$\mathscr{D}^{I_{1}\cup I_{2}}$ & $\mathscr{D}\bigcup \left\{
\epsilon_{it},e_{it}\right\}_{i\in I_{1}\cup I_{2},t\in[T]}$ \\
\midrule$\mathscr{D}^{I_{1}}_{\{e_{is}\}_{s<t}}$ & $\mathscr{D}\bigcup
\{e_{is}\}_{s<t}\bigcup \left\{
\epsilon_{i^{*}t^{*}},e_{i^{*}t^{*}}\right\}_{i^{*}\in I_{1},t^{*}\in[T]}$
\\
\midrule$\mathscr{D}_{e_{i}}^{I_{1}}$ & $\mathscr{D}\bigcup \left\{
e_{it}\right\}_{t\in[T]}\bigcup \left\{ \epsilon
_{i^{*}t^{*}},e_{i^{*}t^{*}}\right\}_{i^{*}\in I_{1},t^{*}\in[T]}$ \\
\midrule$\mathscr{D}_{e_{i}}^{I_{1}\cup I_{2}}$ & $\mathscr{D}\bigcup
\left\{ e_{it}\right\}_{t\in[T]}\bigcup \left\{ \epsilon
_{i^{*}t^{*}},e_{i^{*}t^{*}}\right\}_{i^{*}\in I_{1}\cup I_{2},t^{*}\in[T]}$
\\
\midrule$\mathscr{D}_{e}^{I_{1}\cup I_{2}}$ & $\mathscr{D}\bigcup \left\{
e_{it}\right\}_{i\in [N,]t\in[T]}\bigcup \left\{ \epsilon
_{it},e_{it}\right\}_{i\in I_{1}\cup I_{2},t\in[T]}$ \\
\bottomrule &
\end{tabular}
\end{table}

Let $M$ denote a generic bounded constant that may vary across places. Let $
\mathscr{G}_{i,t-1}$ denote the minimal $\sigma $-field generated by $
\mathscr{D}\cup \{e_{ls}\}_{l\leq i-1,s\in \lbrack T]}\cup \{e_{is}\}_{s\leq
t}\cup \{\epsilon _{ls}\}_{l\leq i-1,s\in \lbrack T]}\cup \{\epsilon
_{is}\}_{s\leq t-1}$. Let $\mathsf{F}_{it}(\cdot )$ and $\mathsf{f}
_{it}(\cdot )$ be the conditional cumulative distribution function (CDF) and
probability density function (PDF) of $\epsilon _{it}$ given $\mathscr{G}
_{i,t-1}$, respectively. Similarly, let $\mathfrak{F}_{it}(\cdot )$ and $
\mathfrak{f}_{it}(\cdot )$ denote the conditional CDF and PDF of $\epsilon
_{it}$ given $\mathscr{D}_{e_{i}};$ $F_{it}(\cdot )$ and $f_{it}(\cdot )$
denote the conditional CDF and PDF of $\epsilon _{it}$ given $\mathscr{D}
_{e}.$ Let $\mathsf{f}_{it}^{\prime }\left( \cdot \right) $, $\mathfrak{f}
_{it}^{\prime }\left( \cdot \right) ,$ and $f_{it}^{\prime }\left( \cdot
\right) $ denotes the first derivative of the density $\mathsf{f}_{it}\left(
\cdot \right) $, $\mathfrak{f}_{it}\left( \cdot \right) ,$ and $f_{it}\left(
\cdot \right) ,$ respectively.

We make the following assumptions.

\begin{ass}
\begin{itemize}
\item[(i)] $\left\{ \epsilon_{it},e_{it}\right\}_{t\in[T]}$ are
conditionally independent across $i$ given $\mathscr{D}$.

\item[(ii)] $\mathbb{E}\left(a_{it}\bigg |\mathscr{D}_{e}\right) =0$.

\item[(iii)] For each $i$, $\left\{ \epsilon_{it},t\geq 1\right\} $ is
strong mixing conditional on $\mathscr{D}_{e_{i}}$, and $\{\left(
\epsilon_{it},e_{it}\right) ,t\geq 1\}$ is strong mixing conditional on $
\mathscr{D}$. Both mixing coefficients are upper bounded by $
\alpha_{i}(\cdot)$ such that $\max_{i\in [N]}\alpha_{i}(z)\leq M\alpha ^{z}$
for some constant $\alpha \in \left(0,1\right) $.

\item[(iv)] $\max_{i\in [N]}\frac{1}{T}\sum_{t\in [ T]}\left\Vert
X_{it}\right\Vert_{2}^{3}\leq M~$a.s., $\max_{t\in[T]}\frac{1}{N_{2}}
\sum_{i\in I_{2}}\left\Vert X_{it}\right\Vert_{2}^{4}\leq M~$a.s., \newline
$\max_{\left(i,t\right) \in [N]\times [T]}\mathbb{E}\left[ \left\Vert
X_{it}\right\Vert_{2}^{3}\bigg|\mathscr{D}\right] \leq M$a.s., $
\max_{i\in [N]}\sqrt{\frac{1}{T}\sum_{t\in[T]}\left[ \mathbb{E}
\left(\epsilon_{it}^{2}\bigg|\mathscr{D}_{e_{i}}\right) \right] ^{2}}$ $\leq
M~a.s.$, and \newline
$\max_{\left(i,t\right) \in [N]\times [T]}\mathbb{E}\left(\left\Vert
X_{it}\right\Vert_{2}^{2}\bigg|\mathscr{D}_{\left\{e_{is}\right\}_{s<t}}
\right)$ $\leq M~a.s.$

\item[(v)] For $j\in [p]$, there exists a positive sequence $\xi_{N}$ such
that $\max_{\left(i,t\right) \in [N]\times [T]}\left\vert
X_{j,it}\right\vert \leq \xi_{N}~a.s.$

\item[(vi)] $\min_{\left(i,t\right) \in [N]\times [T]} \mathsf{f}_{it}(0)\geq \underline{\mathsf{f}}>0$ and $
\max_{\left(i,t\right) \in [N]\times [T]}\sup_{\epsilon}\left\vert \mathsf{f}
_{it}^{\prime }\left(\epsilon \right) \right\vert\leq \bar{\mathsf{f}}
^{\prime }$.

\item[(vii)] $\min_{\left(i,t\right) \in [N]\times [T]}\mathfrak{f}_{it}(0)\geq \underline{\mathfrak{f}}>0$ and $
\max_{\left(i,t\right) \in [N]\times [T]}\sup_{\epsilon}\left\vert \mathfrak{
f}_{it}^{\prime }\left(\epsilon \right) \right\vert\leq \bar{\mathfrak{f}}
^{\prime }$.

\item[(viii)] $\min_{\left(i,t\right) \in [N]\times [T]}f_{it}(0)\geq \underline{f}>0$ and $\max_{\left(i,t\right)\in
[N]\times [T]}\sup_{\epsilon }\left\vert f_{it}^{\prime}\left(\epsilon
\right) \right\vert \leq \bar{f}^{\prime }$.

\item[(ix)] $\frac{\xi_{N}^{4}\log (N\vee T)\sqrt{N\vee T}}{N\wedge T}=o(1)$
and $\frac{\left(\frac{N}{T}\vee 1\right) ^{1/2}}{\left(N\wedge T\right) ^{
\frac{1}{4+2\vartheta }}}\left(\log (N\vee T)\right) ^{\frac{3+\vartheta }{
4+2\vartheta }}\xi_{N}^{\frac{5+\vartheta }{2+\vartheta }}=o(1)$ for any $
\vartheta >0$.
\end{itemize}

\label{ass:1}
\end{ass}

Assumptions \ref{ass:1}(i) imposes conditional independence of the error
terms and covariates $X_{j,it}$ given the fixed effects. Assumptions \ref
{ass:1}(ii) imposes the moment condition for QR. Assumptions \ref{ass:1}
(iii) imposes the weak dependence assumption along the time dimension via
the use of the notion of conditional strong mixing. See \cite
{Prakasa_Rao2009} for the definition of conditional strong mixing and \cite
{Su_Chen2013} for an application in the panel setup. Assumptions \ref{ass:1}
(iv)-(v) essentially imposes some conditions on the moments and tail
behavior of the both covariates and errors. Note that we allow $X_{j,it}$ to
have an infinite support. Assumptions \ref{ass:1}(vi)-(viii), which are used
in the proofs of Theorems \ref{Thm1}, \ref{Thm2} and \ref{Thm3},
respectively, specify conditions on the conditional density of $\epsilon
_{it}$ given different $\sigma $-fields. Assumption \ref{ass:1}(ix) imposes
some restrictions on $N$, $T$ and $\xi _{N}$ in order to obtain the error
bound of NNR estimators and to achieve the unbiasedness. It allows not only the case that $N$ and $T$ diverge to infinity at the the same rate, but also the case that $N$ diverges to infinity not too faster than $T$, and vice versa.

\begin{ass}
\label{ass:2} $\Theta _{0}^{0}$ is the fixed effect matrix with fixed rank $
K_{0}$ and $\left\Vert \Theta _{0}^{0}\right\Vert _{\max }\leq M$. For each $
j\in \lbrack p],$ $\Theta _{j}^{0}$ is the slope matrix of regressor $j$
with rank being $K_{j}$ such that $\max_{j\in \lbrack p]}\left\Vert \Theta
_{j}\right\Vert _{\max }\leq M$ and $\max_{j\in \lbrack p]}K_{j}\leq \bar{K}$
for some fixed finite $\bar{K}.$
\end{ass}

Assumption \ref{ass:2} is the low-rank assumption for the intercept and
slope matrices, which is the key assumption for the NNR.
The uniform boundedness of elements of these matrices facilitates the
asymptotic analysis,  but can be relaxed at the cost of more lengthy argument.
See \cite{ma2020detecting} for a similar condition.

\begin{ass}
\label{ass:3} There exist some constants $C_{\sigma}$ and $c_{\sigma}$ such
that
\begin{align*}
\infty>C_{\sigma}\geq \lim\sup_{N,T} \max_{j\in [p]\cup\{0\}}
\sigma_{1,j}\geq \lim\inf_{N,T} \min_{j\in [p]\cup\{0\}}
\sigma_{K_{j},j}\geq c_{\sigma}>0.
\end{align*}
\end{ass}

Assumption \ref{ass:3} imposes some conditions on the singular values of the
coefficient matrices. It implies that we only allow pervasive factors when
these matrices are written as a factor structure. Such an assumption is
common in the literature; see, e.g., Assumption 3 in \cite{ma2020detecting}.

To introduce the next assumption, we need some notation. Let $\Theta
_{j}^{0}=R_{j}\Sigma _{j}S_{j}^{\prime }$ be the SVD for $\Theta _{j}^{0}$.
Further decompose $R_{j}=\left( R_{j,r},R_{j,0}\right) $ with $R_{j,r}$
being the singular vectors corresponding to the nonzero singular values, $
R_{j,0}$ being the singular vectors corresponding to the zero singular
values. Decompose $S_{j}=\left( S_{j,r},S_{j,0}\right) $ with $S_{j,r}$ and $
S_{j,0}$ defined analogously. For any matrix $W\in \mathbb{R}^{N\times T}$,
we define
\begin{equation*}
\mathcal{P}_{j}^{\bot }\left( W\right) =R_{j,0}R_{j,0}^{\prime
}WS_{j,0}S_{j,0}^{\prime },\quad \mathcal{P}_{j}\left( W\right) =W-\mathcal{P
}_{j}^{\bot }\left( W\right) ,
\end{equation*}
where $\mathcal{P}_{j}\left( W\right) $ and $\mathcal{P}_{j}^{\bot }\left(
W\right) $ are the linear projection of matrix $W$ onto the low-rank space
and its orthogonal space, respectively. Let $\Delta _{\Theta _{j}}=\Theta
_{j}-\Theta _{j}^{0}$ for any $\Theta _{j}$. With some positive constants $C_{1}$ and $C_{2}$, we define the
following cone-like restricted set:
\begin{equation*}
\mathcal{R}(C_{1},C_{2}):=\left\{ \left(\{\Delta _{\Theta
_{j}}\}_{j=0}^{p}\right):\sum_{j=0}^{p}\left\Vert \mathcal{P}_{j}^{\bot }(\Delta
_{\Theta _{j}})\right\Vert _{\ast }\leq C_{1}\sum_{j=0}^{p}\left\Vert
\mathcal{P}_{j}(\Delta _{\Theta _{j}})\right\Vert _{\ast
},\,\sum_{j=0}^{p}\left\Vert \Delta _{\Theta _{j}}\right\Vert _{F}^{2}\geq
C_{2}\sqrt{NT}\right\} .
\end{equation*}


\begin{ass}
\label{ass:4} Let $C_{2}>0$ be a sufficiently large but fixed constant.
There are constants $C_{3},C_{4}$, such that, uniformly over $(\{\Delta
_{\Theta _{j}}\}_{j=0}^{p})\in \mathcal{R}(3,C_{2})$, we have
\begin{equation*}
\left\Vert \Delta _{\Theta _{0}}+\sum_{j=1}^{p}\Delta _{\Theta _{j}}\odot
X_{j}\right\Vert _{F}^{2}\geq C_{3}\sum_{j=0}^{p}\left\Vert \Delta _{\Theta
_{j}}\right\Vert _{F}^{2}-C_{4}(N+T)~w.p.a.1.
\end{equation*}
The same condition holds when $\Theta_j^0$ is replaced by $\{\Theta_{j,it}^0\}_{i \in I_a, t \in [T]}$ for $a = 1,2,3.$
\end{ass}

Assumption \ref{ass:4} parallels the restricted strong convexity (RSC)
condition in Assumption 3.1 of \cite{chernozhukov2019inference} who also
provide some sufficient primitive conditions.

For any $j\in \left\{ 0,\cdots ,p\right\} $, define $\tilde{\Delta}_{\Theta
_{j}}=\tilde{\Theta}_{j}-\Theta _{j}^{0}$ and $\tilde{\Delta}_{\Theta
_{j}}^{(1)}=\tilde{\Theta}_{j}^{(1)}-\Theta _{j}^{0,(1)}$, where $\Theta
_{j}^{0,(1)}=\left\{ \Theta _{j,it}^{0}\right\} _{i\in I_{1},t\in \lbrack
T]} $. The following theorem establishes the convergence rates of the NNR
estimators of the coefficient matrices.

\begin{theorem}
\label{Thm1} If Assumptions \ref{ass:1}-\ref{ass:4} hold, for $\forall
j\in\left\{0,\cdots,p\right\}$, we have

\begin{itemize}
\item[(i)] $\frac{1}{\sqrt{NT}}\left\Vert \tilde{\Delta}_{\Theta_{j}}\right
\Vert_{F}=O_{p}\left(\sqrt{\frac{\log (N\vee T)}{N\wedge T}}
\xi_{N}^{2}\right) $, $\frac{1}{\sqrt{NT}}\left\Vert \tilde{\Delta}
_{\Theta_{j}}^{(1)}\right\Vert_{F}=O_{p}\left(\sqrt{\frac{\log (N\vee T)}{
N\wedge T}}\xi_{N}^{2}\right) $,

\item[(ii)] $\max_{k\in [K_{j}]}\left\vert \tilde{\sigma}_{k,j}-\sigma_{k,j}
\right\vert =O_{p}\left(\sqrt{\frac{\log (N\vee T)}{N\wedge T}}
\xi_{N}^{2}\right) $, $\max_{k\in [K_{j}]}\left\vert \tilde{\sigma}
_{k,j}^{(1)}-\sigma_{k,j}\right\vert =O_{p}\left(\sqrt{\frac{\log (N\vee T)}{
N\wedge T}}\xi_{N}^{2}\right) $,

\item[(iii)] $\frac{1}{\sqrt{T}}\left\Vert V_{j}^{0}-\tilde{V}
_{j}O_{j}\right\Vert_{F}=O_{p}\left(\sqrt{\frac{\log (N\vee T)}{N\wedge T}}
\xi_{N}^{2}\right) $, $\frac{1}{\sqrt{T}}\left\Vert V_{j}^{0}-\tilde{V}
_{j}^{(1)}O_{j}^{(1)}\right\Vert_{F}=O_{p}\left(\sqrt{\frac{\log \,N\vee T}{
N\wedge T}}\xi_{N}^{2}\right) $,
\end{itemize}

where $O_{j}$ and $O_{j}^{(1)}$ are some orthogonal rotation matrices
defined in the proof.
\end{theorem}

\textbf{Remark 1.} Theorem \ref{Thm1}(i) reports the \textquotedblleft
rough\textquotedblright\ convergence rates of the NNR estimators of the
coefficient matrices in terms of Frobenius norm for both the full-sample and
sub-sample estimators. Unlike the traditional $\left( N\wedge T\right)
^{-1/2}$-rate in the least squares framework, NNR estimators' convergence rates in the PQR
framework usually have an additional $\sqrt{\log (N\vee T)}$ term due to the use
of some exponential inequalities. The extra term $\xi
_{N}^{2}$\ in our rate is due to the upper bound of $|X_{j,it}|$, and it
disappears in case $X_{j,it}$'s are uniformly bounded. Theorem \ref{Thm1}(ii)-(iii) report the convergence rates for the estimators
of the factors and factor loadings of $\Theta _{j}^{0}$, which are inherited from those in Theorem \ref{Thm1}(i).  To derive these results, we establish the symmetrization inequality and contraction principle for the sequential symmetrization developed by \cite{rakhlin2015sequential}. See Lemmas \ref{Lem:symmetrization} and \ref{Lem:contraction} in the online supplement for more detail.






\subsection{Second Stage Estimator}

To study the asymptotic properties of the second-stage estimators, we add some notation. Define
\begin{equation*}
\Phi_{i}=\frac{1}{T}\sum_{t=1}^{T}\Phi_{it}^{0}\Phi_{it}^{0\prime}\quad \text{and }\Psi_{t}=\frac{1}{N_{2}}\sum_{i\in I_{2}}\Psi_{it}^{0}\Psi_{it}^{0\prime},
\end{equation*}
where $\Phi_{it}^{0}=(v_{t,0}^{0\prime},v_{t,1}^{0\prime}X_{1,it},\cdots,v_{t,p}^{0\prime}X_{p,it})^{\prime}$ and $\Psi_{it}^{0}=(u_{i,0}^{0\prime},u_{i,1}^{0\prime}X_{1,it},\cdots,u_{i,p}^{0\prime}X_{p,it})^{\prime}.$ Let $K=\sum_{j=0}^{p}K_{j}.$ Note that $\Phi_{i}$ and $\Psi_{t}$ are $K\times K$ matrices. We add the following two assumptions.

\begin{ass}
\label{ass:5} There exist constants $C_{\phi }$ and $c_{\phi }$ such that
a.s.
\begin{align*}
\infty & >C_{\psi }\geq \limsup_{T}\max_{t\in[T]}\lambda_{\max}(\Psi_{t})\geq \liminf_{T}\min_{t\in[T]}\lambda_{\min }(\Psi_{t})\geq c_{\psi }>0, \\
\infty & >C_{\phi }\geq \limsup_{N}\max\limits_{i\in I_{2}}\lambda_{\max}(\Phi_{i})\geq \liminf_{N}\min_{i\in I_{2}}\lambda_{\min}(\Phi_{i})\geq c_{\phi }>0.
\end{align*}
\end{ass}

Assumption \ref{ass:5} is similar to Assumption 8 in \cite{ma2020detecting}.
To introduce Theorem \ref{Thm2}, we define
\begin{align*}
& \dot{\varpi}_{it}=\left(\dot{v}_{t,0}^{(1)\prime},\dot{v}
_{t,1}^{(1)\prime}X_{1,it},\cdots ,\dot{v}_{t,p}^{(1)\prime}X_{p,it}\right)
^{\prime}, \\
& \varpi_{it}^{0}=\left(\left(O_{0}^{(1)}v_{t,0}^{0}\right)
^{\prime},\left(O_{1}^{(1)}v_{t,1}^{0}\right) ^{\prime}X_{1,it},\cdots
,\left(O_{p}^{(1)}v_{t,p}^{0}\right) ^{\prime}X_{p,it}\right) ^{\prime}, \\
& u_{i}^{0}=\left(u_{i,0}^{0\prime},\cdots ,u_{i,p}^{0\prime}\right)
^{\prime},\quad \dot{\Delta}_{t,j}=O_{j}^{(1)\prime}\dot{v}
_{t,j}^{(1)}-v_{t,j}^{0},\quad \dot{\Delta}_{t,v}=\left(\dot{\Delta}
_{t,0}^{\prime},\cdots ,\dot{\Delta}_{t,p}^{\prime}\right) ^{\prime}, \\
& \dot{\Delta}_{i,j}=O_{j}^{(1)\prime}\dot{u}_{i,j}^{(1)}-u_{i,j}^{0},\quad
\dot{\Delta}_{i,u}=\left(\dot{\Delta}_{i,0}^{\prime},\cdots ,\dot{\Delta}
_{i,p}^{\prime}\right) ^{\prime}, \\
& D_{i}^{I}=\frac{1}{T}\sum_{t=1}^{T}\mathfrak{f}_{it}(0)\varpi_{it}^{0}
\varpi_{it}^{0\prime},\quad D_{i}^{II}=\frac{1}{T}\sum_{t=1}^{T}\left[ \tau -
\mathbf{1}\left\{ \epsilon_{it}\leq 0\right\} \right] \varpi_{it}^{0}, \\
& \mathbb{J}_{i}\left(\left\{ \dot{\Delta}_{t,v}\right\}_{t\in[T]}\right) =
\frac{1}{T}\sum_{t=1}^{T}\left[ \mathbf{1}\left\{ \epsilon_{it}\leq
0\right\} -\mathbf{1}\left\{ \epsilon_{it}\leq \dot{\Delta}
_{t,v}^{\prime}\Psi_{it}^{0}\right\} \right] \varpi_{it}^{0}.
\end{align*}
Theorem \ref{Thm2} below gives the uniform convergence rate and linear
expansion of the factor loading estimators from second stage estimation.

\begin{theorem}
\label{Thm2} Suppose Assumptions \ref{ass:1}-\ref{ass:5} hold. Then for each
$j\in \left\{ 0,\cdots ,p\right\} $, we have

\begin{itemize}
\item[(i)] $\operatornamewithlimits{\max}\limits_{i\in I_{2}\cup I_{3}}
\left\Vert\dot{u}_{i,j}^{(1)}-O_{j}^{(1)}u_{i,j}^{0}\right\Vert_{2}=O_{p}
\left(\sqrt{\frac{\log(N \vee T)}{N\wedge T}}\xi_{N}^{2 }\right)$,

\item[(ii)] $\operatornamewithlimits{\max}\limits_{t\in[T]} \left\Vert\dot{v}
_{t,j}^{(1)}-O_{j}^{(1)}v_{t,j}^{0}\right\Vert_{2}=O_{p}\left(\sqrt{\frac{
\log(N \vee T)}{N\wedge T}}\xi_{N}^{2 }\right)$,

\item[(iii)] $\dot{\Delta}_{i,u}=\left[ D_{i}^{I}\right] ^{-1}\left[
D_{i}^{II}+\mathbb{J}_{i}\left( \left\{ \dot{\Delta}_{t,v}\right\} _{t\in
\lbrack T]}\right) \right] +o_{p}\left( \left( N\vee T\right) ^{-1/2}\right)
$ uniformly over $i\in I_{3}$.
\end{itemize}
\end{theorem}

\textbf{Remark 2.} Theorem \ref{Thm2}(i) reports the uniform convergence
rate for the factor loading estimators of $\Theta _{j}^{0}$ for $i\in
I_{2}\cup I_{3};$ Theorem \ref{Thm2}(ii) reports the uniform convergence
rate for the factor estimators of $\Theta _{j}^{0}$ for $t\in \left[ T\right]
$; Theorem \ref{Thm2}(iii) reports the linear expansion for the factor
loading estimators of $\Theta _{j}^{0}$ for $i\in I_{3}$. However, the $\mathbb{J}_{i}\left( \left\{ \dot{\Delta}_{t,v}\right\} _{t\in
\lbrack T]}\right)$ term is not mean-zero and represents the bias induced by the first stage NNR. In the third stage below, we aim to remove such a bias from the linear expansion.


\subsection{Third Stage Estimator}

In the debiasing stage, we first apply PCA to all independent variables $X_{j,it}$, and then run the row- and column-wise QRs to obtain the final
estimators. Below we give Assumptions \ref{ass:7}-\ref{ass:9} for the PCA
procedure and establish the asymptotic linear expansions of PCA estimates in
the online supplement. Theorem \ref{Thm3} below gives the asymptotic
distribution of our final factor and factor loading estimates.

\begin{ass}
\label{ass:7} For all $j\in[p]$, there exists a constant $M>0$ such that
\begin{itemize}
\item[(i)] $\mathbb{E}\left(e_{j,it}|\mu_{j,it} \right)=0$,

\item[(ii)] $\mathbb{E}\left[\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left[
e_{j,it}e_{j,is}-\mathbb{E}(e_{j,it}e_{j,is})\right]\right]^{2}\leq M, $

\item[(iii)] for all $i\in[N]$, $\frac{1}{T}\sum_{t=1}^{T}\sum_{s=1}^{T}
\left\vert\mathbb{E}(e_{j,it}e_{j,is})\right\vert\leq M$,

\item[(iv)] $\max_{t\in[T]}\frac{1}{N\sqrt{T}}\left\Vert
e_{j,t}^{\prime}E_{j}\right\Vert_{2}=O_{p}\left(\frac{\log N\vee T}{N\wedge T
}\right)$, and for $\max_{i\in [N]}\frac{1}{T\sqrt{N}}\left\Vert
e_{j,i}^{\prime}E_{j}^{\prime}\right\Vert_{2}=O_{p}\left(\frac{\log N\vee T}{
N\wedge T}\right)$, where $e_{j,i}=\left(e_{j,i1},...,e_{j,iT}\right)
^{\prime}$, $e_{j,t}=\left(e_{j,1t},...,e_{j,Nt}\right) ^{\prime}$, and $
E_{j}=\{e_{j,it}\}_{i \in [N],t \in [T]}$.
\end{itemize}
\end{ass}

\begin{ass}
\label{ass:8} For all $j\in[p]$,

\begin{itemize}
\item[(i)] recall that $L_{j}^{0}=\left(l_{j,1}^{0},\cdots,l_{j,N}^{0}
\right) ^{\prime}$ and $W_{j}^{0}=\left(w_{j,1}^{0},\cdots,w_{j,T}^{0}
\right) ^{\prime}$. $\lim_{N\rightarrow \infty }\frac{L_{j}^{0
\prime}L_{j}^{0}}{N}=\Sigma_{L_{j}}>0$ and $\lim_{T\rightarrow\infty }\frac{
W_{j}^{0\prime}W_{j}^{0}}{T}=\Sigma_{W_{j}}>0$,

\item[(ii)] the $r_{j}$ eigenvalues of $\Sigma_{L_{j}}\Sigma_{W_{j}}$ are
distinct.
\end{itemize}
\end{ass}

\begin{ass}
\label{ass:9} For all $j\in[p]$, there exists a constant $M>0$ such that

\begin{itemize}
\item[(i)] $\max_{t\in[T]}\mathbb{E}\left\Vert \frac{1}{\sqrt{N}}
\sum_{i=1}^{N}l_{j,i}^{0}e_{j,it}\right\Vert_{2}^{2}\leq M$ and $\max_{t\in[T
]}\frac{1}{NT}e_{j,t}^{\prime}E_{j}^{\prime}L_{j}^{0}=O_{p}\left(\frac{\log
\left(N\vee T\right) }{N\wedge T}\right) $,

\item[(ii)] $\max_{i\in [N]}\mathbb{E}\left\Vert \frac{1}{\sqrt{T}}
\sum_{t=1}^{T}w_{j,t}^{0}e_{j,it}\right\Vert_{2}^{2}\leq M$ and $\max_{i\in[N
]}\frac{1}{NT}e_{j,t}^{\prime}E_{j}^{\prime}W_{j}^{0}=O_{p}\left(\frac{\log
\left(N\vee T\right) }{N\wedge T}\right) $.
\end{itemize}
\end{ass}

\begin{ass}
\label{ass:10} For $\forall j\in[p]$,

\begin{itemize}
\item[(i)] $\mathbb{E}\left[f_{it}(0)e_{j,it} \bigg|\mathscr{D}\right]=0$,

\item[(ii)] for each $i\in[N]$ and $j\in[p]$, $\left
\{f_{it}(0),f_{it}(0)e_{j,it}\right\}$ is stationary strong mixing across $t$
conditional on $\mathscr{D}$.
\end{itemize}
\end{ass}

Assumptions \ref{ass:7}-\ref{ass:9} are stronger than those in \cite
{bai2020simpler} because we strengthen their Assumptions A1(c) and A3 to hold uniformly. Assumption \ref{ass:10} imposes some moment
and mixing conditions. Even though $f_{it}(\cdot )$ (the PDF of $\epsilon _{it}$ given $\mathscr{D}
_{e}$) is a function of $
\left\{ e_{j,it}\right\} _{j\in \lbrack p],i\in \lbrack N],t\in \lbrack T]}$
, we can show that Assumption \ref{ass:10} holds under some reasonable
conditions. For example, we consider the location scale model:
\begin{equation*}
Y_{it}=\beta _{0,it}+\sum_{j\in \lbrack p]}X_{j,it}\beta _{j,it}+\left(
\gamma _{0,it}+\sum_{j\in \lbrack p]}X_{j,it}\gamma _{j,it}\right)
u_{it},\quad \text{with}\quad X_{j,it}=l_{j,i}^{0\prime
}w_{j,t}^{0}+e_{j,it},
\end{equation*}
where $u_{it}$ is independent of $\{w_{j,t}^{0},e_{j,it}\}_{j\in \lbrack
p],t\in \lbrack T]}$ and $l_{j,i}^{0}$ and $\beta _{j,it}$ are nonrandom. In
this case, $\Theta _{j,it}^{0}=\beta _{j,it}+\gamma _{j,it}\mathscr{Q}_{\tau
}(u_{it})$, $\epsilon _{it}=\left( \gamma _{0,it}+\sum_{j\in \lbrack
p]}X_{j,it}\gamma _{j,it}\right) \left[ u_{it}-\mathscr{Q}_{\tau }(u_{it})
\right] $, where $\mathscr{Q}_{\tau }(u_{it})$ is the $\tau $-quantile of $
u_{it}$. It is clear that $f_{it}(\cdot )$ is the function of $\left\{
e_{j,it}\right\} _{j\in \lbrack p],i\in \lbrack N],t\in \lbrack T]}$ and all
factors. However, if $u_{it}$ is independent of sequence $\left\{
e_{j,it}\right\} _{j\in \lbrack p],t\in \lbrack T]}$, we observe that $
f_{it}(0)$ is the PDF of $u_{it}-\mathscr{Q}_{\tau }(u_{it})$ evaluated at
zero point, which is independent of $\left\{ e_{j,it}\right\} _{j\in \lbrack
p],t\in \lbrack T]}$. Therefore, Assumption \ref{ass:10}
(i) holds under mild conditions that $u_{it}$ is independent of the sequence $
\left\{ e_{it}\right\} _{t\in \lbrack T]}$ and $\mathbb{E}\left( e_{it}\big|
\mathscr{D}\right) =0$.

Define
\begin{align*}
& \hat{V}_{u_{j},i}=\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[
f_{it}(0)e_{j,it}^{2}\big|\mathscr{D}\right] v_{t,j}^{0}v_{t,j}^{0\prime},
\quad V_{u_{j},i}=\mathbb{E}\left(\hat{V}_{u_{j}}\right) \\
& \Omega_{u_{j},i}=Var\left[ \frac{1}{\sqrt{T}}
\sum_{t=1}^{T}e_{j,it}v_{t,j}^{0}(\tau -\mathbf{1}\left\{ \epsilon_{it}\leq
0\right\} )\right] , \\
& \hat{V}_{v_{j},t}^{(3)}=\frac{1}{N_{3}}\sum_{i\in
I_{3}}f_{it}(0)e_{j,it}^{2}u_{i,j}^{0}u_{i,j}^{0\prime},\quad V_{v_{j}}=
\mathbb{E}\left(\hat{V}_{v_{j},t}^{(3)}\right),
\quad\Omega_{v_{j}}^{(3)}=\tau \left(1-\tau \right) \frac{1}{N_{3}}
\sum_{i\in I_{3}}\mathbb{E}\left(e_{j,it}^{2}u_{i,j}^{0}u_{i,j}^{0\prime}
\right) .
\end{align*}
Let $\Sigma_{u_{j},i}=O_{j}^{(1)}V_{u_{j},i}^{-1}
\Omega_{u_{j},i}V_{u_{j},i}^{-1}O_{j}^{(1)\prime}$, $
\Sigma_{v_{j}}^{(3)}=O_{j}^{(1)}\left(V_{v_{j}}^{(3)}\right)^{-1}
\Omega_{v_{j}}\left(V_{v_{j}}^{(3)}\right)^{-1}O_{j}^{(1)\prime}$, $
b_{j,it}^{0}=e_{j,it}v_{t,j}^{0}(\tau -\mathbf{1}\left\{ \epsilon_{it}\leq
0\right\} )$ and $\xi_{j,it}^{0}=e_{j,it}u_{i,j}^{0}\left(\tau -\mathbf{1}
\left\{ \epsilon_{it}\leq 0\right\} \right) $. The following theorem establishes
the asymptotic properties of the third-stage estimators.

\begin{theorem}
\label{Thm3} Suppose that Assumptions \ref{ass:1}-\ref{ass:10} hold. Suppose
that Assumption \ref{ass:15} in Appendix B.3 of the online supplement hold.
Let $O_{u,j}^{(1)}$ be the bounded matrix defined in the appendix that is
related to rotation matrix $O_{j}^{(1)}$. Then we have that $\forall j\in
\lbrack p]$,

\begin{itemize}
\item[(i)] $\hat{u}_{i,j}^{(3,1)}-O_{u,j}^{(1)}u_{i,j}^{0}=O_{j}^{(1)}\hat{V}
_{u_{j},i}^{-1}\frac{1}{T}\sum_{t=1}^{T}b_{j,it}^{0}+\mathcal{R}_{i,u}^{j}$
and $\sqrt{T}\left(\hat{u}_{i,j}^{(3,1)}-O_{u,j}^{(1)}u_{i,j}^{0}\right)
\rightsquigarrow \mathcal{N}\left(0,\Sigma_{u_{j},i}\right) $ $\forall i\in
I_{3}$,

\item[(ii)] $\hat{v}_{t,j}^{(3,1)}-\left( O_{u,j}^{(1)\prime }\right)
^{-1}v_{t,j}^{0}=O_{j}^{(1)}\left( \hat{V}_{v_{j},t}^{(3)}\right) ^{-1}\frac{
1}{N_{3}}\sum_{i\in I_{3}}\xi _{j,it}^{0}+\mathcal{R}_{t,v}^{j}$ and $\sqrt{
N_{3}}\left( \hat{v}_{t,j}^{(3,1)}-O_{v,j}^{(1)}v_{t,j}^{0}\right)
\rightsquigarrow \mathcal{N}\left( 0,\Sigma _{v_{j}}^{(3)}\right) $ $\forall
t\in \lbrack T]$,\newline
where $\max_{i\in I_{3}}\left\vert \mathcal{R}_{i,u}^{j}\right\vert
=o_{p}\left( \left( N\vee T\right) ^{-1/2}\right) $, and $\max_{t\in \lbrack
T]}\left\vert \mathcal{R}_{t,v}^{j}\right\vert =o_{p}\left( \left( N\vee
T\right) ^{-\frac{1}{2}}\right) $.
\end{itemize}
\end{theorem}

\textbf{Remark 3.} Theorem \ref{Thm3} reports the linear expansions for the
factor and factor loading estimators for each slope matrix obtained in Step
3. Compared with \cite{chernozhukov2019inference}, Theorem \ref{Thm3}
obtains the uniform convergence rate rather than the point-wise result for
the reminder terms $\mathcal{R}_{i,u}^{j}$ and $\mathcal{R}_{t,v}^{j}.$ In
addition, since the regressors in the debiasing step are obtained from Step
2 instead of Step 1, we don't have independence between the regressors and
error terms, which makes the proof more complex than that in \cite
{chernozhukov2019inference}. See the proof in the appendix on how to handle
the dependence. Assumption \ref{ass:15} in the online supplement is a regularity condition on the density of $\epsilon_{it}$.

Following Theorem \ref{Thm3} and estimators defined in Table \ref
{tab:estimator symbol}, we have that $\forall j\in \left[ p\right] $, $
\forall i\in \lbrack N]$ and $\forall t\in \left[ T\right] $,
\begin{align*}
& \hat{u}_{i,j}^{(a,b)}-O_{u,j}^{(b)}u_{i,j}^{0}=O_{j}^{(b)}\hat{V}
_{u_{j},i}^{-1}\frac{1}{T}\sum_{t=1}^{T}v_{t,j}^{0}b_{j,it}^{0}+\mathcal{R}
_{i,u}^{j}, \\
& \hat{v}_{t,j}^{(a,b)}-\left( O_{u,j}^{(b)\prime }\right)
^{-1}v_{t,j}^{0}=O_{j}^{(b)}\left( \hat{V}_{v_{j},t}^{(a)}\right) ^{-1}\frac{
1}{N_{a}}\sum_{i\in I_{a}}u_{i,j}^{0}\xi _{j,it}^{0}+\mathcal{R}_{t,v}^{j},
\end{align*}
where $\hat{V}_{v_{j},t}^{(a)}=\frac{1}{N_{a}}\sum_{i\in
I_{a}}f_{it}(0)e_{j,it}^{2}u_{i,j}^{0}u_{i,j}^{0\prime }$, $a\in \lbrack 3]$
and $b\in \lbrack 3]\setminus \{a\}$.

Given the above estimates for the factors and factor loadings, we can
estimate $\Theta _{j,it}^{0}$ by
\begin{equation*}
\hat{\Theta}_{j,it}=\frac{1}{2}\sum_{a\in \lbrack 3]}\sum_{b\in \lbrack
3]\setminus \{a\}}\left\{ \hat{u}_{i,j}^{(a,b)\prime }\hat{v}
_{t,j}^{(a,b)}\right\} \mathbf{1}_{ia}
\end{equation*}
where $\mathbf{1}_{ia}=\mathbf{1}\left\{ i\in I_{a}\right\} $ for $i\in
\lbrack N]$. Let $\Xi _{j,it}^{0}=\frac{1}{T}v_{t,j}^{0\prime }\Sigma
_{u_{j},i}v_{t,j}^{0}+\sum_{a=1}^{3}\frac{1}{N_{a}}\mathbf{1}
_{ia}u_{i,j}^{0\prime }\Sigma _{v_{j}}^{a}u_{i,j}^{0}$. The following
proposition studies the asymptotic properties of $\hat{\Theta}_{j,it}$.

\begin{proposition}
\label{Pro4} Under Assumptions \ref{ass:1}-\ref{ass:10} and Assumption \ref
{ass:15}, $\forall j\in \lbrack p]$ we have

\begin{itemize}
\item[(i)] $\hat{\Theta}_{j,it}-\Theta_{j,it}^{0}=\sum_{a=1}^{3}u_{i,j}^{0
\prime}\left(\hat{V}_{v_{j},t}^{(a)}\right) ^{-1}\frac{1}{N_{a}}
\sum_{i^{*}\in I_{a}}\xi_{j,i^{*}t}\mathbf{1}_{i^{*}a}+v_{t,j}^{0\prime}\hat{
V}_{u_{j}}^{-1}\frac{1}{T}\sum_{t^{*}=1}^{T}b_{j,it^{*}}^{0}+\mathcal{R}
_{it}^{j}$, where $\max_{i\in I_{3},t\in[T]}\left\vert \mathcal{R}_{it}^{j}\right\vert
=o_{p}\left(\left(N\vee T\right) ^{-1/2}\right)$,

\item[(ii)] $\max_{i\in[N],t\in[T]}\left\vert \hat{\Theta}
_{j,it}-\Theta_{j,it}^{0}\right\vert =O_{p}\left(\sqrt{\frac{\log N\vee T}{
N\wedge T}}\right)$,

\item[(iii)] $\left(\Xi_{j,it}^{0}\right) ^{-1/2}\left(\hat{\Theta}
_{j,it}-\Theta_{j,it}^{0}\right) \rightsquigarrow \mathcal{N}
\left(0,1\right) $.
\end{itemize}

\end{proposition}

\textbf{Remark 4.} Proposition \ref{Pro4} establishes the distribution
theory for the slope estimators. Recall that we remove the principle
component from the independent variables $X_{j,it}$ which is the key point
in the debiasing step and why we don't have the distribution theory
result for the intercept estimates $\hat{\Theta}_{0,it}$ in the current
framework. However, once we have the distribution theory for the slope
estimates, we can follow \cite{chen2021quantile} and obtain a new estimator for  $\Theta_{0,it}^0$
from the smoothed quantile regression and establish its distribution theory. We leave this for the further research.

To make inference for $u_{i,j}^{0},$ $v_{t,j}^{0},$ and $\Theta _{j,it}^{0},$
one needs to estimate their asymptotic variances $\Sigma _{u_{j},i},$ $\Sigma
_{v_{j}}$ and $\Xi _{j,it}^{0}$ consistently. Let $k(\cdot )$ be a PDF-type
kernel function and $K(\cdot )$ be its survival function such that $\int
k(u)du=1$ and $K(u):=\int_{u}^{\infty }k(v)dv$. Let $h_{N}$ be the bandwidth
such that $h_{N}\rightarrow 0$ with $N\rightarrow \infty $. Define $
K_{h_{N}}(\cdot )=K(\frac{\cdot }{h_{N}})$, $k_{h_{N}}(\cdot )=\frac{1}{h_{N}
}k(\frac{\cdot }{h_{N}})$. Let $\hat{\epsilon}_{it}=Y_{it}-\hat{\Theta}
_{0,it}-\sum_{j\in \lbrack p]}X_{j,it}\hat{\Theta}_{j,it}$, $\hat{\mathrm{v}}
_{t,s,j}=\frac{1}{6}\sum_{a\in \lbrack 3]}\sum_{b\in \lbrack 3]\setminus
\{a\}}\hat{v}_{t,j}^{(a,b)}\hat{v}_{s,j}^{(a,b)\prime }$, and $\hat{\mathrm{u
}}_{i,i,j}=\frac{1}{2}\sum_{a\in \lbrack 3]}\sum_{b\in \lbrack 3]\setminus
\{a\}}\hat{u}_{i,j}^{(a,b)}\hat{u}_{i,j}^{(a,b)\prime }\mathbf{1}_{ia}$.
Define
\begin{align*}
\hat{\mathbb{V}}_{u_{j}}& =\frac{1}{NT}\sum_{i\in \lbrack N]}\sum_{t\in
\lbrack T]}k_{h_{N}}(\hat{\epsilon}_{it})\hat{e}_{j,it}^{2}\hat{\mathrm{v}}
_{t,t,j},\quad \hat{\mathbb{V}}_{v_{j}}=\frac{1}{NT}\sum_{i\in \lbrack
N]}\sum_{t\in \lbrack T]}k_{h_{N}}(\hat{\epsilon}_{it})\hat{e}_{j,it}^{2}
\hat{\mathrm{u}}_{i,i,j}, \\
\hat{\Omega}_{u_{j}}& =\frac{1}{NT}\sum_{i\in \lbrack N]}\left\{ \sum_{t\in
\lbrack T]}\tau (1-\tau )\hat{e}_{j,it}^{2}\hat{\mathrm{v}}
_{t,t,j}+\sum_{t=1}^{T-T_{1}}\sum_{s=t+1}^{t+T_{1}}S_{j,its}+
\sum_{t=1+T_{1}}^{T}\sum_{s=t-T_{1}}^{t-1}S_{j,its}\right\} , \\
\hat{\Omega}_{v_{j}}& =\frac{\tau (1-\tau )}{NT}\sum_{i\in \lbrack
N]}\sum_{t\in \lbrack T]}\hat{e}_{j,it}^{2}\hat{\mathrm{u}}_{i,i,j},\quad
\hat{\Sigma}_{u_{j}}=\hat{\mathbb{V}}_{u_{j}}^{-1}\hat{\Omega}_{u_{j}}\hat{
\mathbb{V}}_{u_{j}}^{-1},\quad \hat{\Sigma}_{v_{j}}=\hat{\mathbb{V}}
_{v_{j}}^{-1}\hat{\Omega}_{v_{j}}\hat{\mathbb{V}}_{v_{j}}^{-1},
\end{align*}
where $S_{j,its}=\hat{e}_{j,it}\hat{e}_{j,is}\hat{\mathrm{v}}_{t,s,j}\left[
\tau -K\left( \frac{\hat{\epsilon}_{it}}{h_{N}}\right) \right] \left[ \tau
-K\left( \frac{\hat{\epsilon}_{is}}{h_{N}}\right) \right] .$ We further
define
\begin{equation*}
\hat{\Xi}_{j,it}=\frac{1}{2}\sum_{a\in \lbrack 3]}\sum_{b\in \lbrack
3]\setminus \{a\}}\left( \frac{1}{T}\hat{v}_{t,j}^{(a,b)\prime }\hat{\Sigma}
_{u_{j}}\hat{v}_{t,j}^{(a,b)}+\frac{1}{N_{a}}\mathbf{1}_{ia}\hat{u}
_{i,j}^{(a,b)\prime }\hat{\Sigma}_{v_{j}}\hat{u}_{i,j}^{(a,b)}\right) .
\end{equation*}
Let $F_{i,ts}(\cdot ,\cdot )$ and $f_{i,ts}(\cdot ,\cdot )$ denote the joint
CDF and PDF of $(\epsilon _{it},\epsilon _{is})$ given $\mathscr{D}_{e},$
respectively. To justify the consistency of the variance estimators, we add
the following assumption.

\begin{ass}
\label{ass:11}

\begin{itemize}
\item[(i)] $\int_{-\infty}^{+\infty}k(u)du=1$, $\int_{-\infty}^{+
\infty}k(u)u^{j}du=0$ for $j\in\{1,\cdots,m-1\}$ and $\int_{-\infty}^{+
\infty}k(u)u^{m}du\ne0$ for $m\geq 1$.

\item[(ii)] $h_{N}\rightarrow 0$ and $\left(\frac{\log (N\vee T)}{N\wedge T}
\right) ^{1/4}\frac{\xi_{N}^{2}}{h_{N}}\rightarrow 0$.

\item[(iii)] $T_{1}\rightarrow \infty $ and $\sqrt{\frac{\log (N\vee T)}{
N\wedge T}}\frac{\xi_{N}^{2 }T_{1}}{h_{N}^{2}}\rightarrow 0$.

\item[(iv)] $f_{it}(c)$ is $m$ times continuously differentiable with
respect to $c$ and $f_{i,ts}(c_{1},c_{2})$ is $m$ times continuously
differentiable with respect to $(c_{1},c_{2})$.

\item[(v)] $\forall i\in \lbrack N]$, $V_{u_{j},i}=V_{u_{j}}$ and $\Omega
_{u_{j},i}=\Omega _{u_{j}}$.

\item[(vi)] $\forall a\in \lbrack 3]$, $V_{v_{j}}^{(a)}=\frac{1}{N}
\operatornamewithlimits{\sum}\limits_{i\in \lbrack N]}\mathbb{E}\left[
f_{it}(0)e_{j,it}^{2}u_{i,j}^{0}u_{i,j}^{0\prime }\right] +o_{p}(1)$ and $
\Omega _{v_{j}}^{(a)}=\frac{\tau \left( 1-\tau \right) }{N}
\operatornamewithlimits{\sum}\limits_{i\in \lbrack N]}\mathbb{E}\left(
e_{j,it}^{2}u_{i,j}^{0}u_{i,j}^{0\prime }\right) +o_{p}(1)$.
\end{itemize}
\end{ass}

Assumption \ref{ass:11}(i)-(iv) are standard for consistent estimation of
the asymptotic variance matrix; see, e.g., \cite{chen2019two} and \cite
{galvao2016smoothed}. Assumption \ref{ass:11}(v) imposes the homogeneity
moment condition across individuals, and Assumption \ref{ass:11}(vi) assumes
the moments calculated from subsamples are close to those from the full
sample given the random splitting. Under Assumption \ref{ass:11}, following
the idea of \cite{chen2019two}, we establish in Lemma \ref{Lem:covhat} of
the online supplement the consistency of $\hat{\Sigma}_{u_{j}}$ and $\hat{
\Sigma}_{v_{j}}$. Similar conclusions hold for the other estimates.


\section{Specification Tests}

In this section, we consider two specification tests under different rank
conditions.

\subsection{Testing for Homogeneity across Individuals or Time}

When $K_{j}=1$ for some $j\in [p]$, it is interesting to test whether the
matrix $\Theta_{j}^{0}$ is homogeneous across individuals (i.e., row-wise) or
across time (i.e., column-wise). For these two cases, we can write factors and
factor loadings as
\begin{equation*}
u_{i,j}^{0}=u_{j}+c_{i,j}^{u}\quad \text{and}\quad v_{t,j}^{0}=v_{j}+c_{t,j}^{v}, \quad \text{respectively},
\end{equation*}
where $u_{j}=\frac{1}{N}\sum_{i=1}^{N}u_{i,j}^{0}$ and $v_{j}=\frac{1}{T}
\sum_{t=1}^{T}v_{t,j}^{0}.$ For the homogeneity across individuals, the null
and alternative hypotheses can be written as
\begin{equation}
H_{0}^{I}:c_{i,j}^{u}=0\quad \forall i\in [N]\quad v.s.\quad
H_{1}^{I}:c_{i,j}^{u}\neq 0\text{ for some }i\in [N].  \label{Hypo1}
\end{equation}
Similarly, for the homogeneity across time, the null and alternative
hypotheses can be written as
\begin{equation}
H_{0}^{II}:c_{t,j}^{v}=0\quad \forall t\in[T]\quad v.s.\quad
H_{1}^{II}:c_{t,j}^{v}\neq 0\text{ for some }t\in[T].  \label{Hypo2}
\end{equation}
Note that we aim to test the two null hypotheses separately. That is, we can
test for homogeneous slope across individuals while allowing for
heterogeneous slopes across time and vice versa. This is different from the
majority of the literature which either tests for slope homogeneity across
individuals while assuming the slopes are homogeneous across time or tests
for structural breaks across time while assuming the slopes are homogeneous
across individuals.


We first consider testing $H_{0}^{I}$. Following the lead of \cite
{castagnetti2015inference}, we define\footnote{
Alternatively, we can also define $S_{u_{j}}^{o}=\max \left(
S_{u_{j}}^{(3,2)},S_{u_{j}}^{(2,1)},S_{u_{j}}^{(1,3)}\right) .$ It is easy
to show that this statistic shares the same asymptotic null distribution as $
S_{u_{j}}.$ But due to the unknown dependence structure between the two, we
cannot take the maximum or the other continuous function of $S_{u_{j}}$ and $
S_{u_{j}}^{o}$ as a new test statistic.}
\begin{align}
& S_{u_{j}}^{(a,b)}=\max_{i\in I_{a}}T(\hat{u}_{i,j}^{(a,b)}-\hat{\bar{u}}
_{j}^{(a,b)})^{\prime }\hat{\Sigma}_{u_{j}}^{-1}(\hat{u}_{i,j}^{(a,b)}-\hat{
\bar{u}}_{j}^{(a,b)})\quad \text{and}  \notag \\
& S_{u_{j}}=\max \left(
S_{u_{j}}^{(3,1)},S_{u_{j}}^{(2,3)},S_{u_{j}}^{(1,2)}\right) ,
\label{statistic for u}
\end{align}
where $\hat{\bar{u}}_{j}^{(a,b)}=\frac{1}{N_{a}}\sum_{i\in I_{a}}\hat{u}
_{i,j}^{(a,b)}$. Similarly, to test for $H_{0}^{II}$, we construct
\begin{equation*}
S_{v_{j}}=\max \left( \tilde{S}_{v_{j}}^{(3,1)},\tilde{S}_{v_{j}}^{(2,3)},
\tilde{S}_{v_{j}}^{(1,3)}\right) ,
\end{equation*}
where
\begin{equation*}
S_{v_{j}}^{(a,b)}=\max_{t\in \lbrack T]}N(\hat{v}_{t,j}^{(a,b)}-\hat{\bar{v}}
_{j}^{(a,b)})^{\prime }\hat{\Sigma}_{v_{j}}^{-1}(\hat{v}_{t,j}^{(a,b)}-\hat{
\bar{v}}_{j}^{(a,b)}),\quad  \tilde{S}_{v_{j}}^{(a,b)}=\frac{1
}{2}S_{v_{j}}^{(a,b)}-\mathsf{b}(T),
\end{equation*}
$\hat{\bar{v}}_{j}^{(a,b)}=\frac{1}{T}\sum_{t=1}^{T}\hat{v}
_{t,j}^{(a,b)},$ and $\mathsf{b}(n)=\log n-\frac{1}{2}\log \log n-\log
\Gamma (\frac{1}{2})$ for $n\in \left\{ N,T,NT\right\} .$


To proceed, we introduce some notation. Recall that $
b_{j,it}^{0}=e_{j,it}v_{t,j}^{0}(\tau -\mathbf{1}\left\{ \epsilon_{it}\leq
0\right\} )$ and $\xi_{j,it}^{0}=e_{j,it}u_{i,j}^{0}\left(\tau -\mathbf{1}
\left\{ \epsilon_{it}\leq 0\right\} \right) $. Define
\begin{equation*}
\mathfrak{b}_{j,it}^{(1)}=\hat{V}_{u_{j}}^{-1}b_{j,it}^{0},\quad \mathfrak{b}
_{j,it}^{(2)}=\left(\hat{V}_{v_{j},t}^{(3)}\right) ^{-1}\xi_{j,it}^{0},\text{
}\mathfrak{b}_{j,it}^{(3)}=\left(\hat{V}_{v_{j},t}^{(2)}\right)^{-1}
\xi_{j,it}^{0},\text{ and }\mathfrak{b}_{j,it}^{(4)}=\left(\hat{V}
_{v_{j},t}^{(1)}\right) ^{-1}\xi_{j,it}^{0}.
\end{equation*}
Let $\mathfrak{B}_{j,t}^{(\ell )}=\left(\mathfrak{b}_{j,1t}^{(\ell
)\prime},\cdots ,\mathfrak{b}_{j,Nt}^{(\ell )\prime }\right) ^{\prime}$ for $
\ell \in \left[ 4\right] .$ Define
\begin{align*}
& \Sigma_{\mathfrak{B},j}^{(1)}=\frac{1}{T}\sum_{t=1}^{T}\sum_{s=1}^{T}
\mathbb{E}\left(\mathfrak{B}_{j,t}^{(1)}\mathfrak{B}_{j,s}^{(1)\prime}
\right) ,\quad \Sigma_{\mathfrak{B},j}^{(2)}=\frac{1}{N_{3}}\sum_{i\in I_{3}}
\mathbb{E}\left(\mathfrak{B}_{j,i}^{(2)}\mathfrak{B}_{j,i}^{(2)\prime}
\right) , \\
& \Sigma_{\mathfrak{B},j}^{(3)}=\frac{1}{N_{2}}\sum_{i\in I_{2}}\mathbb{E}
\left(\mathfrak{B}_{j,i}^{(3)}\mathfrak{B}_{j,i}^{(3)\prime }\right) ,\quad
\text{and }\Sigma_{\mathfrak{B},j}^{(4)}=\frac{1}{N_{1}}\sum_{i\in I_{1}}
\mathbb{E}\left(\mathfrak{B}_{j,i}^{(4)}\mathfrak{B}_{j,i}^{(4)\prime}
\right) .
\end{align*}
We add the following two assumptions. \bigskip

\begin{ass}
\label{ass:12} $\forall j\in \lbrack p]$, we assume
\begin{equation*}
\bar{\lambda}\geq \lambda _{\max }\left( \Sigma _{u_{j}}\right) \geq \lambda
_{\min }\left( \Sigma _{u_{j}}\right) \geq \underline{\lambda }>0,\quad \bar{
\lambda}\geq \lambda _{\max }\left( \Sigma _{v_{j}}\right) \geq \lambda
_{\min }\left( \Sigma _{v_{j}}\right) \geq \underline{\lambda }>0.
\end{equation*}
\end{ass}

\begin{ass}
\label{ass:13}
\begin{itemize}
\item[(i)] There exists a high dimensional Gaussian vector $\mathbb{Z}_{
\mathfrak{B}}^{(1)}\sim N\left(0,\Sigma_{\mathfrak{B},j}^{(1)}\right) $ such
that $\left\Vert \frac{1}{\sqrt{T}}\sum_{t=1}^{T}\mathfrak{B}_{j,t}^{(1)}-
\mathbb{Z}_{\mathfrak{B}}^{(1)}\right\Vert_{\max }=o_{p}(1).$

\item[(ii)] There exists high dimensional Gaussian vectors $\mathbb{Z}_{
\mathfrak{B}}^{(\ell )}\sim N\left(0,\Sigma_{\mathfrak{B},j}^{(\ell)}\right)
$ for $\ell =2,3,4$ such that $\left(\mathbb{Z}_{
\mathfrak{B}}^{(2 )}, \mathbb{Z}_{
\mathfrak{B}}^{(3 )}, \mathbb{Z}_{
\mathfrak{B}}^{(4)}\right)$ are independent,
\begin{align*}
& \left\Vert \frac{1}{\sqrt{N_{3}}}\sum_{i\in I_{3}}\mathfrak{B}_{j,i}^{(2)}-
\mathbb{Z}_{\mathfrak{B}}^{(2)}\right\Vert_{\max }=o_{p}(1),\quad\left\Vert
\frac{1}{\sqrt{N_{2}}}\sum_{i\in I_{2}}\mathfrak{B}_{j,i}^{(3)}-\mathbb{Z}_{
\mathfrak{B}}^{(3)}\right\Vert_{\max }=o_{p}(1),\text{ and} \\
& \left\Vert \frac{1}{\sqrt{N_{1}}}\sum_{i\in I_{1}}\mathfrak{B}_{j,i}^{(4)}-
\mathbb{Z}_{\mathfrak{B}}^{(4)}\right\Vert_{\max }=o_{p}(1).
\end{align*}
\end{itemize}
\end{ass}

Assumption \ref{ass:12} implies that both $\Sigma _{u_{j}}$ and $
\Sigma_{v_{j}}$ are well behaved. Assumption \ref{ass:13} imposes that we
can approximate high dimensional vectors $\frac{1}{\sqrt{T}}\sum_{t\in[T]}
\mathfrak{B}_{j,t}^{(1)}$, $\frac{1}{\sqrt{N_{3}}}\sum_{i\in I_{3}}\mathfrak{
B}_{j,i}^{(2)}$, $\frac{1}{\sqrt{N_{2}}}\sum_{i\in I_{2}}\mathfrak{B}
_{j,i}^{(3)}$ and $\frac{1}{\sqrt{N_{1}}}\sum_{i\in I_{1}}\mathfrak{B}
_{j,i}^{(4)}$ by four Gaussian vectors. Similar conditions have been imposed
in the literature; see, e.g., Assumption SA3 \cite{lu2021uniform}.

The following theorem reports the asymptotic properties of $S_{u_{j}}$ and $
S_{v_{j}}$ under the respective null and alternative hypotheses.

\begin{theorem}
\label{Thm5} Suppose that Assumptions \ref{ass:1}-\ref{ass:13} and
Assumptions \ref{ass:15} in the online supplement hold and $(N,T)\rightarrow \infty $. Then

\begin{itemize}
\item[(i)] Under $H_{0}^{I}$, we have $\mathbb{P}\left(\frac{1}{2}
S_{u_{j}}\leq x+\mathsf{b}(N)\right) \rightarrow e^{-e^{-x}}$; and under $
H_{0}^{II}$, we have $\mathbb{P}\left(S_{v_{j}}\leq x\right) \rightarrow
e^{-3e^{-x}}.$

\item[(ii)] Under $H_{1}^{I}$, if $\frac{T}{\log N}\max_{i\in [N]}\left\Vert
c_{i,j}^{u}\right\Vert_{2}^{2}\rightarrow \infty $, we have $\mathbb{P}
\left(S_{u_{j}}>c_{\alpha ,1\cdot N}\right) \rightarrow 1$ with $c_{\alpha
,1\cdot N}=2\mathsf{b}(N)-\log \left\vert \log \left(1-\alpha \right)
\right\vert ^{2}$ and $\alpha $ is the significance level. Under $H_{1}^{II}$
, if $\frac{N}{\log T}\max_{t\in[T]}\left\Vert
c_{t,j}^{v}\right\Vert_{2}^{2}\rightarrow \infty $, we have $\mathbb{P}
\left(S_{v_{j}}>c_{\alpha,2}\right)\to 1$ with $c_{\alpha,2}=-\log \left(-
\frac{1}{3}\log \left(1-\alpha\right)\right) $.
\end{itemize}
\end{theorem}

\textbf{Remark 5.} Theorem \ref{Thm5} implies that our test statistics
follow the Gumbel distributions asymptotically under the null, are
consistent under the global alternatives, and have non-trivial power against
the local alternatives. The power function of $S_{u_{j}}$ approaches 1 as
long as $\frac{T}{\log N}\max_{i\in [N]}\left\Vert
c_{i,j}^{u}\right\Vert_{2}^{2}$ diverges to infinity as $\left(N,T\right)
\rightarrow \infty .$


\subsection{Test for an Additive Structure}

When $K_{j}=2$ for some $j\in [p]$, it is interesting to test whether $
\Theta_{j,it}^{0}$ exhibits the additive structure which is widely assumed in a
two-way fixed effects model. That is, one may test the following null
hypothesis
\begin{equation}
H_{0}^{III}:\Theta_{j,it}^{0}=\lambda_{j,i}+f_{j,t},~ \forall
\left(i,t\right) \in [N]\times [T],  \label{Hypo3}
\end{equation}
The alternative hypothesis $H_{1}^{III}$ is the negation of $H_{0}^{III}.$

Let $\bar{\Theta}_{j,i\cdot }=\frac{1}{T}\sum_{t\in \lbrack T]}\Theta
_{j,it}^{0}$, $\bar{\Theta}_{j,\cdot t}^{I_{a}}=\frac{1}{N_{a}}\sum_{i\in
I_{a}}\Theta _{j,it}^{0}$, and $\bar{\Theta}_{j}^{I_{a}}=\frac{1}{N_{a}T}
\sum_{i\in I_{a}}\sum_{t\in \lbrack T]}\Theta _{j,it}^{0}$ for $a\in \lbrack
3]$. Define
\begin{equation*}
\Theta _{j,it}^{\ast }=\Theta _{j,it}^{0}-\bar{\Theta}_{j,i\cdot }-\bar{
\Theta}_{j,\cdot t}^{I_{a}}+\bar{\Theta}_{j}^{I_{a}},\quad \forall i\in
I_{a},t\in \lbrack T],j\in \lbrack p].
\end{equation*}
Note that $\Theta _{j,it}^{\ast }=0$ $\forall \left( i,t\right) \in \lbrack
N]\times \lbrack T]$ under $H_{0}^{III}.$ So we can propose a test for $
H_{0}^{III}$ based on estimates of $\Theta _{j,it}^{\ast }.$ Define
\begin{equation*}
\hat{\bar{\Theta}}_{j,i\cdot }=\frac{1}{T}\sum_{t\in \lbrack T]}\hat{\Theta}
_{j,it},\text{ }\hat{\bar{\Theta}}_{j,\cdot t}^{I_{a}}=\frac{1}{N_{a}}
\sum_{i\in I_{a}}\hat{\Theta}_{j,it},\text{ and }\hat{\bar{\Theta}}
_{j}^{I_{a}}=\frac{1}{N_{a}T}\sum_{i\in I_{a}}\sum_{t\in \lbrack T]}\hat{
\Theta}_{j,it}
\end{equation*}
for $a\in \lbrack 3]$. Then, we can define the sample analogue of $\hat{
\Theta}_{j,it}^{\ast }$ as
\begin{equation*}
\hat{\Theta}_{j,it}^{\ast a}=\hat{\Theta}_{j,it}-\hat{\bar{\Theta}}
_{j,i\cdot }-\hat{\bar{\Theta}}_{j,\cdot t}^{I_{a}}+\hat{\bar{\Theta}}
_{j}^{I_{a}},\quad \forall i\in I_{a},t\in \lbrack T],j\in \lbrack p].
\end{equation*}
Its corresponding asymptotic variance can be estimated by $\hat{\Sigma}
_{j,it}^{\ast }$ defined as
\begin{align*}
\hat{\Sigma}_{j,it}^{\ast }& =\frac{1}{2}\sum_{a\in \lbrack 3]}\sum_{b\in
\lbrack 3]\setminus \{a\}}\frac{1}{N_{a}}\left( \hat{u}_{i,j}^{(a,b)}-\hat{
\bar{u}}_{j}^{(a,b)}\right) ^{\prime }\hat{\Sigma}_{v_{j}}\left( \hat{u}
_{i,j}^{(a,b)}-\bar{u}_{j}^{(a,b)}\right) \mathbf{1}_{ia} \\
& +\frac{1}{6}\sum_{a\in \lbrack 3]}\sum_{b\in \lbrack 3]\setminus \{a\}}
\frac{1}{T}\left( \hat{v}_{t,j}^{(a,b)}-\hat{\bar{v}}_{j}^{(a,b)}\right)
^{\prime }\hat{\Sigma}_{u_{j}}\left( \hat{v}_{t,j}^{(a,b)}-\hat{\bar{v}}
_{j}^{(a,b)}\right) ,
\end{align*}
where $\hat{\bar{u}}_{j}^{(a,b)}=\frac{1}{N_{a}}\sum_{i\in I_{a}}\hat{u}
_{i,j}^{(a,b)}$ and $\hat{\bar{v}}_{j}^{(a,b)}=\frac{1}{T}\sum_{t\in \lbrack
T]}\hat{v}_{t,j}^{(a,b)}$. Then, the final test statistic is
$$S_{NT}=\max_{i\in \lbrack N],t\in \lbrack T]}\left( \hat{\Theta}
_{j,it}^{\ast }\right) ^{2}/\hat{\Sigma}_{j,it}^{\ast }.$$

The following theorem studies the asymptotic properties of $S_{NT}$ under
the null and alternatives.

\begin{theorem}
\label{Thm6} Suppose Assumptions \ref{ass:1}-\ref{ass:15} hold and $(N,T)\rightarrow \infty $. Under $H_{0}^{III}$,
\begin{equation*}
\mathbb{P}\left(\frac{1}{2}S_{NT}\leq x+\mathsf{b}(NT)\right) \rightarrow
e^{-e^{-x}};
\end{equation*}
under $H_{1}^{III}$, if $\frac{N\wedge T}{\log NT}\max_{i\in [N],t\in[T]
}\left\vert \Theta_{j,it}^{\ast }\right\vert ^{2}\rightarrow\infty $, then
we have $\mathbb{P}\left(S_{NT}>c_{\alpha ,3\cdot NT}\right) \rightarrow 1$
with $c_{\alpha ,3\cdot NT}=2\mathsf{b}(NT)-\log \left\vert \log
\left(1-\alpha \right) \right\vert ^{2}$.
\end{theorem}

Similar remark after Theorem \ref{Thm5} holds here. In particular, $
S_{NT}$ has the desired asymptotic Gumbel distribution under the null and is
consistent under the global alternative.

\section{Monte Carlo Simulations}

In this section, we conduct a set of Monte Carlo simulations to show the
finite sample performance of our low-rank quantile regression estimates and specification tests.

\subsection{Data Generating Processes}

Below we will consider the following data generating process (DGP):
\begin{equation*}
Y_{it}=\Theta _{0,it}+X_{it}^{\prime }\Theta
_{it}+(1+0.1X_{1,it}+0.1X_{2,it})u_{it},
\end{equation*}
where $X_{it}=(X_{1,it},X_{2,it})^{\prime }$, $\Theta _{it}=(\Theta
_{1,it},\Theta _{2,it})^{\prime }$, $\Theta _{0,it}$ is the intercept term
which will be specified via the IFEs.

First, we consider four DGPs where the rank of each slope matrix is 1:

\begin{itemize}[leftmargin=40pt]
\item[DGP 1:] \textbf{Constant slope with i.i.d. error.} Let $\Theta
_{0,it}=\lambda _{i}f_{t}$, where $\lambda _{i},f_{t}\sim N(2,5)$. Then let $
\Theta _{1,it}=\Theta _{2,it}=2$ $\forall \left( i,t\right) \in \lbrack
N]\times \lbrack T]$, and $X_{j,it}=l_{j,i}^{0}w_{j,t}^{0}+U(0,1)$ for $j\in
\{1,2\}$ with $l_{1,i}^{0}$, $l_{2,i}^{0}$, $w_{1,t}^{0}$ and $w_{2,t}\sim U(0,1)$. $
u_{it}\operatornamewithlimits{\sim}\limits^{i.i.d}\frac{t(3)}{\sqrt{3}}$.

\item[DGP 2:] \textbf{Factor slope with rank 1 and i.i.d. error.} Same as DGP 1 except that the slope coefficients follow the factor structure with one
factor rather than homogeneous across both individuals and time, i.e., $
\Theta _{1,it}=a_{1,i}g_{1,t}$, $\Theta _{2,it}=a_{2,i}g_{2,t}$, where $
a_{1,i}$, $g_{1,t}$, $a_{2,i}$ and $g_{2,t}\sim N(0,2)$. Except these, all
other settings remain the same as in DGP 1.

\item[DGP 3:] \textbf{Constant slope with serial correlation.} Same as DGP 1 except that we set $u_{it}=0.2u_{i,t-1}+\varepsilon _{it}$, $\varepsilon _{it}
\operatornamewithlimits{\sim}\limits^{i.i.d}\frac{t(3)}{\sqrt{3}}$ and all
other settings remain the same.

\item[DGP 4:] \textbf{Factor slope with rank 1 and serial correlation.}
Same as DGP 2 except that we set $u_{it}=0.2u_{i,t-1}+\varepsilon _{it}$, $
\varepsilon _{it}\operatornamewithlimits{\sim}\limits^{i.i.d}\frac{t(3)}{
\sqrt{3}}$ and all other settings remain the same.
\end{itemize}

For the case that the rank of the slope matrix is 2, we consider two DGPs
which have the additive structure for the slope coefficient of one regressor
and the factor structure with two factors for the slope coefficient of
another regressor. Specifically,

\begin{itemize}[leftmargin=40pt]
\item[DGP 5:] \textbf{Additive and factor slopes with i.i.d. error.} $
\Theta_{0,it}=\lambda_{i}f_{t}$, $\Theta_{1,it}=a_{1,i}+g_{1,t}$ and $
\Theta_{2,it}=a_{2,i}^{\prime}g_{2,t}$ such that $
a_{2,i}=(a_{2,i,1},a_{2,i,2})^{\prime}$, $g_{2,t}=(g_{2,t,1},g_{2,t,2})^{
\prime}$, $\lambda_{i},f_{t},a_{1,i},g_{1,i}\sim N(2,5)$ and $
a_{2,i,1},a_{2,i,2},g_{2,i,1},g_{2,i,2}\sim N(0,5)$. Moreover, $
X_{1,it}=l_{1,i}^{0}w_{1,t}^{0}+U(0,4)$, $
X_{2,it}=l_{2,i}^{0}w_{2,t}^{0}+Beta(2,5)$ with $l_{1,i}^{0},w_{1,t}^{0}\sim
U(0,4)$ and $l_{2,i}^{0},w_{2,t}^{0}\sim Beta(2,5)$. $u_{it}
\operatornamewithlimits{\sim}\limits^{i.i.d}\frac{t(3)}{\sqrt{3}}$.

\item[DGP 6:] \textbf{Additive and factor slopes with serial correlation.}
Same as DGP 5 except that the error $u_{it}$ follows AR(1) process like in DGPs 3 and 4.
\end{itemize}



\subsection{Estimation Results}


For $\Theta \in \mathbb{R}^{N\times T}$, define $RMSE(\Theta )=\frac{1}{
\sqrt{NT}}\left\Vert \Theta -\Theta ^{0}\right\Vert _{F}$. Table \ref
{tab:RMSE} shows the RMSEs of the full-sample low rank matrix estimates
under different quantiles for each DGP. As Theorem \ref{Thm1}(i) predicts,
the RMSEs decrease as both $N$ and $T\ $increase. Given the fact that $
N\wedge T=T$ in the simulations, the decrease of the RMSEs is largely driven
by the increase of $T.$

Table \ref{tab:rank} reports the frequency of correct rank estimation by the
singular value thresholding (SVT) approach based on 1000 replications. Note
that the true ranks of the intercept and slope matrices in DGPs 1-4 and
5-6 are 1 and 2, respectively. The results show that the SVT can
accurately determine the correct rank of the coefficient matrices in all
DGPs for all three quantile indices under investigation.


\begin{table}[h]
\caption{RMSEs of low rank estimates in the full sample}
\scriptsize \centering
\label{tab:RMSE}
\scalebox{0.8}{\begin{tabular}{cccccccccccc}
    \toprule
    \toprule
    \multirow{2}[4]{*}{DGP} & \multirow{2}[4]{*}{N} & \multirow{2}[4]{*}{T} & \multicolumn{3}{c}{$\tau=0.25$} & \multicolumn{3}{c}{$\tau=0.50$} & \multicolumn{3}{c}{$\tau=0.75$} \\
\cmidrule{4-12}          &       &       & $\tilde{\Theta}_{0}$ & $\tilde{\Theta}_{1}$ & $\tilde{\Theta}_{2}$ & $\tilde{\Theta}_{0}$ & $\tilde{\Theta}_{1}$ & $\tilde{\Theta}_{2}$ & $\tilde{\Theta}_{0}$ & $\tilde{\Theta}_{1}$ & $\tilde{\Theta}_{2}$ \\
    \midrule
    \multirow{4}[2]{*}{1} & \multirow{2}[1]{*}{75} & 35    & 0.922 & 0.324 & 0.329 & 1.242 & 0.288 & 0.297 & 1.839 & 0.609 & 0.658 \\
          &       & 70    & 0.707 & 0.280 & 0.275 & 0.819 & 0.220 & 0.203 & 1.266 & 0.519 & 0.523 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.012 & 0.337 & 0.340 & 1.099 & 0.258 & 0.262 & 1.932 & 0.661 & 0.623 \\
          &       & 70    & 0.745 & 0.272 & 0.265 & 0.825 & 0.205 & 0.206 & 1.324 & 0.522 & 0.504 \\
    \midrule
    \multirow{4}[2]{*}{2} & \multirow{2}[1]{*}{75} & 35    & 0.871 & 0.521 & 0.505 & 0.881 & 0.704 & 0.680 & 1.278 & 1.055 & 0.970 \\
          &       & 70    & 0.692 & 0.401 & 0.373 & 0.672 & 0.553 & 0.537 & 1.057 & 0.744 & 0.768 \\
          & \multirow{2}[1]{*}{150} & 35    & 0.877 & 0.507 & 0.480 & 1.022 & 0.790 & 0.815 & 1.334 & 1.018 & 1.040 \\
          &       & 70    & 0.703 & 0.374 & 0.373 & 0.689 & 0.531 & 0.538 & 1.059 & 0.829 & 0.787 \\
    \midrule
    \multirow{4}[2]{*}{3} & \multirow{2}[1]{*}{75} & 35    & 0.945 & 0.334 & 0.329 & 1.115 & 0.280 & 0.265 & 1.876 & 0.630 & 0.627 \\
          &       & 70    & 0.682 & 0.286 & 0.279 & 0.809 & 0.230 & 0.214 & 1.244 & 0.486 & 0.492 \\
          & \multirow{2}[1]{*}{150} & 35    & 0.973 & 0.334 & 0.331 & 1.211 & 0.287 & 0.291 & 1.771 & 0.590 & 0.612 \\
          &       & 70    & 0.757 & 0.274 & 0.272 & 0.801 & 0.208 & 0.195 & 1.360 & 0.494 & 0.527 \\
    \midrule
    \multirow{4}[2]{*}{4} & \multirow{2}[1]{*}{75} & 35    & 0.885 & 0.515 & 0.519 & 0.915 & 0.693 & 0.723 & 1.382 & 1.125 & 1.037 \\
          &       & 70    & 0.669 & 0.393 & 0.384 & 0.652 & 0.511 & 0.520 & 1.053 & 0.812 & 0.774 \\
          & \multirow{2}[1]{*}{150} & 35    & 0.889 & 0.513 & 0.483 & 0.905 & 0.761 & 0.686 & 1.409 & 1.118 & 1.133 \\
          &       & 70    & 0.725 & 0.376 & 0.377 & 0.717 & 0.547 & 0.565 & 1.058 & 0.724 & 0.775 \\
    \midrule
    \multirow{4}[2]{*}{5} & \multirow{2}[1]{*}{75} & 35    & 0.218 & 0.268 & 0.450 & 0.307 & 0.308 & 0.606 & 0.844 & 0.466 & 0.936 \\
          &       & 70    & 0.174 & 0.226 & 0.414 & 0.213 & 0.200 & 0.493 & 0.610 & 0.388 & 0.838 \\
          & \multirow{2}[1]{*}{150} & 35    & 0.236 & 0.245 & 0.458 & 0.299 & 0.291 & 0.634 & 1.299 & 0.863 & 1.778 \\
          &       & 70    & 0.174 & 0.214 & 0.423 & 0.216 & 0.203 & 0.450 & 0.629 & 0.377 & 0.679 \\
    \midrule
    \multirow{4}[2]{*}{6} & \multirow{2}[1]{*}{75} & 35    & 0.253 & 0.267 & 0.293 & 0.382 & 0.227 & 0.421 & 1.293 & 0.609 & 0.892 \\
          &       & 70    & 0.207 & 0.239 & 0.278 & 0.261 & 0.192 & 0.366 & 0.576 & 0.287 & 0.415 \\
          & \multirow{2}[1]{*}{150} & 35    & 0.225 & 0.254 & 0.269 & 0.363 & 0.225 & 0.422 & 1.486 & 0.695 & 0.992 \\
          &       & 70    & 0.193 & 0.254 & 0.263 & 0.254 & 0.171 & 0.379 & 0.797 & 0.391 & 0.551 \\
    \bottomrule
    \end{tabular}}
\end{table}

\begin{table}[h]
\caption{Frequency of correct rank estimation via the SVT approach}
\label{tab:rank}
\scriptsize \centering
\scalebox{0.8}{\begin{tabular}{cccccccccccc}
    \toprule
    \toprule
    \multirow{2}[4]{*}{DGP} & \multirow{2}[4]{*}{N} & \multirow{2}[4]{*}{T} & \multicolumn{3}{c}{$\tau=0.25$} & \multicolumn{3}{c}{$\tau=0.50$} & \multicolumn{3}{c}{$\tau=0.75$} \\
\cmidrule{4-12}          &       &       & $\hat{K}_{0}$ & $\hat{K}_{1}$ & $\hat{K}_{2}$ & $\hat{K}_{0}$ & $\hat{K}_{1}$ & $\hat{K}_{2}$ & $\hat{K}_{0}$ & $\hat{K}_{1}$ & $\hat{K}_{2}$ \\
    \midrule
    \multirow{4}[2]{*}{1} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 0.996 & 0.996 & 1.00  & 0.999 & 1.00  & 1.00  & 0.999 & 0.999 \\
          &       & 70    & 1.00  & 0.994 & 0.996 & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 0.994 & 0.995 & 1.00  & 1.00  & 0.999 & 1.00  & 1.00  & 0.999 \\
          &       & 70    & 1.00  & 0.995 & 0.996 & 1.00  & 0.999 & 1.00  & 1.00  & 0.999 & 1.00 \\
    \midrule
    \multirow{4}[2]{*}{2} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 0.993 & 0.999 & 1.00  & 0.999 & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 0.997 & 0.998 & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 0.995 & 0.996 & 1.00  & 0.998 & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 0.998 & 1.00  & 1.00  & 1.00 \\
    \midrule
    \multirow{4}[2]{*}{3} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 0.990 & 0.997 & 1.00  & 0.997 & 0.997 & 1.00  & 0.999 & 0.999 \\
          &       & 70    & 1.00  & 0.994 & 0.994 & 1.00  & 0.999 & 0.999 & 1.00  & 0.999 & 0.998 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 0.999 & 0.992 & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 0.996 & 0.994 & 1.00  & 0.998 & 1.00  & 1.00  & 1.00  & 1.00 \\
    \midrule
    \multirow{4}[2]{*}{4} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 0.992 & 0.991 & 1.00  & 0.999 & 0.999 & 1.00  & 0.999 & 1.00 \\
          &       & 70    & 1.00  & 0.995 & 0.995 & 1.00  & 0.999 & 0.999 & 1.00  & 1.00  & 1.00 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 0.996 & 0.997 & 1.00  & 0.999 & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 0.997 & 0.999 & 1.00  & 0.999 & 1.00  & 1.00  & 1.00  & 1.00 \\
    \midrule
    \multirow{4}[2]{*}{5} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
    \midrule
    \multirow{4}[2]{*}{6} & \multirow{2}[1]{*}{75} & 35    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          & \multirow{2}[1]{*}{150} & 35    & 1.00  & 1.00  & 0.999 & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
          &       & 70    & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00  & 1.00 \\
    \bottomrule
    \end{tabular}}
\end{table}

\subsection{Test Results}

In Section 4, we define $S_{u_{j}}$\ and $S_{v_{j}}$\ as the sup-type test
statistics. Table \ref{tab:test_homo} reports the empirical size and power
at the 5\% nominal level for the null hypothesis that the slope coefficient
is homogeneous across either $i$ or $t$. The results in DGPs 1 and 3 give
the empirical size, and those in DGPs 2 and 4 give the empirical power. As
the results in Table \ref{tab:test_homo} indicate, our tests have reasonable
size despite the fact that they are slightly conservative like most extreme-value
based sup-tests in the literature. In terms of power, out tests have superb
power in both DGPs across all three quantile indices.

Table \ref{tab:test_additive} shows the empirical size and power of our test
for DGPs 5 and 6. The findings are similar to those in Table \ref
{tab:test_homo}. In particular, our tests are a bit conservative under the null. The empirical power tends to 1 quickly as $T$ increases.



\begin{table}[h]
\caption{Empirical size and power of testing slope homogeneity across either
$i$ or $t$ (nominal level: $0.05$)}
\label{tab:test_homo}
\scriptsize \centering
\begin{tabular}{ccccccccccccccc}
\toprule \toprule \multirow{2}[4]{*}{DGP} & \multirow{2}[4]{*}{N} &
\multirow{2}[4]{*}{T} & \multicolumn{4}{c}{$\tau=0.25$} & \multicolumn{4}{c}{
$\tau=0.5$} & \multicolumn{4}{c}{$\tau=0.75$} \\
\cmidrule{4-15} &  &  & $u_{1}$ & $v_{1}$ & $u_{2}$ & $v_{2}$ & $u_{1}$ & $
v_{1}$ & $u_{2}$ & $v_{2}$ & $u_{1}$ & $v_{1}$ & $u_{2}$ & $v_{2}$ \\
\midrule \multirow{4}[2]{*}{DGP 1} & \multirow{2}[1]{*}{75} & 35 & 0.040 &
0.051 & 0.049 & 0.032 & 0.024 & 0.054 & 0.034 & 0.054 & 0.036 & 0.047 & 0.036
& 0.048 \\
&  & 70 & 0.040 & 0.055 & 0.050 & 0.044 & 0.020 & 0.056 & 0.017 & 0.068 &
0.025 & 0.037 & 0.029 & 0.029 \\
& \multirow{2}[1]{*}{150} & 35 & 0.028 & 0.036 & 0.058 & 0.048 & 0.065 &
0.054 & 0.052 & 0.055 & 0.074 & 0.030 & 0.076 & 0.024 \\
&  & 70 & 0.034 & 0.025 & 0.030 & 0.023 & 0.035 & 0.048 & 0.028 & 0.040 &
0.035 & 0.025 & 0.039 & 0.025 \\
\midrule \multirow{4}[2]{*}{DGP 2} & \multirow{2}[1]{*}{75} & 35 & 1.00 &
1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
\\
&  & 70 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
& 1.00 & 1.00 \\
& \multirow{2}[1]{*}{150} & 35 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 &
1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
&  & 70 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
& 1.00 & 1.00 \\
\midrule \multirow{4}[2]{*}{DGP 3} & \multirow{2}[1]{*}{75} & 35 & 0.045 &
0.057 & 0.050 & 0.041 & 0.022 & 0.050 & 0.038 & 0.089 & 0.054 & 0.047 & 0.049
& 0.047 \\
&  & 70 & 0.048 & 0.031 & 0.031 & 0.033 & 0.028 & 0.086 & 0.023 & 0.069 &
0.041 & 0.046 & 0.032 & 0.038 \\
& \multirow{2}[1]{*}{150} & 35 & 0.065 & 0.054 & 0.058 & 0.034 & 0.064 &
0.051 & 0.068 & 0.045 & 0.084 & 0.018 & 0.089 & 0.023 \\
&  & 70 & 0.046 & 0.030 & 0.044 & 0.025 & 0.022 & 0.037 & 0.037 & 0.030 &
0.046 & 0.022 & 0.048 & 0.015 \\
\midrule \multirow{4}[2]{*}{DGP 4} & \multirow{2}[1]{*}{75} & 35 & 1.00 &
1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
\\
&  & 70 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
& 1.00 & 1.00 \\
& \multirow{2}[1]{*}{150} & 35 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 &
1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
&  & 70 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00
& 1.00 & 1.00 \\
\bottomrule &  &  &  &  &  &  &  &  &  &  &  &  &  &
\end{tabular}
\end{table}

\begin{table}[h]
\caption{Empirical size and power for testing additive slopes (nomial level:
$0.05$)}
\label{tab:test_additive}
\scriptsize \centering
\begin{tabular}{ccccccccc}
\toprule \toprule \multirow{2}[4]{*}{DGP} & \multirow{2}[4]{*}{N} &
\multirow{2}[4]{*}{T} & \multicolumn{2}{c}{$\tau=0.25$} & \multicolumn{2}{c}{
$\tau=0.50$} & \multicolumn{2}{c}{$\tau=0.75$} \\
\cmidrule{4-9} &  &  & size & power & size & power & size & power \\
\midrule \multirow{4}[2]{*}{DGP 5} & \multirow{2}[1]{*}{75} & 35 & 0.027 &
0.807 & 0.034 & 1.00 & 0.026 & 0.979 \\
&  & 70 & 0.065 & 1.00 & 0.061 & 1.00 & 0.034 & 1.00 \\
& \multirow{2}[1]{*}{150} & 35 & 0.018 & 1.00 & 0.026 & 1.00 & 0.011 & 1.00
\\
&  & 70 & 0.012 & 1.00 & 0.029 & 1.00 & 0.010 & 1.00 \\
\midrule \multirow{4}[2]{*}{DGP 6} & \multirow{2}[1]{*}{75} & 35 & 0.039 &
0.796 & 0.058 & 1.00 & 0.045 & 0.979 \\
&  & 70 & 0.024 & 1.00 & 0.06 & 1.00 & 0.032 & 1.00 \\
& \multirow{2}[1]{*}{150} & 35 & 0.019 & 1.00 & 0.025 & 1.00 & 0.022 & 1.00
\\
&  & 70 & 0.01 & 1.00 & 0.022 & 1.00 & 0.014 & 1.00 \\
\bottomrule &  &  &  &  &  &  &  &
\end{tabular}
\end{table}

\section{Empirical Study}

In this section we consider two empirical applications: the heterogeneous investment equation and the heterogeneous
quantile effect of foreign direct investment on unemployment.
\subsection{Investment Equation}

In this subsection, we revisit the investment equation. \cite
{fazzari1988financing} point out that investment may show sensitivity to
movements in cash flow when firms face constraints for external finance.
Since \cite{fazzari1988financing}, there has been a large literature on the effect of cash flow on
the corporate investment; see \cite{devereux199011}, \cite
{gilchrist1995evidence}, \cite{kaplan1995financing}, \cite
{cleary1999relationship}, \cite{rauh2006investment}, and \cite
{almeida2007financial}, among others. Using the panel dataset, we consider
the scaled version of the investment equation as follows:
\begin{equation*}
\frac{I_{it}}{K_{i,t-1}}=\Theta _{0,it}+\Theta _{1,it}\frac{CF_{it}}{
K_{i,t-1}}+\Theta _{2,it}q_{i,t-1}+u_{it},
\end{equation*}
where $I$ is the corporate investment, $CF$ is the cash flow, $q$ is the
Tobin's q, $K$ is the capital stock and $u$ is the innovation. $\Theta
_{0,it}$ refers to the fixed effects (FEs). Rather than the mean estimation,
\cite{galvao2015efficient} estimate the effects of the firm's cash flow and
Tobin's q on investment at different quantiles. By using the panel quantile
regression with individual FEs, they show that the slope estimates change
across $\tau $. However, they do not allow the slope
coefficients, $\Theta _{1}$ and $\Theta _{2},$ to change either over $i$ or $
t.$ Inspired by \cite{galvao2015efficient}, we estimate the following model
\begin{equation}
\label{base}
\mathscr{Q}_{\tau }\left( IK_{it}\big |\left\{ CFK_{it},q_{i,t-1}\right\}
_{t\in \lbrack T]},\left\{ \Theta _{j,it}\right\} _{t\in \lbrack T],j\in
\{0,1,2\}}\right) =\Theta _{0,it}(\tau )+\Theta _{1,it}(\tau
)CFK_{it}+\Theta _{2,it}(\tau )q_{i,t-1},
\end{equation}
where $IK_{it}=\frac{I_{it}}{K_{i,t-1}}$, and $CFK_{it}=\frac{CF_{it}}{
K_{i,t-1}}$. Here we don't restrict the specific structure on the FEs and
they can be either additive or interactive.

The data are taken from the China Stock Market \& Accounting Research
(CSMAR) Database. We use quarterly data for 195 manufacturing firms in China
from 2003 to 2020. Based on the model (\ref{base}), we define corporate
investment as $I_{it}=LI_{it}-LI_{i,t-1}$, where $LI_{it}$ is the total
value of long-term corporate investment as the sum of long-term equity
investment, long-term bound investment, fixed assets and immaterial assets.
The investment measures the change of firm's total investment compared to
the last period. All these four variables can be easily obtained from the
balance sheet. We directly use Tobin's q from the CSMAR database, where by
definition $q=\frac{MV}{K}$ and $MV$ is the market value of the firm. We obtain a balanced panel dataset with 195 firms and 72
time periods.
The units of corporate investment, capital and cash flow are measured by billions of
Chinese RMB.

By using the SVT approach, we obtain the estimates of the ranks of $\Theta
_{1}$ and $\Theta _{2}:$ $\hat{r}_{1}=\hat{r}_{2}=1$ for each $\tau
=\{0.25,0.5,0.75\}$. Consequently, we can consider the test that whether $
\Theta _{j,it}$ is constant over $i$ or constant over $t$ for both $j=1,2$.
Specifically, we want to test whether the effect of cash flow and Tobin's q
on the firm's investment is homogeneous over $i$ or across $t$ with market
imperfection. That is, for $j\in \{1,2\},$ we shall test

\begin{itemize}
\item $H_{0}^{a}$: $\Theta _{j,it}$ is a constant over $i$,

\item $H_{0}^{b}$: $\Theta _{j,it}$ is a constant over $t$.
\end{itemize}

Figure \ref{fig:emp_1} shows the estimation results for the factor and
factor loadings of two slope coefficient matrices under different quantiles.
In each sub-figure, the first and second rows report the results for $\Theta
_{1}$ and $\Theta _{2},$ respectively. Specifically, the first row of Figure
\ref{fig:emp_1}(a) gives the plot of $\left\{ \hat{u}_{i,1}\right\} _{i\in
\lbrack N]}$ as a catenation of $\{\hat{u}_{i,1}^{(1,2)},\hat{u}
_{i,1}^{(2,3)},\hat{u}_{i,1}^{(3,1)}\}$ at the left and as a cantenation of $
\{\hat{u}_{i,1}^{(1,3)},\hat{u}_{i,1}^{(2,1)},\hat{u}_{i,1}^{(3,2)}\}$ at
the right in the first row, and similarly the plot of $\left\{ \hat{u}
_{i,2}\right\} _{i\in \lbrack N]}\ $in the second row. Similarly, the first
row of Figure \ref{fig:emp_1}(d) shows $\{\hat{v}_{t,1}^{(a,b)}\}_{t\in
\lbrack T]}$ for $a\in \lbrack 3],$ $b\in \lbrack 3]\setminus \{a\}$ in the
first row and $\{\hat{v}_{t,3}^{(a,b)}\}_{t\in \lbrack T]}$ for $a\in
\lbrack 3],$ $b\in \lbrack 3]\setminus \{a\}$ in the second row.

Table \ref{tab:emp_1} reports the test statistics, critical values, and $p$-values. Tobin's q can measure a firm's investment demand. After controlling
the Tobin's q and the intercept FEs, the coefficient of cash flow captures a
firm's potential for external investment with the variation of internal
finance. It is clear that we can reject the homogeneous hypotheses for both $
i$ and $t$ at the 1\% significance level for each $\tau \in
\{0.25,0.5,0.75\} $. This indicates that with high probability, the slope
coefficient of both $CFK$ and Tobin's q follow the factor structure with one
factor.

The above study shows strong evidence that under imperfect market, the
sensitivity of corporate investment to cash flow exhibits both individual
heterogeneity and time heterogeneity across quantiles. It implies that
neither the usual homogenous panel QR model nor the panel QR model with
either cross-section or time heterogeneity alone in the slope coefficients
fails to fully capture the unobserved heterogeneity in the investment
equation.


\begin{table}[tbph]
\caption{Test results under different quantiles for the investment equation}
\label{tab:emp_1}
\scriptsize \centering
\begin{tabular}{ccccccc}
\toprule \toprule $\tau$ & Test & $S$ & $cv_{\alpha=0.01}$ & $
cv_{\alpha=0.05}$ & $cv_{\alpha=0.1}$ & $p$-value \\
\midrule \multirow{4}[2]{*}{0.25} & $u_{CFK}$ & $1.28\times10^{3}$ &
\multirow{2}[1]{*}{16.94} & \multirow{2}[1]{*}{13.68} &
\multirow{2}[1]{*}{12.24} & 0.00 \\
& $u_{q}$ & $4.16\times10^{4}$ &  &  &  & 0.00 \\
& $v_{CFK}$ & 13.85 & \multirow{2}[1]{*}{5.70} & \multirow{2}[1]{*}{4.07} &
\multirow{2}[1]{*}{3.35} & 0.00 \\
& $v_{q}$ & 870.85 &  &  &  & 0.00 \\
\midrule \multirow{4}[2]{*}{0.50} & $u_{CFK}$ & 148.28 &
\multirow{2}[1]{*}{16.94} & \multirow{2}[1]{*}{13.68} &
\multirow{2}[1]{*}{12.24} & 0.00 \\
& $u_{q}$ & $1.24\times10^{5}$ &  &  &  & 0.00 \\
& $v_{CFK}$ & 49.57 & \multirow{2}[1]{*}{5.70} & \multirow{2}[1]{*}{4.07} &
\multirow{2}[1]{*}{3.35} & 0.00 \\
& $v_{q}$ & 138.83 &  &  &  & 0.00 \\
\midrule \multirow{4}[2]{*}{0.75} & $u_{CFK}$ & 313.21 &
\multirow{2}[1]{*}{16.94} & \multirow{2}[1]{*}{13.68} &
\multirow{2}[1]{*}{12.24} & 0.00 \\
& $u_{q}$ & $2.03\times10^{4}$ &  &  &  & 0.00 \\
& $v_{CFK}$ & 31.50 & \multirow{2}[1]{*}{5.70} & \multirow{2}[1]{*}{4.07} &
\multirow{2}[1]{*}{3.35} & 0.00 \\
& $v_{q}$ & 58.29 &  &  &  & 0.00 \\
\bottomrule &  &  &  &  &  & \\
\multicolumn{7}{p{10cm}}{$Notes$: $S$ is the test statistics for the factor
or factor loadings under different quantiles, $H_{0}^{a}(CFK)$ and $
H_{0}^{a}(q) $ refer to the hypotheses that the slope of CFK and Tobin'q is
homogeneous across $i$, respectively. $H_{0}^{b}(CFK)$ and $H_{0}^{b}(q)$
refer to the the hypotheses that the slope of CFK and Tobin'q is homogeneous
across $t$, respectively. $\text{cv}_{\alpha=a}$ is the critical value under
the significance level a where a=0.1, 0.05, and 0.01.}
\end{tabular}
\end{table}

\begin{figure}[tbp]
\centering
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.25_u.png}
         \caption{factor loading estimates under $\tau=0.25$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.5_u.png}
         \caption{factor loading estimates under $\tau=0.5$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.75_u.png}
         \caption{factor loading estimates under $\tau=0.75$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.25_v.png}
         \caption{factor estimates under $\tau=0.25$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.5_v.png}
         \caption{factor estimates under $\tau=0.5$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp1_tau0.75_v.png}
         \caption{factor estimates under $\tau=0.75$}
     \end{subfigure}
\caption{Factor loading and factor estimates under different quantiles}
\label{fig:emp_1}
\end{figure}

\subsection{Foreign Direct Investment and Unemployment}

Investment is one of the major driving forces for economic growth and
employment. Among the investment, foreign direct investment (FDI) is an
important contributor to the employment. See \cite{craigwell2006foreign},
\cite{aktar2009can}, \cite{karlsson2009foreign}, \cite{mucuk2013effect}, and
\cite{strat2015fdi}, among others. Controversially, \cite{mucuk2013effect}
argue that FDI may have both positive and negative effects on employment. On
the one hand, FDI adds to the net capital and creates jobs through forward and backward linkages and
multiplier effects in local economy. On the other hand, acquisitions may
rely on imports or displacement of existing firms which may result in job
loss.

To study the relationship of FDI, economic growth rate and unemployment at
the country level, we consider the following panel quantile regression
model,
\begin{equation*}
\mathscr{Q}_{\tau }\left( U_{it}\big |\left\{ G_{i,t-1},FDI_{it}\right\}
_{t\in \lbrack T]},\left\{ \Theta _{j,it}\right\} _{t\in \lbrack T],j\in
\{0,1,2\}}\right) =\Theta _{0,it}(\tau )+\Theta _{1,it}(\tau
)G_{i,t-1}+\Theta _{2,it}(\tau )FDI_{it},
\end{equation*}
where $U_{it}$ is the unemployment rate of country $i$ at year $t$, $
G_{i,t-1}$ is the economic growth measured by the growth of real GDP. $
\Theta _{0,it}$ is the FEs of country $i$ and year $t$, $\Theta _{1,it}$ is
the elasticity of the economic growth in the previous year to the
unemployment this year, and $\Theta _{2,it}$ is the elasticity of FDI to the
unemployment.

We draw the data for 126 countries from 1992-2019. The data for the
unemployment rate are taken from International Labor Organization (ILO) and
GDP growth and FDI are from the World Bank Development Indicators (WDI)
historical database. The rank estimation procedure shows that $\hat{r}_{1}=2$
and $\hat{r}_{2}=1$. Consequently, we can test whether the elasticity of FDI
to the unemployment rate is homogeneous across individual countries and over
years 1992-2009, and whether the elasticity of growth rate to unemployment
follows the additive structure, i.e.,

\begin{itemize}
\item $H_{0}^{c}$: $\Theta_{1,it}=\Theta_{1,i}+\Theta_{1,t}$,

\item $H_{0}^{d}$: $\Theta_{2,it}$ is a constant over $i$,

\item $H_{0}^{e}$: $\Theta_{2,it}$ is a constant over $t$.
\end{itemize}

Table \ref{tab:emp_2} reports the test results under quantiles 0.25, 0.5 and
0.75 for the above three null hypotheses. Figure \ref{fig:emp_2} gives the
estimation results for the factor and factor loading estimates of the slope
coefficient $\Theta _{2}$. As Table \ref{tab:emp_2} suggests, we can reject
all the above three null hypotheses safely at the conventional 5\%
significance level. This means that the effect of FDI on the unemployment
rate is different across both countries and time even though the estimated
rank of $\Theta _{2}$ is one, and the effect of economic growth rate on the
unemployment is heterogeneous across both countries and time and it does not
exhibit an additive structure.





\begin{table}[htbp]
  \caption{Test results under different quantiles}
  \label{tab:emp_2}
   \scriptsize \centering
    \begin{tabular}{ccccccc}
    \toprule
    \toprule
    Test  & $\tau$ & $S$   & $cv_{\alpha=0.01}$ & $cv_{\alpha=0.05}$ & $cv_{\alpha=0.10}$ & $p-$value \\
    \midrule
    \multirow{3}[2]{*}{$H_{0}^{c}$} & 0.25  & 38.92 & 35.55 & 32.29 & 30.85 & 0.00 \\
          & 0.50   & 80.84 & 22.29 & 19.03 & 17.59 & 0.00 \\
          & 0.75  & 66.24 & 35.55 & 32.29 & 30.85 & 0.00 \\
    \midrule
    \multirow{3}[2]{*}{$H_{0}^{d}$} & 0.25  & $1.41\times 10^{6}$ & \multirow{3}[2]{*}{16.15} & \multirow{3}[2]{*}{12.89} & \multirow{3}[2]{*}{11.45} & 0.00 \\
          & 0.50   & $6.39\times 10^{6}$ &       &       &       & 0.00 \\
          & 0.75  & $3.07\times 10^{7}$ &       &       &       & 0.00 \\
    \midrule
    \multirow{3}[2]{*}{$H_{0}^{e}$} & 0.25  & 36.36 & \multirow{3}[2]{*}{5.70} & \multirow{3}[2]{*}{4.07} & \multirow{3}[2]{*}{3.35} & 0.00 \\
          & 0.50   & 164.03 &       &       &       & 0.00 \\
          & 0.75  & 5.44  &       &       &       & 0.013 \\
    \bottomrule
    \end{tabular}
\end{table}



\begin{figure}[h]
\centering
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.25_u.png}
         \caption{factor loading estimates under $\tau=0.25$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.5_u.png}
         \caption{factor loading estimates under $\tau=0.5$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.75_u.png}
         \caption{factor loading estimates under $\tau=0.75$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.25_v.png}
         \caption{factor estimates under $\tau=0.25$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.5_v.png}
         \caption{factor estimates under $\tau=0.5$}
     \end{subfigure}
\hfill
\begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth,height=2.5cm]{emp2_tau0.75_v.png}
         \caption{factor estimates under $\tau=0.75$}
     \end{subfigure}
\caption{Factor loading and factor estimates of $\Theta_{2}$ under different
quantiles}
\label{fig:emp_2}
\end{figure}

\section{Conclusion}

This paper considers panel QR model with heterogeneous slopes over both $i$ and $t$. Compared to \cite{chernozhukov2019inference}, to remove the bias from the nuclear norm regularization, we split the full sample into three subsamples. We then use the first subsample to compute initial estimators via NNR, the second sample to refine the convergence rate of the initial estimator, and the last subsample to debias the refined estimator. Our asymptotic theory shows that the factor estimates, factor loading estimates and the slope estimates all follow the normal distributions asymptotically. By constructing the consistent estimator for the asymptotic variance, we also conduct two specification tests: (1) the slope coefficient is constant over time or individuals under the case that true rank of slope matrix equals one and (2) the slope coefficient exhibits the additive structure under the case that true rank of the slope coefficient matrix equals two. Our test statistics are shown to follow the Gumbel distribution asymptotically under the null, consistent under the global alternative and have non-trivial power against local alternatives.  Monte Carlo simulation and empirical studies illustrate the finite sample performance of our algorithm and test statistics.


\bibliography{chapter1}
\newpage