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.
87,585 characters
Panel Stochastic Frontier Models with Latent Group Structures
\def\spacingset#1{\renewcommand{1.2}
{#1}\small\normalsize}
\spacingset{1}
\title{Panel Stochastic Frontier Models with Latent Group Structures\thanks{We thank the editor, Michal Koles\'ar, an Associate Editor, and two
anonymous referees for their helpful comments. We are grateful for
Bin Peng, Valentin Zelenyuk, and the participants of the Econometric Society
Australasian Meeting 2024 and AE$^2$ conference 2025 for their helpful feedback. The first author acknowledges that this research was supported by the Commonwealth through an Australian Government Research Training Program Scholarship \href{https://doi.org/10.82133/C42F-K220}{[DOI: https://doi.org/10.82133/C42F-K220]}, as well as the Japan-IMF Scholarship. The second author gratefully acknowledges the financial support of the
National Natural Science Foundation of China (No. 72394392 and No. 72573158).}}
\author{{Kazuki Tomioka}\thanks{Graduate School of Humanities and Social Sciences, Hiroshima University, Hiroshima, Japan. Email: [email removed]}.}\\
Hiroshima University \and {Thomas T. Yang}\thanks{Corresponding author. Research School of Economics, The Australian
National University, Canberra, ACT 2601, Australia. Email: [email removed]}.}\\
Australian National University \and { Xibin Zhang}\thanks{Department of Econometrics and Business Statistics, Monash University,
Caulfield East, Victoria 3145, Australia. Email: [email removed]}.}\\
Monash University}
\maketitle
\vspace{0.5cm}
\begin{abstract}
\noindent Stochastic frontier models have attracted considerable attention due to the incorporation of an inefficiency term in addition to the conventional error term. In this paper, we propose a general estimation framework for panel stochastic frontier models that accommodates potential heterogeneity through latent group structures. The framework is tailored to the distinctive features of stochastic frontier models and is paired with a practical hybrid estimation procedure that combines individual-level and joint panel estimation. We illustrate the estimation framework using a panel stochastic frontier model that treats the inefficiency term as a random effect, and show that it can be readily extended to a range of fixed effects specifications common in the literature. Simulation studies indicate strong finite-sample performance, and we further demonstrate the practicality of the approach in an empirical application to the cost efficiency of the U.S. commercial banking sector.
\end{abstract}
\begin{description}
\item [{Keywords:}] Classification, Group Structures, Panel Data, Stochastic
Frontier
\item [{\emph{JEL classification:}}] C23, C33, C38, C51
\end{description}
\pagebreak{}
\newpage\spacingset{1.8}
\section{Introduction}
\citet{Aigneretal1977} and \citet{MeeusenvanDenBroeck1977} introduced
the stochastic frontier (SF) model to study the productive (in)efficiency
of firms, and such SF models have since attracted considerable attention.\footnote{We refer readers to \citet{KumbhakarLovell2000} for early developments and \citet{Kumbhakaretal2022a} and \citet{Tsionas2023} for recent advancements and a comprehensive review of the literature.} One distinct feature of a SF model is the decomposition of the error term, typically expressed as $\varepsilon=v-u$, where $v$ is
a random disturbance and $u\geq0$ is the inefficiency term.
Unlike standard regression models, the identification of $v$ and
$u$ is crucial in SF modeling because of their economic interpretations.
In this paper, we develop a general estimation framework for time-varying panel SF models that can accommodate potential heterogeneity across firms through latent group structures. To illustrate our methodology, we consider a model similar to \citet{YaoZhangKum2019}'s that allows for heterogeneous time-varying coefficients:
\begin{align}
\label{EQ:intromodel}
y_{it} &=\alpha_{i}^{0}+\alpha_{i}(\tau_{t})+x_{it}^{\prime}\beta_{i}(\tau_{t})+\varepsilon_{it},\quad\text{with}\quad\varepsilon_{it}=v_{it}-u_{i}, \nonumber\\
&= \alpha_{i}^{0}-u_{i}+\alpha_{i}(\tau_{t})+x_{it}^{\prime}\beta_{i}(\tau_{t})+v_{it},
\end{align}
for firm $i=1,2,\dots,N$ and time $t=1,2,\dots,T$, where $\tau_{t}=t/T\in(0,1]$. In this specification of the model, $\alpha_{i}^{0}$ and $\alpha_{i}(\tau_{t})$ denote the constant and time-varying intercepts of the frontier, respectively. The term $\beta_{i}(\tau_{t})$ denotes a non-random, time-varying coefficient vector, $v_{it}$ is a zero-mean random error term, and $u_{i}\geq0$ represents a firm-specific random inefficiency term. In addition to allowing for heterogeneity in the efficient frontier, we permit the variance of $v_{it}$ to vary across firms, and allow for the possibility that $u_{i}$ follows a mixture distribution.
Our study is motivated by the role of heterogeneity in the measurement of the efficient frontier and the inefficiency in panel SF models. Heterogeneity can either shift the efficient frontier or distort the location and scale of inefficiency estimates \citep[see, for example,][]{Galan2014}. According to \citet{Greene2005a}, the true underlying frontier may include unmeasured firm-specific characteristics that reflect the technology in use. \citet{Greene2005a} was instrumental in expanding SF models to incorporate firm heterogeneity by allowing for either true fixed effects or true random effects. Building on this, heterogeneity was further explored in subsequent work by allowing the inefficiency term to be purely transient or to contain both transient and persistent components \citep[see, for example][]{Colombi2014,Kumbhakar2014,Tsionas2014}.
These approaches help disentangle unobserved heterogeneity from inefficiency. However, a homogeneous frontier implicitly assumes that all firms operate under the same technology, so any systematic deviation is attributed to inefficiency. As emphasized in \citep{Greene2005a, Greene2005}, failing to distinguish the frontiers and inefficiency leads to a structural misinterpretation: persistent heterogeneity is absorbed into the inefficiency term, resulting in biased efficiency estimates.
The framework we propose addresses this issue by introducing a latent group structure for frontier parameters. Firms are partitioned into a small number of groups, with each group sharing a common frontier that reflects a particular technological regime (e.g., distinct business models or scale-specific production technologies). The finite collection of latent frontiers provides a middle ground between a fully homogeneous frontier and unrestricted firm-specific frontiers.
Formally, we assume that firms can be classified into one of $K^{*}\geq1$ groups. Within each group, firms share a common set of parameters $\left\{ \alpha(\tau_{t}),\beta(\tau_{t}),\text{Var}(v_{t})\right\}$. We show that the classification step is consistent, ensuring that firms are benchmarked against the correct regime with probability approaching one.
Importantly, we do not impose the same group structure for the distribution of the inefficiency term, $u$ or on the constant term, $\alpha^{0}$. We find that defining group membership for $\{\alpha^{0},u\}$ in the
same way as for the other parameters is inappropriate, the reasons of which we discuss in detail
in Section~\ref{SEC:inefficiency}. Instead, we account for potential heterogeneity
in $\alpha^{0}-u$ by modeling it as a mixture of distributions.
The idea of uncovering latent group structures in panel data
models has been extensively studied in recent years. Existing methods can be broadly categorized into two main approaches. The first approach relies on a two-step procedure in which individual-level estimation is followed by applying a clustering algorithm such as K-means \citep{LinNg2012, BonhommeManresa, AndoBai2016}, or Hierarchical Agglomerative Clustering (HAC) \citep{Chen2019}.
The second approach, introduced by
\citet{SuEtal2016} performs simultaneous estimation and classification by penalizing the joint mean squared errors. Further developments along this line include \citet{SuEtal2019}, \citet{HuangEtal2020} and \citet{WangSu2021}. Conceptually, the latter approach pools all observations
across firms for joint estimation and classification, whereas the
former relies on individual (non-pooled) firm-specific estimates as the basis of classification.
These existing approaches cannot be applied directly in our setting for the following reason. Estimating the parameters that govern the underlying distribution of $u$ requires pooling observations across firms, typically via maximum likelihood based on the likelihood function in (\ref{EQ:likelihood_appro}). However, since the log-likelihood function is complex and highly nonlinear, applying the method of \citet{SuEtal2016} that requires simultaneous estimation and classification for SF models is nontrivial. This leads to a dilemma of whether to pool or not to pool observations for inference. To resolve this, we propose a new hybrid approach that combines pooled and non-pooled estimation methods that is flexible enough to be applied to a broad class of SF models.
Our paper makes several contributions to the literature on SF models
and classification methods for panel data models. First, we estimate
a robust panel SF model that incorporates heterogeneity across firms.
We carefully set out the framework, clarifying what is feasible and what is not, and provide detailed explanations of the underlying rationale.\footnote{Previous work on robust estimation of the SF model has primarily focused on semi-parametric or non-parametric specifications for either the efficient frontier \citep{ParkSimar1994, YaoZhangKum2019} or the error term distribution \citep{Greene2005, LaiKumbhakar2023}. Additionally, there is a growing literature on specification tests for the distribution of inefficiency \citep{Chengetal2024}.} Second, to the best of our knowledge, this is one of the few papers in the literature to model and allow the variance of the error term to exhibit latent group structures. A notable
study in this area is \citet{LoyoBoot}, where they focus on modeling the variance of the error term, primarily for efficiency gains and/or the specific role played by the variance. In contrast, we aim to uncover heterogeneity in both the error and inefficiency terms, for the economic interpretations given in SF models. We emphasize the importance of modeling heterogeneity in both the variance of the error term, $v$, and the distribution of inefficiency, $u$, since inefficiency estimates depend on the variance of $v$ in the widely used \citet{Jondrowetal1982}'s point estimator. Given the importance of the joint modeling of both the variance of $v$ and the distribution of $u$ in SF models, we consider this to be a substantial contribution. Moreover, we propose a hybrid estimation procedure that combines firm-level and joint panel estimation to address the potential heterogeneity present in both the frontier and error components.
As the description of the SF model in (\ref{EQ:intromodel}) alluded, we present the hybrid approach in the main text using a random effects specification of the SF model, and show in Appendix \ref{APP:FE} how it extends to other variants, including the fixed effects SF models studied by \citet{Greene2005a, Greene2005, Chenetal2014}, and \citet{ZhouEtal2020}. Finally, the separation between the frontier and inefficiency may be complicated by model misspecification. We rule out this possibility by assuming correct model specification. Developing procedures that are more robust to misspecification is an important direction for future research.
The rest of this paper is organized as follows. In Section~\ref{SEC:model}, we introduce the estimation procedure and propose information criteria to determine the number of groups and whether $u$ follows a unique or a mixture of distributions. Section~\ref{SEC:Theoretical} examines the theoretical properties of our procedure, specifying the conditions required for tuning parameters. In Section~\ref{SEC:simulations}, we evaluate the small-sample performance of our method, providing practical recommendations for tuning parameters that meet the conditions and perform well in simulations. In Section \ref{SEC:application}, we apply our procedure to a dataset of large U.S. commercial banks, using the recommended tuning parameters from the simulations. Our findings confirm the presence of heterogeneity in the frontiers and a mixture distribution for the inefficiency term. The Online Appendix includes several important supplementary results. Appendix \ref{APP:FE1} extends our approach to fixed effects SF models under distributional assumptions. Appendix \ref{APP:FE2} further generalizes this to fixed effects SF models without any distributional restrictions, allowing for more flexible conditions on the inefficiency term. Appendix \ref{APP:innocuous} discusses a relaxation of the normalization condition. Appendix \ref{APP:ZISF} addresses the possibility that the inefficiency term equals zero with positive probability. Appendix \ref{APP:more_mix} considers the mixture distribution with more than two components. Other parts of the Appendix support the results in the main body of the paper. Specifically, Appendices \ref{APP:likeli} and \ref{APP:HAC} describe, respectively, the approximate likelihood function of the model and the clustering method used. Proofs of the theorems and propositions from Section \ref{SEC:Theoretical} are provided in Appendix~\ref{APP:mainproofs}, with technical lemmas presented in Appendix \ref{APP:lemmaProofs}. Additional simulation and application results are included in Appendix~ \ref{APP:Tables}, and Appendix~\ref{App:derivatives} further discusses the properties of the information matrix in support of Appendix~\ref{APP:mainproofs}.
\section{Model and Estimation}
\label{SEC:model}
This section presents the model and the procedure. To facilitate exposition,
we assume $\alpha_{i}^{0}=\alpha^{0}$ and $\text{Var}(u_{i}) = \sigma_{ui}^{2} = \sigma_{u}^{2}$ for all $i$ in Sections \ref{sec:The-Model-More} to \ref{SEC:Estimation}. These restrictions are relaxed to the general case subsequently. The technical conditions and main theorems are postponed to the next section.
\subsection{The Model}\label{sec:The-Model-More}
We illustrate our approach using the panel SF model with random effects (RE) and relegate other variants of the model to Appendix \ref{APP:FE}. The particular model we consider modifies \citet{YaoZhangKum2019}'s by incorporating heterogeneous, time-varying coefficients:
\begin{align}
y_{it} &=\alpha^{0} + \alpha_{i}\left(\tau_{t}\right)+x_{it}^{\prime}\beta_{i}\left(\tau_{t}\right)+\varepsilon_{it}\nonumber \\
&=\alpha^{0}-u_{i}+\alpha_{i}\left(\tau_{t}\right)+\sum_{l=1}^{p}x_{itl}\beta_{il}\left(\tau_{t}\right)+v_{it},\label{EQ:model}
\end{align}
noting the temporal restriction $\alpha_{i}^{0}=\alpha^{0}$,
for expositional purposes to be generalized subsequently.
The subscript $l$ denotes the $l$-th element of a vector $x_{it}$, which
is $p$ by $1$. In this specification, $\varepsilon_{it}=v_{it}-u_{i}$,
where $v_{it}$ is a mean-zero random error term, and $u_{i}$ is
a non-negative term capturing the inefficiency of firm $i$.\footnote{We describe the model, estimation method, and simulations in terms
of the production frontier model, but use the cost frontier model
for our application. The only difference is that the inefficiency
term, $u_{i}$, enters the model negatively (production frontier)
or positively (cost frontier). This distinction is minor, and one
can let $\varepsilon_{it}\equiv v_{it}+u_{i}$ for cost frontiers.} We let $\alpha_{i}(\cdot)$ and
$\beta_{i}(\cdot)$ evolve smoothly in $\tau_{t}$ to capture the gradual technological drift, regulatory cycles, and business-model adjustments that are pervasive in long panels (for example, banking costs). Smoothness avoids implausible jumps while allowing flexible, low-frequency changes that pooled static frontiers cannot capture.
As in \citet{YaoZhangKum2019}, we assume that
\begin{equation}
\label{why.RE}
v_{it}\sim N(0,\sigma_{vi}^{2}), \quad u_{i}\sim|N(0,\sigma_{u}^{2})|, \quad\text{and} \quad v_{it}\perp u_{i}\perp x_{it}.
\end{equation}
We restrict $\sigma_{u}^{2}$ to be identical for all $i$ (no group structure) momentarily, again to facilitate exposition. We note that the two parameters where homogeneity is imposed, $\alpha^{0}$ and $\sigma_{u}^{2}$, exhibit distinctive features compared to other parameters. Given the importance of these parameters in SF models, we devote a separate section to discuss these differences in detail in Section \ref{SEC:inefficiency}.
The assumption given by (\ref{why.RE}) implies that inefficiency arises from managerial, organizational or behavioral factors that are unrelated to observed input or output variables in the frontier. The resulting panel SF model is RE in the sense of \citet{Greene2005a}. In panel SF models, RE is attractive because fixed effects (FE) type likelihoods face an incidental-parameters problem in finite $T$. RE avoids this by treating unit effects probabilistically, which \citet{Tsionas2014} emphasize when arguing for a fully likelihood-based/Bayesian route for panel SF model analysis with multiple error components. We adopt this setup primarily to illustrate our proposed approach, although the methodology can be extended to alternative model specifications, such as the four-component panel stochastic frontier model introduced by \citet{Tsionas2014} and \citet{LaiKumbhakar2023}. FE-type models are discussed separately in Appendices \ref{APP:FE1} and \ref{APP:FE2}, as their treatment is relatively straightforward given the procedure developed in the main body of the paper.
We assume that there are $K^{\ast}\geq1$ groups of parameters, and
each firm's parameters belong to one of these groups. Mathematically,
\begin{equation}
\left\{ \alpha_{i}\left(\tau_{t}\right),\beta_{i}\left(\tau_{t}\right),\sigma_{vi}\right\} =\sum_{k=1}^{K^{\ast}}\left\{ \alpha_{(k)}^{\ast}\left(\tau_{t}\right),\beta_{(k)}^{\ast}\left(\tau_{t}\right),\sigma_{v(k)}^{\ast}\right\} \boldsymbol{1}(i\in G_{k}),\label{EQ:group_para}
\end{equation}
where $\boldsymbol{1}(\cdot)$ is the indicator function, equaling
1 if $(\cdot)$ is true and 0 otherwise. Additionally, parameters from
different groups are distinct, meaning
$
\left\{ \alpha_{(k)}^{\ast}\left(\tau_{t}\right), \beta_{(k)}^{\ast}\left(\tau_{t}\right),\sigma_{v(k)}^{\ast}\} \neq\{ \alpha_{(j)}^{\ast}\left(\tau_{t}\right),\beta_{(j)}^{\ast}\left(\tau_{t}\right),\sigma_{v(j)}^{\ast}\right\} ,
$ for $j\neq k$. The group membership sets satisfy
$G_{j}\cap G_{k}=\emptyset\text{ and }\bigcup_{k=1}^{K^{\ast}}G_{k}=\left\{ 1,2,\dots,N\right\}$.
It is worth noting that we impose the following assumption on each $\alpha_{(k)}^{*}(s)$:
\begin{equation}
\int_{0}^{1}\alpha_{(1)}^{\ast}(s)\,\textrm{d}s = \int_{0}^{1}\alpha_{(2)}^{\ast}(s)\,\textrm{d}s =\ldots=\int_{0}^{1}\alpha_{(K^{\ast})}^{\ast}(s)\, \textrm{d}s,\label{eq:restriction_level}
\end{equation}
although we generalize it to allow them to differ in Appendix \ref{APP:innocuous}.\footnote{We explain in detail why we impose this condition, how we adjust the procedure without it, and when we recommend relaxing the condition in Appendix \ref{APP:innocuous}.} The normalization adopted here is $\int_{0}^{1}\alpha_{(k)}^{\ast}(s)\,\textrm{d}s=0$.
Clearly, when $K^{\ast}=1$ (the homogeneous case), this normalization
is innocuous; however, it is not in the general case due to the restriction in (\ref{eq:restriction_level}). When the intercept term does not vary over time, $\alpha_{i}(s)=0$. This normalization ensures that $\alpha_{i}(s)$ captures the time-varying component of the intercept term.
\subsection{Approximation of \texorpdfstring{$\alpha(\cdot)$}{a} and \texorpdfstring{$\beta(\cdot)$}{b} }\label{SEC:serial_approximation}
The approximations we adopt are standard in the literature. Let $L^{2}\left[0,1\right]=\{f\left(s\right):\int_{0}^{1}f^{2}\left(s\right)\text{d}s<\infty\}$
represent the space of square-integrable functions. The inner product equipped on this space is defined as $\left\langle f_{1},f_{2}\right\rangle \equiv\int_{0}^{1}f_{1}\left(s\right)f_{2}\left(s\right)\text{d}s$,
and the induced norm is $\left\Vert f\right\Vert =\left\langle f,f\right\rangle ^{1/2}$. Following \citet{DongLinton2018} and \citet{Ataketal}, we use cosine functions as basis functions. In particular, $B_{0}\left(s\right)=1$ and $B_{j}\left(s\right)=\sqrt{2}\cos(j\pi s)$
for $j\geq1$. The set $\{B_{j}\left(s\right)\}_{j=0}^{\infty}$ then forms
an orthonormal basis for the Hilbert space $L^{2}\left[0,1\right]$,
such that $\left\langle B_{i},B_{j}\right\rangle =\delta_{ij}$, where
$\delta_{ij}$ is the Kronecker delta.
Suppose $f\in L^{2}\left[0,1\right]$ is $\kappa$-th order continuously differentiable. Then, we have
\begin{align*}
f\left(s\right) &= \sum_{j=0}^{\infty}B_{j}\left(s\right)v_{j}^{0}=\sum_{j=0}^{m-1}B_{j}\left(s\right)v_{j}^{0}+\sum_{j=m}^{\infty}B_{j}\left(s\right)v_{j}^{0}\\
&= \sum_{j=0}^{m-1}B_{j}\left(s\right)v_{j}^{0}+O\left(m^{-\kappa}\right)\equiv\mathbb{B}^{m}\left(s\right)^{\prime}v^{0}+O\left(m^{-\kappa}\right),
\end{align*}
where $\mathbb{B}^{m}\left(s\right)\equiv\left(B_{0}\left(s\right),B_{1}\left(s\right),\dots,B_{m-1}\left(s\right)\right)^{\prime}$,
$v_{j}^{0}=\left\langle f,B_{j}\right\rangle $, and $v^{0}=\left(v_{0}^{0},v_{1}^{0},\dots,v_{m-1}^{0}\right)^{\prime}$.
Here, $\sum_{j=m}^{\infty}B_{j}\left(s\right)v_{j}^{0}$ is the bias
term from using only the first $m-1$ terms contained in $\mathbb{B}^{m}\left(s\right)$ to approximate $f\left(s\right)$. If $f\left(s\right)$ is $\kappa$-th order differentiable,
the bias term is $\sum_{j=m}^{\infty}B_{j}\left(s\right)v_{j}^{0} = O\left(m^{-\kappa}\right)$. When $\int_{0}^{1}f\left(s\right)\text{d}s=0$ is imposed, we approximate $f$ using $\mathbb{B}_{-0}^{m}\left(s\right)\equiv\left(B_{1}\left(s\right),\dots,B_{m-1}\left(s\right)\right)^{\prime}$,
since $B_{0}(s)=1$ and $\int_{0}^{1}B_{j}\left(s\right)\text{d}s=0$ for $j\geq1$. Thanks to this property, it is more convenient than using alternative bases such as B-splines. Similarly,
\[
f\left(s\right)=\mathbb{B}_{-0}^{m}\left(s\right)v_{-0}^{0}+O\left(m^{-\kappa}\right),
\]
for some $v_{-0}^{0}=\left(v_{1}^{0},v_{2}^{0},\dots,v_{m-1}^{0}\right)^{\prime}$.
We apply this approximation to our case. For each firm $i$, we have
\begin{align}
y_{it} & =\alpha^{0}-u_{i}+\alpha_{i}\left(\tau_{t}\right)+\sum_{l=1}^{p}x_{itl}\beta_{il}\left(\tau_{t}\right)+v_{it}\nonumber \\
& \approx\alpha^{0}-u_{i}+\mathbb{B}_{-0}^{m}\left(\tau_{t}\right)^{\prime}\pi_{i0}^{0}+\sum_{l=1}^{p}x_{itl}\mathbb{B}^{m}\left(\tau_{t}\right)^{\prime}\pi_{il}^{0}+v_{it}\nonumber \\
& \equiv\alpha^{0}-u_{i}+\left[\mathbb{B}_{-0}^{m}\left(\tau_{t}\right)^{\prime},\left(x_{it}\otimes\mathbb{B}^{m}\left(\tau_{t}\right)\right)^{\prime}\right]\pi_{i}^{0}+v_{it}\nonumber \\
& \equiv\alpha^{0}-u_{i}+z_{it}^{\prime}\pi_{i}^{0}+v_{it}\equiv\tilde{z}_{it}^{\prime}\tilde{\pi}_{i}^{0}+v_{it},\label{EQ:approx}
\end{align}
where the terms $\mathbb{B}_{-0}^{m}\left(\tau_{t}\right)^{\prime}\pi_{i0}^{0}$
and $\mathbb{B}^{m}\left(\tau_{t}\right)^{\prime}\pi_{il}^{0}$ represent
approximations of $\alpha_{i}\left(\tau_{t}\right)$ and $\beta_{il}\left(\tau_{t}\right)$,
respectively, and
\begin{align*}
x_{it}\otimes\mathbb{B}^{m}\left(\tau_{t}\right) & \equiv\left(x_{it1}B_{0}\left(\tau_{t}\right),\dots,x_{it1}B_{m-1}\left(\tau_{t}\right),\dots,x_{itp}B_{0}\left(\tau_{t}\right),\dots,x_{itp}B_{m-1}\left(\tau_{t}\right)\right)^{\prime},\\
\pi_{i}^{0} & \equiv\left(\pi_{i0}^{0\prime},\pi_{i1}^{0\prime},\dots,\pi_{ip}^{0\prime}\right)^{\prime},\quad\tilde{\pi}_{i}^{0}\equiv\left(\alpha^{0}-u_{i},\pi_{i0}^{0\prime},\pi_{i1}^{0\prime},\dots,\pi_{ip}^{0\prime}\right)^{\prime},\\
z_{it} & \equiv\left[\mathbb{B}_{-0}^{m}\left(\tau_{t}\right)^{\prime},\left(x_{it}\otimes\mathbb{B}^{m}\left(\tau_{t}\right)\right)^{\prime}\right]',\quad\text{and}\quad\tilde{z}_{it}\equiv\left(1,z_{it}^{\prime}\right)'.
\end{align*}
The last two lines of equation (\ref{EQ:approx}) represent three
equivalent ways of expressing the approximation.
\subsection{The Estimation}\label{SEC:Estimation}
The variance of inefficiency, $\sigma_{u}^{2}$ is identified
through the skewness in the distribution of $\alpha^{0}-u_{i}$ across
$i$. As a result, $\sigma_{u}^{2}$ cannot be identified or estimated
without pooling observations across different $i$. However, pooling observations for estimation introduces challenges for numerical optimization, since the log-likelihood functions in (\ref{EQ:likelihood_appro}) and (\ref{EQ:log_approx_p}) are complex and highly nonlinear. This issue exacerbates as the number of unknown parameters increases and becomes particularly pronounced in pooled
estimation with classification methods. This creates a dilemma regarding
whether to pool observations for estimation.
To address these challenges, we propose a hybrid procedure that combines
estimations with and without pooling. We present the detailed steps of the procedure below.
\subsubsection*{Step 1: Individual Estimation}
Using the approximation from (\ref{EQ:approx}), we regress $y_{it}$
against $\mathbb{B}^{m}\left(\tau_{t}\right)$ and $x_{it}\otimes\mathbb{B}^{m}\left(\tau_{t}\right)$
for $t=1,2,\ldots,T$ to obtain $\widehat{\tilde{\pi}}_{i}$. From
this we obtain $\hat{\sigma}_{vi}^{2}$ as the sample variance of
the regression residuals.
Specifically, consider the following expression: $Z_{im}\equiv (\tilde{z}_{i1},...,\tilde{z}_{iT})'$, a $T\times m(p+1)$ vector. Then the OLS estimator for firm $i$ is given by
\begin{equation}
\widehat{\tilde{\pi}}_{i}=\left(Z_{im}^{\prime}Z_{im}\right)^{-1}Z_{im}^{\prime}y_{i},\label{eq:pihat}
\end{equation}
with $y_{i}=\left(y_{i1},\ldots,y_{iT}\right)^{\prime}$. We then
obtain an estimate of $\sigma_{vi}^{2}$ as $\hat{\sigma}_{vi}^{2}=\frac{1}{T-1}\sum_{t=1}^{T}\left(y_{it}-\tilde{z}_{it}'\widehat{\tilde{\pi}}_{i}\right)^{2}.$ Excluding the first element in $\widehat{\tilde{\pi}}_{i}$, we let
$\hat{\pi}_{i}$ denote the estimated coefficients associated with
$z_{it}$. The estimates $\hat{\pi}_{i}$ and $\hat{\sigma}_{vi}$
are collected to form an estimate of $\vartheta_{i}$:
$\hat{\vartheta}_{i}=\left(\hat{\pi}_{i}^{\prime},\hat{\sigma}_{vi}\right)^{\prime},$ based on which we form groups.
\subsubsection*{Step 2: Classification}
Having obtained $\hat{\vartheta}_{1},\hat{\vartheta}_{2},\ldots,\hat{\vartheta}_{N}$
from Step 1, we use the $L_{2}$ norm to measure the distance between
$\hat{\vartheta}_{i}$ and $\hat{\vartheta}_{j}$. Based on this distance
measure, we then apply the classical HAC algorithm to the estimates of each firm's functional coefficient to determine group memberships. The HAC is a widely used algorithm for clustering, and several variants of it are employed in heterogeneous panel data models (see, for e.g., \citet{Chen2019}). Details of the HAC method are provided in Appendix \ref{APP:HAC} and refer the readers to \citet{Everittetal} for a comprehensive treatment. Given a value
for $K$, we apply the HAC to obtain an estimate of the group membership, denoted as $\left(\hat{G}_{1|K},\hat{G}_{2|K},\ldots,\hat{G}_{K|K}\right),$ which forms a partition of the set $\left\{ 1,2,\ldots,N\right\} $.
\subsubsection*{Step 3: Post-Classification Estimation and Determination of $K^{\ast}$}
Within each group, we now have significantly more observations available
for pooling. Recognizing this, we set the number of sieve terms to
$\underline{m}$, which is substantially larger than $m$. Within each estimated group, $\hat{G}_{k|K}$ for $1\leq k\leq K$, we conduct post-classification estimation using standard within-panel data estimation methods. Let $\underline{z}_{it}=\left[\mathbb{B}_{-0}^{\underline{m}}\left(\tau_{t}\right)^{\prime},\left(x_{it}\otimes\mathbb{B}^{\underline{m}}\left(\tau_{t}\right)\right)^{\prime}\right]^{\prime}$ denote the new regressors. At this stage, we do not consider the inefficiency term $\alpha^{0}-u_{i}$. The group specific coefficient is given by
\[
\hat{\pi}_{(k|K)}=\arg\min_{\pi}\sum_{i\in\hat{G}_{k|K}}\sum_{t=1}^{T}\left(\ddot{y}_{it}-\underline{\ddot{z}}_{it}^{\prime}\pi\right)^{2},
\]
where $\ddot{y}_{it}=y_{it}-\frac{1}{T}\sum_{t=1}^{T}y_{it}$, and $\underline{\ddot{z}}_{it}=\underline{z}_{it}-\frac{1}{T}\sum_{t=1}^{T}\underline{z}_{it}.$ The estimate of the variance of $v_{it}$ for group $\hat{G}_{k|K}$ is
\[
\hat{\sigma}_{v(k|K)}^{2}=\frac{1}{N_{k}(T-1)}\sum_{i\in\hat{G}_{k|K}}\sum_{t=1}^{T}\left(\ddot{y}_{it}-\underline{\ddot{z}}_{it}^{\prime}\hat{\pi}_{(k|K)}\right)^{2},
\]
where $N_{k}=\sharp\{\hat{G}_{k|K}\}$ is the number of elements in
$\hat{G}_{k|K}$. For simplicity, we do not explicitly distinguish
between $\hat{N}_{k}=\sharp\{\hat{G}_{k|K^{*}}\}$ and $N_{k}=\sharp\{G_{k|K^{*}}\}$.
Similarly, $\hat{\vartheta}_{(k|K)}=\left(\hat{\pi}_{(k|K)}^{\prime},\hat{\sigma}_{v(k|K)}\right)^{\prime}.$
Inspired by the pseudo log-likelihood, we construct an information
criterion to determine the optimal number of groups as follows:
\begin{align}
\text{IC}(K,\lambda_{NT}) & =\sum_{k=1}^{K}\left\{ N_{k}T\log(\hat{\sigma}_{v(k|K)})+\sum_{i\in\hat{G}_{k|K}}\sum_{t=1}^{T}\frac{\left(\ddot{y}_{it}-\underline{\ddot{z}}_{it}^{\prime}\hat{\pi}_{(k|K)}\right)^{2}}{\hat{\sigma}_{v(k|K)}^{2}}\right\} +\lambda_{NT}K\nonumber \\
& =\sum_{k=1}^{K}\left\{ N_{k}T\log(\hat{\sigma}_{v(k|K)})+N_{k}(T-1)\right\} +\lambda_{NT}K,\label{eq:IC_group}
\end{align}
where $\lambda_{NT}$ is a suitable penalty term. The optimal number of groups is the minimizer of (\ref{eq:IC_group})
\[
\hat{K}(\lambda_{NT})=\arg\min_{K=1,2,\ldots,\bar{K}}\text{IC}(K,\lambda_{NT}),
\]
given a suitable $\bar{K}$. For brevity, we henceforth refer to this as
$\hat{K}$. The final group estimates are then given by
\[
\hat{\vartheta}_{(k|\hat{K})}=\left(\hat{\pi}_{(k|\hat{K})},\hat{\sigma}_{v(k|\hat{K})}\right),\quad k=1,2,\ldots,\hat{K}.
\]
\subsubsection*{Step 4: Estimation of $\alpha^{0}$ and $\sigma_{u}^{2}$}
We estimate $\alpha^{0}$ and $\sigma_{u}^{2}$ pooling all observations
via maximum likelihood estimation (MLE). Specifically, the estimate is given by
\[
\left(\hat{\alpha}^{0},\hat{\sigma}_{u}^{2}\right)=\arg\max_{(s,\delta_{u}^{2})}\sum_{k=1}^{\hat{K}}\sum_{i\in\hat{G}_{k|\hat{K}}}\log f\left(y_{i}\mid x_{i};s,\delta_{u}^{2},\hat{\vartheta}_{(k|\hat{K})}\right),
\]
where $y_{i}=\left(y_{i1},\ldots,y_{iT}\right)^{\prime}$, $x_{i}=\left(x_{i1},\ldots,x_{iT}\right)^{\prime}$,
and $f\left(y_{i}\mid x_{i};s,\delta_{u}^{2},\hat{\vartheta}_{(k|\hat{K})}\right)$
is defined in (\ref{EQ:likelihood_appro}), noting that we use the post-classification
estimates of $\vartheta$. Since only two parameters, $\alpha^{0}$
and $\sigma_{u}^{2}$, are being estimated at this stage, the numerical
optimization is straightforward.
\subsection{Inefficiency Term}\label{SEC:inefficiency}
We now consider the general case where $\alpha^{0}$ can be heterogeneous
across $i$, and we will use $\alpha_{i}^{0}$ from this point onward.
Modeling the underlying structure of $\alpha_{i}^{0}-u_{i}$ differs
from that of $\left\{\alpha_{i}\left(\tau_{t}\right),\beta_{i}\left(\tau_{t}\right),\sigma_{vi}^{2}\right\}$
because $u_{i}$ is assumed to be random effects, and $\alpha_{i}^{0}-u_{i}$
naturally varies across $i$, even when $\alpha_{i}^{0}$ is identical.
For this reason, we focus on identifying the distribution of $\alpha_{i}^{0}-u_{i}$
rather than the actual values. While we can uncover the underlying
distribution, consistently estimating group membership remains challenging.
Consider the following example to illustrate this point. Suppose we
have two random variables $\varepsilon_{1}$ and $\varepsilon_{2}$,
with $\varepsilon_{1}\sim0-\left|N(0,2)\right|$ and $\varepsilon_{2}\sim1-\left|N(0,1)\right|$.
If we mix i.i.d. realizations of $\varepsilon_{1}$ and $\varepsilon_{2}$,
such as $\left\{ \varepsilon_{11},\varepsilon_{12},\ldots,\varepsilon_{1n},\varepsilon_{21},\varepsilon_{22},\ldots,\varepsilon_{2n}\right\} $,
it is likely that many $\varepsilon_{1i}$ and $\varepsilon_{2j}$ values lie very close to one another. For example, a small-scale Monte Carlo experiment
with $n=100$ suggests that about 28\% of $\varepsilon_{1i}$
have at least one $\varepsilon_{2j}$ within a radius of 0.01. In
such cases, swapping their memberships would likely have a minimal impact
on the likelihood function, complicating their distinct identification
from the data.
Misclassification of group memberships can have a serious impact on
the inefficiency term, unlike parameters at the frontiers, where only
similar frontiers can be misclassified together due to low power or
minor estimation errors. Continuing the previous example, suppose
$\varepsilon_{1i}=0$ (highly efficient with $u_{1i}=0$) is misclassified
as $\varepsilon_{2}$, then the inefficiency term for $\varepsilon_{1i}$
would be calculated as 1 (indicating inefficiency). Conversely, if
$\varepsilon_{2j}=0$ (originally not efficient with $u_{2j}=1$)
is misclassified as $\varepsilon_{1}$, then the inefficiency term
for $\varepsilon_{2j}$ would be calculated as 0 (indicating high
efficiency).
Given these challenges and the serious implications of misclassification,
we adopt a mixture distribution approach as follows. Suppose there exist an integer $\mathcal{K}^{*}\geq1$, such
that with probability $\tau_{j}^{0}$, it is distributed as $\alpha_{(j)}^{0}-\left|N(0,\sigma_{u(j)}^{2})\right|$
for $j=1,2,...,\mathcal{K}^{*}-1$, and with probability $\tau_{\mathcal{K}^{*}}^{0}=1-\tau_{1}^{0}-...-\tau_{\mathcal{K}^{*}-1}^{0}$,
as $\alpha_{(\mathcal{K}^{*})}^{0}-\left|N(0,\sigma_{u(\mathcal{K}^{*})}^{2})\right|$,
where $\left(\alpha_{(j)}^{0},\sigma_{u(j)}^{2}\right)$, $j=1,2,...,\mathcal{K}^{*}$,
are distinct vectors, $0<\tau_{j}^{0}<1,$ $j=1,2,...,\mathcal{K}^{*}-1,$
and $1-\tau_{1}^{0}-...-\tau_{\mathcal{K}^{*}-1}^{0}>0$. When $\mathcal{K}^{*}=1$,
the error distribution is reduced to that of a
unique distribution. The mixture distribution on the inefficiency
term is similar in spirit to the latent class model in \citet{Greene2005}.
However, the latent class model is only a small part of \citet{Greene2005} and so the treatment is very brief. We examine this issue in depth by proposing an information criterion to determine the
number of components, rigorously establish its theoretical properties,
and assess its small-sample performance through simulation studies.
An alternative approach is to model the composite term $\alpha_i^{0}-u_i$ by assuming that they follow a mixture distribution as a whole. However, without additional identifying restrictions, $\alpha_i^{0}$ and $u_i$ cannot be separately identified. As emphasized in the Introduction, separating these components is central to stochastic frontier analysis. For this reason, we keep our current specification.
The potential presence of a mixture distribution significantly alters
the interpretation of the results. With uniquely distributed inefficiency term, we
can remove the subscript $i$ from $\alpha_{i}^{0}$ because $\alpha_{i}^{0}-u_{i}\overset{d}{\sim}\alpha^{0}-\left|N\left(0,\sigma_{u}^{2}\right)\right|$.
A point estimate of $\alpha_{}^{0}-u_{i}$ is
\[
\widehat{\alpha^{0}-u_{i}}=\frac{1}{T}\sum_{t=1}^{T}\left(y_{it}-z_{it}'\hat{\pi}_{i}\right),
\]
where $\hat{\pi}_{i}$ is a sub-vector of $\hat{\tilde{\pi}}_{i}$
defined in (\ref{eq:pihat}) in Step 1. Thus, $\alpha^{0}-u_{i}$ can be estimated consistently as $T\to\infty$, allowing us to rank firms according to inefficiency
because $\alpha^{0}$ is identical across $i$. However, in the case
of a mixture distribution, although we can still consistently estimate
$\widehat{\alpha_{i}^{0}-u_{i}}$ (similar to the above), it is not
possible to rank firms as in the former scenario. This limitation
arises because memberships, or equivalently, the values of $\alpha_{i}^{0}$, cannot
be identified. This observation aligns with the findings for the cross-sectional
case discussed in \citet{Greene2005}.
The presence of a mixture distribution in the distribution of $\alpha_{i}^{0}-u_{i}$
does not impact the estimation of $\vartheta_{i}$ given independence
among $u_{i}$, $v_{it}$, and $x_{it}$. Consequently, Steps 1, 2,
and 3 remain unchanged. Details of the revised Step 4, now referred
to as Step 4', are provided below.
\subsubsection*{Step 4': Estimation of $\alpha_{(j)}^{0},\sigma_{u(j)}^{2},\text{ and }\tau_{j}^{0}$}
We adopt the mixture distribution for $\alpha_{i}^{0}-u_{i}$ as previously
described. Assuming that the inefficiency terms come from $\mathcal{K\geq}1$
distributions, we obtain an estimate of $\left(\alpha_{(1)}^{0},\sigma_{u(1)}^{2},...,\alpha_{(\mathcal{K})}^{0},\sigma_{u(\mathcal{K})}^{2},\tau_{1}^{0},...,\tau_{\mathcal{K}-1}^{0}\right)$
using MLE as follows:
\begin{align*}
& \left(\hat{\alpha}_{(1)}^{0},\hat{\sigma}_{u(1)}^{2},...,\hat{\alpha}_{(\mathcal{K})}^{0},\hat{\sigma}_{u(\mathcal{K})}^{2},\hat{\tau}_{1},...,\hat{\tau}_{\mathcal{K}-1}\right)\\
= & \arg\max_{(s,\delta_{u}^{2},\tau)}\sum_{k=1}^{\hat{K}}\sum_{i\in\hat{G}_{k|\hat{K}}}\log\tilde{f}\left(y_{i}\left\vert x_{i};s_{(1)},\delta_{u\left(1\right)}^{2},...,s_{(\mathcal{K})},\delta_{u\left(\mathcal{K}\right)}^{2},\tau_{1},...,\tau_{\mathcal{K}-1},\hat{\vartheta}_{(k|\hat{K})}\right.\right)
\end{align*}
where $\tilde{f}$ is the likelihood function defined in (\ref{EQ:log_approx_p}),
and we incorporate estimates from Step 3, as detailed in Section \ref{SEC:Estimation}.
\subsubsection*{Step 5: Determination of the Distributional Structures of the Inefficiency
Term}
To determine the optimal number of mixtures, we introduce a new information criterion for this task:
\begin{equation}
\widetilde{\mathrm{IC}}(\mathcal{K},\tilde{\lambda}_{NT})=-\sum_{k=1}^{\hat{K}}\sum_{i\in\hat{G}_{k|\hat{K}}}\log\tilde{f}\left(y_{i}\left\vert x_{i};\hat{\alpha}_{(1)}^{0},\hat{\sigma}_{u(1)}^{2},...,\hat{\alpha}_{(\mathcal{K})}^{0},\hat{\sigma}_{u(\mathcal{K})}^{2},\hat{\tau}_{1},...,\hat{\tau}_{\mathcal{K}-1},\hat{\vartheta}_{(k|\hat{K})}\right.\right)+\mathcal{K}\tilde{\lambda}_{NT},\label{eq:IC_tide2}
\end{equation}
where $\tilde{\lambda}_{NT}$ is a suitable penalty term, and the
estimates are as obtained from Steps 3 and 4'.
The optimal number of mixtures is the minimizer of (\ref{eq:IC_tide2})
\[
\hat{\mathcal{K}}(\tilde{\lambda}_{NT})=\arg\min_{\mathcal{K}=1,2,\ldots,\bar{\mathcal{K}}}\widetilde{\mathrm{IC}}(\mathcal{K},\tilde{\lambda}_{NT}),
\]
and we write $\hat{\mathcal{K}}$ for short. Finally, the estimated parameters are
\[
\left(\hat{\alpha}_{(1)}^{0},\hat{\sigma}_{u(1)}^{2},...,\hat{\alpha}_{(\mathcal{\hat{\mathcal{K}}})}^{0},\hat{\sigma}_{u(\mathcal{\hat{\mathcal{K}}})}^{2},\hat{\tau}_{1},...,\hat{\tau}_{\mathcal{\hat{\mathcal{K}}}-1}\right).
\]
\subsection{A Summary of the Estimation Procedure}
The outline of the estimation procedure is as follows:
\begin{enumerate}
\item Conduct estimations of the frontiers for each firm as described in
Step 1.
\item Apply the HAC algorithm using the individual estimations
from Step 1 for $K=1,2,\ldots,\bar{K}$.
\item Use the information criterion in (\ref{eq:IC_group}) from Step 3
to determine the optimal number of groups and group memberships.
Obtain an estimate of the frontier with the determined groups.
\item Based on the group assignments from Step 3, perform joint estimation
for the inefficiency term as in Step 4'.
\item Use the information criterion in (\ref{eq:IC_tide2}) to determine
the distribution of $\alpha_{i}^{0}-u_{i}$, as in Step~5.
\item Collect results. The estimates of the frontiers and the distribution
of the inefficiency term are derived from Steps~3 and 5, respectively.
\end{enumerate}
\section{Asymptotic Properties} \label{SEC:Theoretical}
We examine classification consistency in Section \ref{SEC:Classification}.
Subsequently, we discuss the large sample properties of the post-classification
estimators in Section \ref{SEC:post}.
\subsection{Classification}
\label{SEC:Classification}
\begin{assumption}\label{A:mixing} The process $\{(x_{it}',v_{it}),t=1,\ldots,T\}$
is strong mixing with a mixing coefficient $\alpha(j)$ that satisfies
$\alpha(j)\leq C_{\alpha}\rho^{j}$ for some positive $C_{\alpha}$
and $0<\rho<1$, and this holds for $i=1,\ldots,N$. \end{assumption}
\begin{assumption}\label{A:moment} $\max_{1\leq i\leq N,1\leq t\leq T}\text{\emph{E}}\|x_{it}\|^{q}\leq\bar{C}_{x}<\infty,$
and $\max_{1\leq i\leq N,1\leq t\leq T}\text{\emph{E}}|v_{it}|^{q}\leq\bar{C}_{v}<\infty,$
for some $q>4.$ \end{assumption}
\begin{assumption}\label{A:rank} Denote $\tilde{x}_{it}\equiv(1,x_{it}^{\prime})^{\prime}$. Let
$\mu_{\min}$ and $\mu_{\max}$ denote the minimum and maximum eigenvalues
of a matrix, respectively. There exist $\underline{C}_{xx}$ and $\bar{C}_{xx}$ with $0<\underline{C}_{xx} \leq \bar{C}_{xx} <\infty$ such that
\[
0<\underline{C}_{xx}\leq\min_{1\leq i\leq N,1\leq t\leq T}\mu_{\min}[\text{\emph{E}}(\tilde{x}_{it}\tilde{x}_{it}')]\leq\max_{1\leq i\leq N,1\leq t\leq T}\mu_{\max}[\text{\emph{E}}(\tilde{x}_{it}\tilde{x}_{it}^{\prime})]\leq\bar{C}_{xx}<\infty.
\]
\end{assumption}
\begin{assumption}\label{A:coef} For $k=1,2,\ldots,K^{*}$, $\alpha_{(k)}^{*}(s)$,
$\beta_{(k)1}^{*}(s)$, ..., and $\beta_{(k)p}^{*}(s)$ belong to
$L^{2}\left[0,1\right]$ and are $\kappa$ times continuously differentiable.
\end{assumption}
\begin{assumption}\label{A:group_diff} There exists a positive $\underline{C}^{*}$,
such that
\[
\min_{1\leq j\neq k\leq K^{*}}\left\{ \|\alpha_{(j)}^{*}-\alpha_{(k)}^{*}\|+\sum_{l=1}^{p}\|\beta_{(j)l}^{*}-\beta_{(k)l}^{*}\|+|\sigma_{v(j)}^{*}-\sigma_{v(k)}^{*}|\right\} \geq\underline{C}^{*}>0,
\]
and $\min_{1\leq k\leq K^{*}}\sigma_{v(k)}^{*2}>0.$ \end{assumption}
\begin{assumption}\label{A:tuningPara} (i) $N$ can either be a finite,
or divergent. If $N$ is divergent, $N=O(T^{C})$ for some positive $C$ as $T\rightarrow\infty$. (ii) $m\rightarrow\infty$ as $T\rightarrow\infty.$
$Nm^{q/2+2}(\log N)^{2q}/T^{q/2-1}\rightarrow0,$ with $q$ in
Assumption \ref{A:moment}. \end{assumption}
Assumption \ref{A:mixing} imposes a condition of weak dependence
across $t$, noting that independence across $i$ is not required
for classifications. Assumption \ref{A:moment} requires that $x$ and
$v$ have finite $q$-th moment. Assumption \ref{A:rank} is the classic
full rank condition. Assumption \ref{A:coef} stipulates that the
coefficients are $\kappa$-th order differentiable, a standard condition
for nonparametric or semiparametric estimation. Assumption \ref{A:group_diff}
requires that at least one of the coefficients, including the variance
of $v$, must differ across groups.
Assumption \ref{A:tuningPara} specifies that the moment conditions
must be sufficiently large or that $T$ grows fast enough. $N$ can
be fixed. If $N$ diverges, it cannot be too fast, e.g., at the rate
of $\exp(T)$. In the empirical application, $(N,T)=(466,80)$
and $466\approx80^{1.4}$, thus any $C\geq1.4$ works in (i). The
most stringent requirements arise from the estimation of the ``design''
matrix $\frac{1}{T}Z_{im}^{\prime}Z_{im}$ with diverging dimensions,
used in $\widehat{\tilde{\pi}}_{i}$ (see (\ref{eq:pihat})); a similar
condition was imposed in \citet{Chen2019}. We take a logarithm of
all covariates before estimation, and $q$ can be reasonably considered
large, e.g., $q\geq8$. If we set $m=T^{1/5}$, Assumption \ref{A:tuningPara}
is satisfied. We do not have the usual bias and variance tradeoff
here, as explained below. The results of Theorem \ref{TH:classify}
are underpinned by the uniform convergence of $\widehat{\tilde{\pi}}_{i}$
without any rate requirement. As a result, it is not necessary to
consider the trade-off between the bias and variance of the estimates
for this aspect, when deciding $m$. Of course, we do need $m\rightarrow\infty$ to ensure the uniform convergence. However, to achieve the optimal convergence rate for the post-classification estimates, this consideration becomes crucial, as reflected in Assumption \ref{A:tuning2} (ii) in the subsequent section. With the aforementioned technical conditions, we demonstrate the consistency
of the classification.
\begin{theorem}\label{TH:classify} Suppose Assumptions \ref{A:mixing}
through \ref{A:tuningPara} hold. Then:
\noindent (i) For any small positive $\epsilon$,
\[
\Pr\left(\max_{i=1,2,\ldots,N}\left\Vert \hat{\vartheta}_{i}-\vartheta_{i}\right\Vert >\epsilon\right)=o(1);
\]
(ii) Assuming $K=K^{*}$, denote the event
\[
\mathcal{M}\equiv\left\{ \left(\hat{G}_{1|K^{*}},\hat{G}_{2|K^{*}},\ldots,\hat{G}_{K^{*}|K^{*}}\right)=\left(G_{1|K^{*}},G_{2|K^{*}},\ldots,G_{K^{*}|K^{*}}\right)\right\} ,
\]
then $\Pr(\mathcal{M})\rightarrow1$. \end{theorem}
Theorem \ref{TH:classify} (i) establishes the uniform convergence
of $\hat{\vartheta}_{i}$ for $i=1,2,\ldots,N$ provided $m\to\infty$. Building on this,
Theorem \ref{TH:classify} (ii) demonstrates that the probability
of correct classification approaches 1, provided that $K=K^{*}$.
In the next section, we will argue that the $K$ we choose converges
to $K^{*}$ with probability approaching 1, and we discuss the asymptotic
properties of the post-classification estimates.
We note that significantly fewer assumptions are required for consistency in classification compared to those needed for post-classification and determining the number of groups. For instance, we do not need independence or weak
dependence across $i$, nor do we require specific distributional assumptions on $v_{it}$ and $u_{i}$.
\subsection{The Number of Groups and Post-Classification Estimator}
\label{SEC:post}
In this section, we address the question regarding the choice of $K$ and the post-classification estimation. We begin by presenting additional assumptions necessary for this analysis.
\begin{assumption}\label{A:additional-1} $(x_{i},\varepsilon_{i})$,
for $i=1,2,\ldots,N$ are independent across $i$.
\end{assumption}
\begin{assumption}\label{A:group} $N_{k}\propto N$ for each $k=1,2,\ldots,K^{*}$.
\end{assumption}
\begin{assumption}\label{A:stochasticF} $v_{it}\overset{iid}{\sim}N\left(0,\sigma_{vi}^{2}\right)$ across $t.$ There exist an integer $\mathcal{K}^{*}\geq1$, such that,
\[
\alpha_{i}^{0}-u_{i}\overset{d}{\sim}\alpha_{(j)}^{0}-\left\vert N\left(0,\sigma_{u(j)}^{2}\right)\right\vert \textrm{ with probability }\tau_{j}^{0},j=1,2,...,\mathcal{K}^{*},
\]
where $0<\tau_{j}^{0}<1$ and $\sigma_{u(j)}^{2}>C>0$ for $j=1,2,...,\mathcal{K}^{*},$
$\tau_{1}^{0}+\tau_{2}^{0}+...+\tau_{\mathcal{K}^{*}}^{0}=1,$ $\left(\alpha_{(j)}^{0},\sigma_{u(j)}^{2}\right)$
differ across $j=1,2,...,\mathcal{K}^{*}$. Lastly, the sequences
$\{v_{it}\}_{t=1}^{T}$, $u_{i}$, and $\{x_{it}\}_{t=1}^{T}$ are
mutually independent. \end{assumption}
\begin{assumption}\label{A:tuning2} (i) $T\propto N^{C_{*}}$ for
some positive $C_{*}$ as $N\to\infty$; (ii) $\underline{m}\to\infty$
as $T,N\to\infty$. Additionally, $\underline{m}/T\to0$, $\underline{m}^{q/2+2}(\log N)^{2q}/(NT)^{q/2-1}\to0$,
and $NT/\underline{m}^{1+2\kappa}\to0$; (iii) $C_{*}>1/(2\kappa)$
and $(q-2)\kappa>3$. \end{assumption}
Assumption \ref{A:additional-1} further imposes independence across
$i.$ Assumption \ref{A:group} says the number of members in each
group is proportional to $N$. This condition is not necessary, but
it facilitates expositions. Assumption \ref{A:stochasticF} specifies the
distributional conditions on the error terms. As explained in Section \ref{SEC:inefficiency}, these conditions are essential for the identification of $\sigma_{u}^{2}$, and common in the literature, (see, for e.g., \citet{YaoZhangKum2019}).
Assumption \ref{A:tuning2} places
restrictions on the rate of growth of $T$ relative to $N$, and the rate of growth of the tuning parameter $\underline{m}$. The condition in (iii) is set to ensure that the set of $\underline{m}$ that satisfies (ii) is
not empty. We need $\underline{m}/T\rightarrow0$ so that the ``design''
matrix $\text{E}\left(\tilde{z}_{it}\tilde{z}_{it}'\right)$ is still
well-behaved, as required in Lemma \ref{LE:splines}. However, condition \textit{$\underline{m}/T\rightarrow0$}
can be restrictive when $N$ is much larger than $T$. $\left.\underline{m}^{q/2+2}\left(\log N\right)^{2q}\right/\left(NT\right)^{q/2-1}\rightarrow0$
is assumed to ensure the sample version of the design matrix, namely $Q_{(k),zz}$
defined in (\ref{eq:Qkzz}), is of full rank with very high probability.
Note that the dimension of $Q_{(k),zz}$ is diverging, so the consistency
of this matrix requires uniform convergence of all elements and hence
this restriction. $NT\left/\underline{m}^{1+2\kappa}\right.\rightarrow0$
ensures the bias term is asymptotically negligible. In the special
case where $\kappa\geq2,$ we need $C_{*}>1/4,$ and $\left(q-2\right)\kappa>3$
is satisfied due to $q>4$ in Assumption \ref{A:moment}. One can
then set, for example, $\underline{m}=\left(NT\right)^{1/4.8}$, which
satisfies condition (ii).
We show the asymptotic properties of our estimators for the case where
$\mathcal{K}^{*}\geq2$. The case in which $\alpha_{i}^{0}-u_{i}$ comes
from a unique distribution is straightforward given this result.
Recall that $\mathbb{B}_{-0}^{m}\left(\tau_{t}\right)'\pi_{i0}^{0}$
and $\mathbb{B}^{m}\left(\tau_{t}\right)^{\prime}\pi_{il}^{0}$ represent
the approximations of $\alpha_{i}\left(\tau_{t}\right)$ and $\beta_{il}\left(\tau_{t}\right).$
For the coefficients on the frontiers, let $\theta\left(s\right)\equiv\left(\alpha\left(s\right),\beta\left(s\right)'\right)',$
and correspondingly $\hat{\theta}\left(s\right)=\left(\mathbb{B}_{-0}^{\underline{m}}\left(\tau_{t}\right)'\hat{\pi}_{0},\mathbb{B}^{\underline{m}}\left(\tau_{t}\right)^{\prime}\hat{\pi}_{1},...,\mathbb{B}^{\underline{m}}\left(\tau_{t}\right)^{\prime}\hat{\pi}_{p}\right)'$.
For the parameters in the distribution of $\alpha_{i}^{0}-u_{i}$,
we denote
\[
\varrho^{0}\equiv\left(\alpha_{(1)}^{0},\sigma_{u(1)}^{2},...,\alpha_{(\mathcal{K^{*}})}^{0},\sigma_{u(\mathcal{K}^{*})}^{2},\tau_{1}^{0},...,\tau_{\mathcal{\mathcal{K}^{*}}-1}^{0}\right),
\]
and correspondingly
\[
\hat{\varrho}=\left(\hat{\alpha}_{(1)}^{0},\hat{\sigma}_{u(1)}^{2},...,\hat{\alpha}_{(\mathcal{K}^{*})}^{0},\hat{\sigma}_{u(\mathcal{\mathcal{K}^{*}})}^{2},\hat{\tau}_{1},...,\hat{\tau}_{\mathcal{\mathcal{K}^{*}}-1}\right).
\]
The following notations are used to characterize the asymptotic distribution. Denote
\[
\mathbb{M_{B}}\left(s\right)\equiv\left(\begin{array}{cccc}
\mathbb{B}_{-0}^{\underline{m}}\left(s\right)' & 0 & \cdots & 0\\
0 & \mathbb{B}^{\underline{m}}\left(s\right)' & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & \mathbb{B}^{\underline{m}}\left(s\right)'
\end{array}\right)_{\left(p+1\right)\times\left(\underline{m}-1+\underline{m}p\right)},
\]
\begin{equation}
Q_{(k),zz}=\frac{1}{N_{k}T}\sum_{i\in G_{k|K^{*}}}\sum_{t=1}^{T}\underline{\ddot{z}}_{it}\underline{\ddot{z}}_{it}',\label{eq:Qkzz}
\end{equation}
and
\[
\mathbb{S}_{(k)}\left(s\right)=\frac{\sigma_{v\left(k\right)}^{*2}}{\underline{m}}\mathbb{M_{B}}\left(s\right)Q_{(k),zz}^{-1}\mathbb{M_{B}}\left(s\right)',
\]
noting that $\mathbb{S}_{(k)}\left(s\right)$ is positive definite with very high probability; which we show it in equation (\ref{eq:var(Ak2)}) of Appendix \ref{APP:lemmaProofs}. For notation convenience, write
\[
\tilde{f}_{i\left(k\right)}\left(\varrho\right)\equiv\tilde{f}\left(y_{i}\left\vert x_{i};\varrho,\vartheta_{\left(k|K^{*}\right)}\right.\right),
\]
and
\[
\mathbb{I}\equiv\left.-\text{E}\left[\frac{1}{N}\sum_{k=1}^{K^{*}}\sum_{i\in G_{k}}\frac{\partial^{2}}{\partial\varrho\partial\varrho'}\log\tilde{f}_{i\left(k\right)}\left(\varrho\right)\right]\right|_{\varrho=\varrho^{0}},
\]
with $\mathbb{I}^{1/2}$ denoting the matrix such that $\mathbb{I}^{1/2}\mathbb{I}^{1/2\prime}=\mathbb{I}$.
$\mathbb{I}$ is a positive definite matrix with finite eigenvalues,
as shown in Appendix \ref{App:derivatives}. As before, we
show the asymptotic property of $\hat{\theta},\hat{\sigma}_{v}^{2}$ and
$\hat{\varrho}$ pretending that we know $K^{*}$ and $\mathcal{K}^{*}.$
We then show that $\hat{K}$ and $\hat{\mathcal{K}}$ converge to
$K^{*}$ and $\mathcal{K}^{*}$, respectively, with probability approaching
1.
\begin{theorem}\label{TH:post-estimation} Suppose Assumptions \ref{A:mixing}
through \ref{A:tuning2} hold, $\hat{K}(\lambda_{NT})=K^{*}$ and
$\hat{\mathcal{K}}\left(\tilde{\lambda}_{NT}\right)=\mathcal{K}^{*}$.
Let $I_{l}$ denote the $l\times l$ identity matrix. Then, for each
$k=1,2,\ldots,K^{*}$,
\noindent (i)
\[
\sqrt{\frac{N_{k}T}{\underline{m}}}\mathbb{S}_{(k)}^{-1/2}(s)\left(\hat{\theta}_{(k|K^{*})}(s)-\theta_{(k)}^{*}(s)\right)\overset{d}{\rightarrow}N(0,I_{p+1});
\]
(ii)
\[
\sqrt{N_{k}T}\left(\hat{\sigma}_{v(k|K^{*})}^{2}-\sigma_{v(k)}^{*2}\right)\overset{d}{\rightarrow}N\left(0,\text{\emph{Var}}(v_{it}^{2}|i\in G_{k|K^{*}})\right);
\]
(iii)
\[
\sqrt{N}\mathbb{I}^{-1/2}(\hat{\varrho}-\varrho^{0})\overset{d}{\rightarrow}N(0,I_{3\mathcal{K}^{*}-1}).
\]
\end{theorem}
This theorem establishes the asymptotic properties of the post-classification estimators. As expected, $\hat{\theta}$ converges at a nonparametric
rate, while $\hat{\sigma}_{v}^{2}$ converges at a parametric rate.
The convergence rate of $\hat{\varrho}$ is $\sqrt{N}$ and does not
depend on $T$. This result may appear odd, but there is a simple explanation. Note that
$\varrho$ collects only the parameters that govern the distribution of $u_{i}$. The best scenario of estimating $\varrho$ is that we observe $u_{1},u_{2},...,u_{N}$ directly, in which case the rate of convergence of $\hat{\varrho}$ is $\sqrt{N}$. In theory, the value of $T$ does not impact the convergence rate of $\hat{\varrho}$. However, in finite samples, large $T$ can potentially ensure a more precise estimation of $u_{i}$, and thus can possibly improve the finite-sample performance of $\hat{\varrho}$.
In addition, the validity of the proposed information criteria relies
on these properties, as they depend on the accuracy and consistency
of the post-classification estimators, as demonstrated above. For example,
$\lambda_{NT}$ depends on both $N$ and $T$ (due to the rates of
$\hat{\theta}$ and $\hat{\sigma}_{v}^{2}$), while $\tilde{\lambda}_{NT}$
depends only on $N$ (due to the rate of $\hat{\varrho}$).
\begin{proposition}\label{Prop:classify} Suppose Assumptions \ref{A:mixing} through \ref{A:tuning2} hold.
\noindent (i) Select a value of $\lambda_{NT}$ such that $(NT)^{-1/2}\lambda_{NT}\rightarrow\infty$
and $(NT)^{-1}\lambda_{NT}\rightarrow0$. Then,
\[
\Pr\left(\hat{K}\left(\lambda_{NT}\right)=K^{*}\right)\rightarrow1.
\]
\noindent (ii) Select a value of $\tilde{\lambda}_{NT}$ such that $\tilde{\lambda}_{NT}\rightarrow\infty$
and $N^{-1}\tilde{\lambda}_{NT}\rightarrow0$. Then,
\[
\Pr\left(\hat{\mathcal{K}}\left(\tilde{\lambda}_{NT}\right)=\mathcal{K}^{*}\right)\rightarrow1.
\]
\end{proposition}
Proposition \ref{Prop:classify} presents conditions under which the
information criteria are valid, focusing on the tuning parameters
$\lambda_{NT}$ and $\tilde{\lambda}_{NT}$. As highlighted earlier,
selecting the correct range for $\lambda_{NT}$ and $\tilde{\lambda}_{NT}$
is crucial. In the subsequent section, we will evaluate specific values
for $\lambda_{NT}$ and $\tilde{\lambda}_{NT}$, identify those that
perform well in simulations, and recommend practical choices.
\section{Monte Carlo Simulations}
\label{SEC:simulations}
\subsection{Simulation Designs}
Heterogeneity in panel SF models arises from various
sources. Thus, we design three different Monte Carlo experiments that
allow us to examine the finite-sample performance of the proposed
method and its ability to identify sources of heterogeneity. In
the first design, we study the classification for the case with heterogeneity from frontiers yet with constant variances of $v_{it}$. In the second design, we study the case with heterogeneous variances of $v_{it}$, yet homogeneous frontiers. In the third design, we check the performance of our methods in a more general and much more complicated scenario. In all three designs, we consider two sub-cases where
the term $\alpha^{0}-u$ comes from either unique or from mixture
distribution, which previous methods did not consider. Due to the similarity and the space constraint, we defer the description and simulation results of Designs 1 and 2 to Appendix \ref{APP:Tables} and only present Design 3 here.
\textbf{Design 3:} In our third design (consisting of DGP3U and DGP3M),
we study the performance of our method in a setting similar to those
in \citet{YaoZhangKum2019}, where there are three groups for both
the frontiers and variances, with two regressors. The DGP is
\[
y_{it}=\alpha_{i}^{0}-u_{i}+\alpha_{i}(\tau_{t})+x_{it1}\beta_{i1}(\tau_{t})+x_{it2}\beta_{i2}(\tau_{t})+v_{it},
\]
where $x_{itl}\sim N(1,0.5^{2})$ for both regressors $l=1,2$.
Group 1 frontiers and error term are specified as $\alpha_{(1)}(s)=-\frac{1}{1+3s}-\varpi_{1}$,
$\beta_{(1)1}(s)=2s^{3}$, $\beta_{(1)2}(s)=\ln(5s)$, $v_{it}\overset{iid}{\sim}N(0,\sigma_{v(1)}^{2})$
with $\sigma_{v(1)}=0.75$ and $\varpi_{1}$ is a mean of $-\frac{1}{1+3s}$.
Group 2 frontiers and error term are specified as $\alpha_{(2)}(s)=-\cos(4s)-\varpi_{2}$,
$\beta_{(2)1}(s)=\sin(4s)$, $\beta_{(2)2}(s)=\ln(\frac{s}{1-s})$,
$v_{it}\overset{iid}{\sim}N(0,\sigma_{v(2)}^{2})$ with $\sigma_{v(2)}=1.25$
and $\varpi_{2}$ is a mean of $-\cos(4s)$. Group 3 frontiers and
error term are specified as $\alpha_{(3)}(s)=5s^{2}-s+1-\varpi_{3}$,
$\beta_{(3)1}(s)=\exp{(-s)}+\sin(5s)$, $\beta_{(3)2}(s)=-5\sin(s)\cos(5s)+1$,
$v_{it}\overset{iid}{\sim}N(0,\sigma_{v(3)}^{2})$, with $\sigma_{v(3)}=1.25$
and $\varpi_{3}$ is a mean of $5s^{2}-s+1$. We consider two
sub-cases of $\alpha^{0}-u$, which we denote
them as DGP3U and DGP3M. In DGP3U, $\alpha^{0}-u$
comes from a unique distribution, with $\alpha^{0}=0.5$ and $u_{i}\overset{iid}{\sim}|N(0,\sigma_{u}^{2})|$,
where $\sigma_{u}=1$. In DGP3M, we let $\alpha^{0}-u$ to come from $\alpha_{(1)}^{0}-|N(0,\sigma_{u(1)}^{2})|$
with probability $\tau^{0}$ and $\alpha_{(2)}^{0}-|N(0,\sigma_{u(2)}^{2})|$
with probability $1-\tau^{0}$, where $\alpha_{(1)}^{0}=1$, $\alpha_{(2)}^{0}=-1$,
$\sigma_{u(1)}=0.75$, $\sigma_{u(2)}=1.25$ and $\tau^{0}=0.5$. It is important to note that the mixture structure of $\alpha^{0}-u$ is
independent of the grouping structure.
We evaluate the performance of each model and the case for any combination
of $N=100,250,$ or $500$ and $T=50,75,$ or $100$. Thus,
there are $3\times2\times9=54$ different cases. We assess the finite
sample properties of our method with 500 MC replications.
Note $(N,T)=(466,80)$ in the empirical application of the paper,
so our simulations, including the recommended tuning parameters in
the next section, offer meaningful and practical guidance.
\subsection{Choices of Tuning Parameters} \label{SEC:tuning}
We set $m=\left\lfloor T^{1/5}\right\rfloor $, where $\left\lfloor \cdot\right\rfloor$ denotes the integer part and $\underline{m}=\left\lfloor \left(N_{k}T\right)^{1/4.8}\right\rfloor $
for each group $k$. The value of the two tuning parameters align with standard choices in the literature.
Theoretically, the valid ranges for $\lambda_{NT}$ and $\tilde{\lambda}_{NT}$
are quite broad. Based on simulation evidence, we recommend setting
$\lambda_{NT}=\left(c_{\lambda}\sqrt{NT}\log(NT)\right)/2$ and $\tilde{\lambda}_{NT}=\left(\tilde{c}_{\lambda}\sqrt{N}\log N\right)/8$,
where $c_{\lambda}$ and $\tilde{c}_{\lambda}$ are constants. These
values of $\lambda_{NT}$ and $\tilde{\lambda}_{NT}$ meet the conditions
specified in Proposition \ref{Prop:classify}. The constants $c_{\lambda}$
and $\tilde{c}_{\lambda}$ serve as sensitivity parameters, over which
we conduct sensitivity analyses for the choice of $\lambda_{NT}$
and $\tilde{\lambda}_{NT}$. Specifically, we test values of $c_{\lambda}$
and $\tilde{c}_{\lambda}$ in $\{3/2,1,3/4\}$, with $c_{\lambda}=\tilde{c}_{\lambda}=1$
as the benchmark setting.
In finding the number of components in the mixture distribution, we restrict our attention to one or two components, i.e., $\mathcal{K}=1,2$. We find that the small-sample properties of mixtures with more than two components perform poorly in simulations. We conjecture that this is due to the complicated log-likelihood functions
(see, e.g., equations (\ref{EQ:likelihood_appro}) and (\ref{EQ:log_approx_p})), which cannot effectively
handle mixtures with more than two components. In Appendix \ref{APP:more_mix}, we propose an alternative method to
identify the mixture structure with potentially more than two components.
We also conduct simulations to evaluate its finite-sample performance.
The simulation results for this alternative method for allowing mixtures with more than two components suggest that the method performs well in determining the correct number of components, but the parameter estimates can be off. Although not perfect, these findings provide a foundation for future research in this direction.
\subsection{Simulation Results}
We report results for DGP3M, the most complex model, featuring three
groups and a mixture distribution structure in $\alpha_{i}^{0}-u_{i}$. Results for the remaining DGPs are reported in Appendix \ref{APP:Tables}.
We first report the performance of the IC in Step 3 for coefficient
groups and Step 5 for the mixture distribution structure in Table \ref{tab:IC_DGP3M} for the benchmark specification with $c_{\lambda}=\tilde{c}_{\lambda}=1$. Additionally, Table \ref{tab:IC_DGP3M} includes the classification errors, denoted as $\bar{\text{Pr}}(\bar{F})$ . It is defined as the average percentage of observations misclassified to other groups across 500 replications.
The performance of the IC in Step 3 for selecting the correct number
of groups, $K^{*}=3$, is reasonable. For $N=500$, the classification
error in Step 3 is less than 1 percent for each $T=50,75,100$. The
performance of the Step 5 IC is also strong for DGP3M, choosing the
correct specification (mixture distribution) with a probability close
to 1. For DPG3U, where $\alpha_{i}^{0}-u_{i}$ comes from a unique
distribution, Table \ref{tab:IC_DGP3U} in Appendix \ref{APP:Tables}
shows that the probability of Step 5 IC selecting the correct distribution (unique distribution) quickly approaches 1 as $N$ increases. Sensitivity analyses for both ICs in Steps 3 and 5, shown in Tables \ref{tab:Sensitivity_Step3_IC_DGP1U}
- \ref{tab:Sensitivity_Step5_IC_DGP3M}, demonstrate that the results
are robust to the selected range of tuning parameters.
We assess the accuracy of the estimates of $\{\sigma_{v}\text{s},\alpha_{1}^{0},\sigma_{u(1)},\alpha_{2}^{0},\sigma_{u(2)},\tau^{0}\}$
using two measures: (i) bias (BIAS) and (ii) root mean squared errors
(RMSE). The reported values in Table \ref{tab:DGP3M_parameters},
obtained by averaging over 500 MC iterations, are reasonable.
Illustrated in Figures \ref{fig:Grouped_frontiers_DGP6_N500_T50},
\ref{fig:Grouped_frontiers_DGP6_N500_T75} and \ref{fig:Grouped_frontiers_DGP6_N500_T100}
in Appendix \ref{APP:Tables} are the estimates of time-varying frontiers
for $(N,T)=(500,50),(500,75),\text{ and }(500,100)$. Black solid lines
depict the true time-varying frontier, dotted lines show the mean
of the estimated grouped frontiers averaged over 500 MC iterations,
and the gray shaded region depicts the 90th percentile of the estimates.
It is clear from Figure \ref{fig:Grouped_frontiers_DGP6_N500_T50}
that, while the mean over MC iterations is reasonably close to the
true frontiers, the 90th percentile bands are wide, suggesting possible
classification errors between neighboring groups. Figures \ref{fig:Grouped_frontiers_DGP6_N500_T75}
and \ref{fig:Grouped_frontiers_DGP6_N500_T100} show that the accuracy
of frontier grouping improves as $T$ increases. This is consistent
with the theory developed, since the Step 2 classification using HAC
is based on $\hat{\vartheta}_{i}=\left(\hat{\pi}_{i}^{\prime},\hat{\sigma}_{vi}\right)^{\prime}$
obtained using $T$ observations.
\begin{table}[!htb]
\centering
\caption{Performance of ICs for DGP3M}
\vspace{0.2cm}
\begin{tabularx}{\textwidth}{X XXXX X XX}
\toprule
$(N,T)$ & $K=1$ & $K=2$ & $K=3$ & $K=4$ & $\bar{\text{Pr}}(\bar{F})$ & $\alpha^0-u$ uni & $\alpha^0-u$ mix \\
\midrule
(100,50) & 0.000 & 0.352 & 0.648 & 0.000 & 0.106 & 0.008 & 0.992 \\
(100,75) & 0.000 & 0.078 & 0.922 & 0.000 & 0.024 & 0.000 & 1.000 \\
(100,100) & 0.000 & 0.000 & 1.000 & 0.000 & 0.024 & 0.000 & 1.000 \\
(250,50) & 0.000 & 0.086 & 0.914 & 0.000 & 0.026 & 0.002 & 0.998 \\
(250,75) & 0.000 & 0.000 & 1.000 & 0.000 & 0.026 & 0.000 & 1.000 \\
(250,100) & 0.000 & 0.000 & 1.000 & 0.000 & 0.026 & 0.000 & 1.000 \\
(500,50) & 0.000 & 0.006 & 0.994 & 0.000 & 0.002 & 0.004 & 0.996 \\
(500,75) & 0.000 & 0.000 & 1.000 & 0.000 & 0.002 & 0.000 & 1.000 \\
(500,100) & 0.000 & 0.000 & 1.000 & 0.000 & 0.002 & 0.000 & 1.000 \\
\bottomrule
\end{tabularx}
\begin{tablenotes}
\footnotesize
\item
\emph{Note}: Results for the baseline case $c_{\lambda}=\tilde{c}_{\lambda}=1$. Reported numbers are probabilities across replications.
\end{tablenotes}
\label{tab:IC_DGP3M}
\end{table}
\begin{table}[!htb]
\centering
\caption{BIAS and RMSE over 500 MC iterations for DGP3M}
\vspace{0.2cm}
\scalebox{0.71}{
\begin{tabularx}{1.4\textwidth}{X c XX c XX c XX c XX c XX c XX c XX c XX}
\toprule
&& \multicolumn{2}{c}{$\hat{\sigma}_{v(1)}$} && \multicolumn{2}{c}{$\hat{\sigma}_{v(2)}$} && \multicolumn{2}{c}{$\hat{\sigma}_{v(3)}$} && \multicolumn{2}{c}{$\hat{\tau}$} && \multicolumn{2}{c}{$\hat\alpha^0_{(1)}$} && \multicolumn{2}{c}{$\hat{\sigma}_{u(1)}$} && \multicolumn{2}{c}{$\hat\alpha^0_{(2)}$} && \multicolumn{2}{c}{$\hat{\sigma}_{u(2)}$} \\
\cmidrule{3-4} \cmidrule{6-7} \cmidrule{9-10} \cmidrule{12-13} \cmidrule{15-16} \cmidrule{18-19} \cmidrule{21-22} \cmidrule{24-25}
$(N,T)$ && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE && BIAS & RMSE\\
\midrule
(100,50) && 0.268 & 0.463 && 0.127 & 0.215 && 0.367 & 0.594 && 0.021 & 0.029 && 0.057 & 0.074 && 0.116 & 0.152 && 0.120 & 0.157 && 0.135 & 0.171 \\
(100,75) && 0.082 & 0.214 && 0.074 & 0.174 && 0.154 & 0.372 && 0.015 & 0.027 && 0.042 & 0.052 && 0.103 & 0.145 && 0.100 & 0.131 && 0.126 & 0.162 \\
(100,100) && 0.024 & 0.099 && 0.027 & 0.093 && 0.053 & 0.196 && 0.013 & 0.018 && 0.040 & 0.052 && 0.090 & 0.116 && 0.091 & 0.118 && 0.108 & 0.137 \\
(250,50) && 0.204 & 0.373 && 0.138 & 0.243 && 0.327 & 0.562 && 0.013 & 0.018 && 0.036 & 0.049 && 0.077 & 0.107 && 0.080 & 0.102 && 0.089 & 0.112 \\
(250,75) && 0.034 & 0.129 && 0.035 & 0.116 && 0.069 & 0.241 && 0.009 & 0.012 && 0.026 & 0.033 && 0.062 & 0.079 && 0.062 & 0.076 && 0.070 & 0.089 \\
(250,100) && 0.005 & 0.033 && 0.008 & 0.032 && 0.014 & 0.064 && 0.008 & 0.011 && 0.024 & 0.030 && 0.057 & 0.075 && 0.059 & 0.074 && 0.076 & 0.095 \\
(500,50) && 0.145 & 0.307 && 0.110 & 0.218 && 0.242 & 0.484 && 0.011 & 0.013 && 0.026 & 0.035 && 0.056 & 0.074 && 0.053 & 0.067 && 0.069 & 0.087 \\
(500,75) && 0.009 & 0.064 && 0.009 & 0.044 && 0.018 & 0.100 && 0.007 & 0.009 && 0.018 & 0.022 && 0.042 & 0.055 && 0.042 & 0.053 && 0.053 & 0.068 \\
(500,100) && 0.002 & 0.003 && 0.004 & 0.005 && 0.007 & 0.009 && 0.006 & 0.008 && 0.017 & 0.021 && 0.041 & 0.054 && 0.040 & 0.053 && 0.056 & 0.070 \\
\bottomrule
\end{tabularx}}
\label{tab:DGP3M_parameters}
\begin{tablenotes}
\tiny
\item
\end{tablenotes}
\end{table}
\section{Application to the U.S. Commercial Banking Sector}
\label{SEC:application}
In this section, we apply the developed method for stochastic cost
frontier model to analyze the cost efficiency of the U.S. large commercial banks in presence of a series of gradual deregulation that allowed banks to increase their capacity of operation. We use the same dataset used by \citet{Fengetal2017}, and focus our analysis on a sample
of banks that operate continuously over the period 1986 to 2005 (thereby
mitigating the impact of entry and exit) with assets of at least \$1
billion in 1986 dollars. Data supporting the findings of this
study are available upon request. As briefly explained in \citet{Fengetal2017}, the banking sector over this period saw a number of gradual deregulation that allowed banks to increase the capacity of operation. In particular, the exact timing of the deregulation varied at the state level, and it was not until June 1997 that banks were allowed to operate across states as a result of the Riegle-Neal Interstate Banking and Branching Efficiency Act of 1994.\footnote{See \citet{Fengetal2017} and \citet{JayaratneStrahan1997} for more
detailed discussion of the history of deregulation in the banking
sector.} Given this context, our method that allows us to group banks based on the time-varying frontiers is well suited to capture the effect
of gradual deregulation, as well as to analyze the inefficiency
of banks in presence of such deregulation.
To set the stage, let $i=1,2,\ldots,N$ denote the banks, $t=1,2,\ldots,T$
denote the time periods. The data is recorded in quarterly frequency,
over 1986 to 2005, with $T=80$ and consists of $N=466$ banks. We
assume that banks use three inputs to generate three outputs. Specifically, the inputs used are: (i) price of labor, $W_{it1}$, (ii) price of purchased funds, $W_{it2}$, and (iii) price of core deposits, $W_{it3}$. Generated outputs are: (i) consumer loans, $Y_{it1}$, (ii) non-consumer loans, $Y_{it2}$, consisting of industrial, commercial, and real estate loans, and (iii) securities, $Y_{it3}$, which includes all non-loan financial assets. Summary statistics of these variables are
reported in Table \ref{tab:Summary_Statistics} in Appendix \ref{APP:Tables}.
We estimate a cost frontier $C(Y_{it}, W_{it})$. Accordingly, $Y_{it}$ are output quantities (loan categories, securities) and $W_{itj}$ are input prices. In particular, $W_{it2}$ is the \emph{price of purchased funds}, not an output; loans are treated as outputs under the intermediation view of banking services. This mapping is consistent with cost duality and with our specification, in which the frontier uses \emph{grouped, smoothly time-varying coefficients} to capture heterogeneous technological regimes under staggered state-level deregulation.
The particular variant of the panel SF model we study is a panel stochastic cost frontier model adapted from \citet{Greene2005}:
\begin{align}
\log c_{it}^{*}=\alpha_{i}^{0} & +u_{i} +\alpha_{i}(\tau_{t})+\beta_{i1}(\tau_{t})\log w_{it1}+\beta_{i2}(\tau_{t})\log w_{it2}\label{eqn:application_model}\nonumber \\
& +\beta_{i3}(\tau_{t})\log y_{it1}+\beta_{i4}(\tau_{t})\log y_{it2}+\beta_{i5}(\tau_{t})\log y_{it3}+v_{it},
\end{align}
where linear homogeneity is imposed in input prices by the normalizations: $c_{it}^{*}=C_{it}/W_{it3}$, $w_{itl}=W_{itl}/W_{it3}$ for $l=1,2$ and $y_{itl}=Y_{itl}/W_{it3}$ for $l=1,2,3$. The inefficiency term $u_{i}\geq0$ enters the model positively as cost frontier models are derived from the dual cost minimization problem of the firm where the cost function is assumed to be Cobb-Douglas.
In a cost frontier, interest-rate conditions primarily operate through $W_{it}$, while local demand affects the output $Y_{it}$.
Our specification already conditions on $(Y_{it},W_{it})$ which are allowed to vary smoothly over time and across latent regimes.
Adding rate or local controls would double-count channels captured by $(Y_{it},W_{it})$.
We estimate the model in \eqref{eqn:application_model} using the
method described in the previous section, setting the tuning parameters
as in Section \ref{SEC:tuning}. We also check the sensitivity of
parameter $c_{\lambda}$ in Step 3 and $\tilde{c}_{\lambda}$ in Step
5 of the proposed method as in Section \ref{SEC:tuning}. Different
values of $c_{\lambda}$ and $\tilde{c}_{\lambda}$ deliver the same
classification results.
\begin{figure}[!htb]
\caption{Scatter plot of elements in $\hat{\vartheta}_{i} = (\hat{\pi}_{i}',\hat{\sigma}_{vi})'$}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Grouped_Frontiers_pihat.pdf}
\label{fig:Grouped_Frontiers_pihat}
\begin{tablenotes}
\vspace{-24pt}
\footnotesize
\item
\emph{Note}: Estimates classified as group 1 are plotted as blue dots, while that for group 2 are plotted as red dots.
\end{tablenotes}
\end{figure}
As in the simulations, we set $\bar{K}=4$. The information criteria
in step 3 selects the optimal number of group for the banks to be
two, splitting $N=466$ banks into $(N_{1},N_{2})=(113,353)$. Figure
\ref{fig:Grouped_Frontiers_pihat} depict the scatter plots of the
elements in $\hat{\vartheta}_{i}=(\hat{\pi}_{i}',\hat{\sigma}_{vi})'$,
that collects the parameters obtained from individual level estimation
in step 1 for classification in step 2. Note $m=\left\lfloor T^{1/5}\right\rfloor =2$,
so each $\hat{\vartheta}_{i}$ is a 12 by 1 vector. Individual estimates classified as groups 1 and 2 are depicted as blue and red dots, respectively. Panel (a) depicts the scatter plot of the estimates $\hat{\pi}_{i1}$ against $\hat{\vartheta}_{i12}$, while panels (b)--(f) depict the respective coefficients on the inputs/outputs $(w_{itl},y_{itl})$ against $(w_{itl}B_{1}(\tau_{t}),y_{itl}B_{1}(\tau_{t}))$. We discuss what drives the classification in Appendix \ref{app:read-groups-lean}.
\begin{figure}[!htb]
\caption{Grouped Frontiers of the U.S. Large Commercial Banks}
\centering
\includegraphics[width=0.8\textwidth]{Figures/Grouped_Frontiers_Application_SE.pdf}
\label{fig:Grouped_Frontiers_Application_SE}
\begin{tablenotes}
\vspace{-24pt}
\footnotesize
\item
\emph{Note}: Top row depict the group 1 time-varying cost frontiers while the bottom row depict that of group 2. Solid lines are the point estimates and the shaded regions are 95\% confidence interval.
\end{tablenotes}
\end{figure}
Figure \ref{fig:Grouped_Frontiers_Application_SE} depicts the frontiers.
The top row depicts the time-varying frontiers of group 1, while the bottom row depicts that of group 2. Solid lines in blue and red are the point estimates for group 1 and group 2 respectively, and the shaded regions depict the 95\% confidence interval. It is evident that there
are substantial time-variations in the estimates, which may be a result
of increasing the capacity of operation as a result of deregulation
in the banking sector. Figure \ref{fig:ES_Application} shows the
estimated economies of scale experienced by two groups of banks, $k=1,2$
defined by the inverse of the sum of elasticities of output, $1/(\hat{\beta}_{(k)3}(\tau_{t})+\hat{\beta}_{(k)4}(\tau_{t})+\hat{\beta}_{(k)5}(\tau_{t}))-1$.
The estimates on economies of scale are comparable to the ones found
in \citet{Greene2005} and suggest some considerable time-variations
for both groups, with group 2 banks enjoying larger economies of scale.
\begin{figure}[ht]
\caption{Estimates of Economy of Scale }
\centering
\includegraphics[width=0.6\textwidth]{Figures/ES_Application.pdf}
\label{fig:ES_Application}
\begin{tablenotes}
\footnotesize
\item
\emph{Note}: Estimates of the economies of scale for each groups are calculated using the point estimates $\hat{\beta}_{(k)l}(\tau_{t})$ for $l=3,4,5$ and $k=1,2$.
\end{tablenotes}
\end{figure}
The results from Step 5 of the proposed method suggest that intercept and idiosyncratic random
effects inefficiency terms, $\alpha^{0}+u$, possess a mixture distribution structure. This result indicates that not only do frontiers form two distinct groups, but so do the level terms that represent the inefficiency of individual banks. The estimated values of the parameters along with the standard errors are presented in Table \ref{tab:Application}. The results
suggest that there are no substantial differences in the standard deviation of random noise, $\hat{\sigma}_{v}$s, although there are significant differences in the standard deviation of the inefficiency terms.
As we briefly mentioned in Section \ref{SEC:inefficiency}, since $\alpha_{i}^{0}$ differs across $i$, we cannot make a valid ranking of the inefficiencies. Luckily, $\hat{\alpha}_{\left(1\right)}^{0}$
and $\hat{\alpha}_{\left(2\right)}^{0}$ are not statistically different (by the likelihood-ratio test) and as such we can view them the same and construct a ranking of the inefficiencies. Recall that we were estimating cost frontiers. From the estimation results, we can view
\begin{equation}
\alpha_{i}^{0}+u_{i}\overset{d}{\sim}\begin{cases}
\begin{array}{c}
\hat{\alpha}^{*}+\left|N\left(0,\hat{\sigma}_{u\left(1\right)}^{2}\right)\right|\\
\hat{\alpha}^{*}+\left|N\left(0,\hat{\sigma}_{u\left(2\right)}^{2}\right)\right|
\end{array} & \begin{array}{c}
\textrm{with probability }\hat{\tau}\\
\textrm{with probability }1-\hat{\tau}
\end{array},\end{cases}\label{eq:mixture_distribution}
\end{equation}
where $\hat{\alpha}^{*}=\hat{\tau}\hat{\alpha}_{(1)}^{0}+\left(1-\hat{\tau}\right)\hat{\alpha}_{(2)}^{0}.$ We compute $\widehat{\textrm{E}\left(\alpha_{i}^{0}+u_{i}|\varepsilon_{i1},...,\varepsilon_{iT}\right)}$ using (\ref{eq:inefficiency_post}).
We compare the ranking in the homogeneous case where the frontiers and
the variances of $v_{it}$ are assumed the same across firms and the
inefficiency term comes from one distribution. The result of top 60 is reported in Figure \ref{fig:Inefficiency_Ranking_Ranking} in Appendix \ref{APP:application_fig}. We can see that the two rankings differ greatly after the top 3. This highlights the importance of classification to ensure valid inference of inefficiency term.
\begin{table}[H]
\centering
\caption{Estimates of $\hat{\sigma}_{v}$s and $\hat{\varrho}$}
\vspace{0.2cm}
\begin{tabularx}{\textwidth}{X X X X X XX }
\toprule
$\hat{\sigma}_{v(1)}$ & $\hat{\sigma}_{v(2)}$ & $\hat{\tau}$ & $\hat\alpha^0_{(1)}$ & $\hat{\sigma}_{u(1)}$ & $\hat\alpha^0_{(2)}$ & $\hat{\sigma}_{u(2)}$ \\
\midrule
0.0862 & 0.0855 & 0.8748 & 0.0157 & 0.4426 & 0.6161 & 0.7756\\
(0.0041) & (0.0008) & (0.1017) & (0.3960) & (0.0362) & (0.1708) & (0.0235) \\
\bottomrule
\end{tabularx}
\label{tab:Application}
\begin{tablenotes}
\footnotesize
\item
\emph{Note}: Reported in parentheses are the standard errors.
\end{tablenotes}
\end{table}
\section{Conclusion}
In this paper, we develop a general framework for panel SF models with latent
group structures. A natural concern is whether allowing for multiple frontiers weakens the interpretation of inefficiency. Our results suggest the opposite. By accounting for latent technological regimes, we prevent unobserved heterogeneity from being mistakenly absorbed into the inefficiency term. Inefficiency in our framework is always measured relative to the appropriate group frontier. This distinction is crucial in empirical applications, such as our U.S. banking study, where ignoring heterogeneity would miscalculate inefficiency.
Two extensions are worth mentioning. First, our framework cannot be
directly generalized to endogenous cases where covariates, $x$ are correlated with the error term, $v$. Extending the framework to accommodate endogeneity is an important direction for future work. Second, it would be valuable to explore a one-step HAC algorithm that avoids the use of information criteria, as proposed by \citet{Mugnier2025}.
\begin{thebibliography}{99}
\harvarditem[Aigner et al.]{Aigner et al.}{1977}{Aigneretal1977}
\textsc{Aigner, D., C. A. K. Lovell, and P. Schmidt} (1977): ``Formulation
and Estimation of Stochastic Frontier Production Function Models,''\ \emph{Journal
of Econometrics,} 6, 21-37.
\harvarditem[Ando and Bai]{Ando and Bai}{2016}{AndoBai2016}
\textsc{Ando, T., and J. Bai} (2016): ``Panel Data Models with Grouped
Factor Structure under Unknown Group Membership,''\ \emph{Journal
of Applied Econometrics,} 31, 163-191.
\harvarditem[Atak et al.]{Atak et al.}{2025}{Ataketal} \textsc{Atak,
A., T. Yang, Y. Zhang, and Q. Zhou} (2025): ``Specification Tests
for Time-Varying Coefficient Panel Data Models,'' \textit{Econometric
Theory}, 41 (1), 123-170.
\harvarditem[Bonhomme and Manresa]{Bonhomme and Manresa}{2015}{BonhommeManresa}\textsc{Bonhomme,
S., and E. Manresa} (2015): ``Grouped Patterns of Heterogeneity in
Panel Data,'' \textit{Econometrica}, 83, 1147-1184.
\harvarditem[Chen]{Chen}{2019}{Chen2019} \textsc{Chen J.}
(2019): ``Estimating Latent Group Structure in Time-Varying Coefficient
Panel Data Models,'' \textit{Econometrics Journal}, 22, 223-240.
\harvarditem[Chen et al.]{Chen et al.}{2014}{Chenetal2014}
\textsc{Chen, Y. Y., P. Schmidt, and H. J. Wang} (2014): ``Consistent Estimation
of the Fixed Effects Stochastic Frontier Model,'' Journal of Econometrics,
181(2) 65-76.
\harvarditem[Cheng et al.]{Cheng et al.}{2024}{Chengetal2024}
\textsc{Cheng, M., S. Wang, L. Xia, and X. Zhang } (2024): ``Testing
Specification of Distribution in Stochastic Frontier Analysis,''\ \emph{Journal
of Econometrics,} 239.
\harvarditem[Colombi et al.]{Colombi et al.}{2014}{Colombi2014}
\textsc{Colombi, R., S. C. Kumbhakar, G. Martini, and G. Vittadini}
(2018): ``Closed-skew Normality in Stochastic Frontiers with Individual
Effects and Long/Short-run Efficiency,''\ \emph{Journal of Productivity
Analysis,} 42, 123-136.
\harvarditem[Dong and Linton]{Dong and Linton}{2018}{DongLinton2018}
\textsc{Dong, C., and O. Linton} (2018): ``Additive Nonparametric
Models with Time Variable and Both Stationary and Nonstationary Regressors,''\ \emph{Journal
of Econometrics,} 207, 212-236.
\harvarditem[Everitt et al.]{Everitt et al.}{2011}{Everittetal}
\textsc{Everitt, B. S., S. Landau, M. Leese, and D. Stahl }(2011).
Cluster Analysis. 5th ed., Wiley, Wiley Series in Probability and
Statistics.
\harvarditem[Feng et al.]{Feng et al.}{2017}{Fengetal2017}
\textsc{Feng, G., J. Gao, B. Peng, and X. Zhang} (2017): ``A Varying-Coefficient
Panel Data Model with Fixed Effects: Theory and an Application to
US commercial Banks,''\ \emph{Journal of Econometrics,} 6, 68-82.
\harvarditem[Galan et al.]{Galan et al.}{2014}{Galan2014}
\textsc{Galán, J. E., and H. Veiga, and M. P. Wiper} (2014): ``Bayesian
Estimation of Inefficiency Heterogeneity in Stochastic Frontier Models,''
\textit{Journal of Productivity Analysis}, 42, 85-101.
\harvarditem[Greene]{Greene}{2005a}{Greene2005a} \textsc{Greene
W.} (2005a): ``Fixed and Random Effects in Stochastic Frontier Models,''
\textit{Journal of Productivity Analysis}, 23, 7-32.
\harvarditem[Greene]{Greene}{2005b}{Greene2005} \textsc{Greene
W.} (2005b): ``Reconsidering Heterogeneity in Panel Data Estimators
of the Stochastic Frontier Model,'' \textit{Journal of Econometrics},
126, 269-303.
\harvarditem[Huang et al.]{Huang et al.}{2020}{HuangEtal2020}
\textsc{Huang, W., S. Jin, and L. Su} (2020): ``Identifying Latent
Grouped Patterns in Cointegrated Panels,'' \textit{Econometric Theory},
36(3), 410-456.
\harvarditem[Jayaratne and Strahan]{Jayaratne and Strahan}{1997}{JayaratneStrahan1997}
\textsc{Jayaratne, J., and P. E. Strahan} (1997): ``The Benefits
of Branching Deregulation,''\ \emph{Economic Policy Review,} 3(4),
13-29.
\harvarditem[Jondrow et al.]{Jondrow et al.}{1982}{Jondrowetal1982} \textsc{Jondrow, J., C. A. K. Lovell, I. M. Materov, and P. Schmidt} (1982): ``On the Estimation of Technical Inefficiency in the Stochastic Frontier Production Function Model,''\ \emph{Journal of Econometrics,} 19(2-3), 233-238.
\harvarditem[Kumbhakar et al.]{Kumbhakar et al.}{2014}{Kumbhakar2014}
\textsc{Kumbhakar, S. C., G. Lien, and J. B. Hardaker} (2014): ``Technical
Efficiency in Competing Panel Data Models: A study of Norwegian Grain
Farming,''\ \emph{Journal of Productivity Analysis,} 41, 321-337.
\harvarditem[Kumbhakar and Lovell]{Kumbhakar and Lovell}{2000}{KumbhakarLovell2000}
\textsc{Kumbhakar, S. C., and C. A. K. Lovell}(2000). Stochastic Frontier
Analysis, Cambridge University Press.
\harvarditem[Kumbhakar et al.]{Kumbhakar et al.}{2022}{Kumbhakaretal2022a}
\textsc{Kumbhakar, S. C., C. Parmeter, and V. Zelenyuk}(2022): ``Stochastic
Frontier Analysis: Foundations and Advances I,'' \textit{Handbook
of Production Economics} ed. by S. C. Ray, R. G. Chambers, and S.
C. Kumbhakar, Springer, Chap. 8, pp. 331-370.
\harvarditem[Lai and Kumbhakar]{Lai and Kumbhakar}{2023}{LaiKumbhakar2023}
\textsc{Lai H. P., and S. C. Kumbhakar} (2023): ``Panel Stochastic
Frontier Model With Endogenous Inputs and Correlated Random Components,''
\textit{Journal of Business \& Economic Statistics}, 41:1, 80-96.
\harvarditem[Lin and Ng]{Lin and Ng}{2012}{LinNg2012} \textsc{Lin,
C.C. and S. Ng} (2012): ``Estimation of Panel Data Models with Parameter
Heterogeneity when Group Membership is Unknown,'' \textit{Journal
of Econometric Methods}, 1(1):42-55.
\harvarditem[Loyo and Boot]{Loyo and Boot}{2024}{LoyoBoot}
\textsc{Loyo, J. A., and T. Boot} (2024): ``Grouped Heterogeneity
in Linear Panel Data Models with Heterogeneous Error Variances,''
\textit{Journal of Business \& Economic Statistics}, 1-13.
\harvarditem[Meeusen and van Den Broeck]{Meeusen and van Den Broeck}{1977}{MeeusenvanDenBroeck1977}
\textsc{Meeusen, W., and J. {van Den Broeck} }(1977): ``Efficiency
Estimation from Cobb-Douglas Production Functions with Composed Error,''
\ \emph{International Economic Review, } 18, 435-444.
\harvarditem[Mugnier]{Mugnier}{2025}{Mugnier2025}\textsc{Mugnier,
M.} (2025): ``A Simple and Computationally Trivial Estimator for Grouped
Fixed Effects Models,'' \textit{Journal of Econometrics}, 250, 106011.
\harvarditem[Park and Simar]{Park and Simar}{1994}{ParkSimar1994}
\textsc{Park, B. U., and L. Simar} (1994): ``Efficient Semiparametric
Estimation in a Stochastic Frontier Model,''\ \emph{Journal of the
American Statistical Association}, 89(427), 929-936.
\harvarditem[Su et al.]{Su et al.}{2016}{SuEtal2016} \textsc{Su,
L., Z. Shi, and P. C. B. Phillips} (2016): ``Identifying Latent Structures
in Panel Data,''\ \emph{Econometrica}, 84(6), 2215-2264.
\harvarditem[Su et al.]{Su et al.}{2019}{SuEtal2019} \textsc{Su,
L., X. Wang, and S. Jin} (2019): ``Sieve Estimation of Time-Varying
Panel Data Models With Latent Structures,''\ \emph{Journal of Business
\& Economic Statistics}, 37(2), 334-349.
\harvarditem[Tsionas and Kumbhakar]{Tsionas and Kumbhakar}{2014}{Tsionas2014}
\textsc{Tsionas, E. G and S. C. Kumbhakar } (2014): ``Firm Heterogeneity,
Persistent and Transient Technical Inefficiency: A Generalized True
Random-Effects Model,''\ \emph{Journal of Applied Econometrics},
29(1), 110--132.
\harvarditem[Tsionas et al.]{Tsionas et al.}{2023}{Tsionas2023}
\textsc{Tsionas, M., F. C. Parmeter, and V. Zelenyuk} (2023): ``Bayesian
artificial neural networks for frontier efficiency analysis,'' \textit{Journal
of Econometrics}, 236(2), 105491.
\harvarditem[Wang and Su]{Wang and Su}{2021}{WangSu2021}
\textsc{Wang, W. and L. Su} (2021): ``Identifying Latent Group Structures
in Nonlinear Panels,'' \textit{Journal of Econometrics}, 220(2),
272-295.
\harvarditem[Yao et al.]{Yao et al.}{2019}{YaoZhangKum2019}
\textsc{Yao F., F. Zhang, and S. C. Kumbhakar} (2019): ``Semiparametric
Smooth Coefficient Stochastic Frontier Model With Panel Data,''\ \emph{Journal
of Business \& Economic Statistics}, 37(3), 556-572.
\end{thebibliography}
\newpage{}