Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
87,750 characters · 22 sections · 48 citation commands
Panel Stochastic Frontier Models with Latent Group Structures
\def\spacingset#1{ {#1}} \spacingset{1}
\\ Hiroshima University \and {Thomas T. Yang}.}\\ Australian National University \and { Xibin Zhang}.}\\ Monash University}
\spacingset{1.8}
Aigneretal1977 and 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 KumbhakarLovell2000 for early developments and Kumbhakaretal2022a and 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 YaoZhangKum2019's that allows for heterogeneous time-varying coefficients:
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 Galan2014. According to Greene2005a, the true underlying frontier may include unmeasured firm-specific characteristics that reflect the technology in use. 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 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 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). 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 LinNg2012, BonhommeManresa, AndoBai2016, or Hierarchical Agglomerative Clustering (HAC) Chen2019. The second approach, introduced by SuEtal2016 performs simultaneous estimation and classification by penalizing the joint mean squared errors. Further developments along this line include SuEtal2019, HuangEtal2020 and 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)). However, since the log-likelihood function is complex and highly nonlinear, applying the method of 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 ParkSimar1994, YaoZhangKum2019 or the error term distribution Greene2005, LaiKumbhakar2023. Additionally, there is a growing literature on specification tests for the distribution of inefficiency 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 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 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)) alluded, we present the hybrid approach in the main text using a random effects specification of the SF model, and show in Appendix (ref) how it extends to other variants, including the fixed effects SF models studied by Greene2005a, Greene2005, Chenetal2014, and 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), 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) examines the theoretical properties of our procedure, specifying the conditions required for tuning parameters. In Section (ref), 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), 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) extends our approach to fixed effects SF models under distributional assumptions. Appendix (ref) further generalizes this to fixed effects SF models without any distributional restrictions, allowing for more flexible conditions on the inefficiency term. Appendix (ref) discusses a relaxation of the normalization condition. Appendix (ref) addresses the possibility that the inefficiency term equals zero with positive probability. Appendix (ref) 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) and (ref) describe, respectively, the approximate likelihood function of the model and the clustering method used. Proofs of the theorems and propositions from Section (ref) are provided in Appendix (ref), with technical lemmas presented in Appendix (ref). Additional simulation and application results are included in Appendix (ref), and Appendix (ref) further discusses the properties of the information matrix in support of Appendix (ref).
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) to (ref). These restrictions are relaxed to the general case subsequently. The technical conditions and main theorems are postponed to the next section.
We illustrate our approach using the panel SF model with random effects (RE) and relegate other variants of the model to Appendix (ref). The particular model we consider modifies YaoZhangKum2019's by incorporating heterogeneous, time-varying coefficients:
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 YaoZhangKum2019, we assume that
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).
The assumption given by ((ref)) 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 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 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 Tsionas2014 and LaiKumbhakar2023. FE-type models are discussed separately in Appendices (ref) and (ref), 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,
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)$:
although we generalize it to allow them to differ in Appendix (ref).\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).} 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)). 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.
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 DongLinton2018 and 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
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
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
The last two lines of equation ((ref)) represent three equivalent ways of expressing the approximation.
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)) and ((ref)) 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.
Using the approximation from ((ref)), 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
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.
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., Chen2019). Details of the HAC method are provided in Appendix (ref) and refer the readers to 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\} $.
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:
where $\lambda_{NT}$ is a suitable penalty term. The optimal number of groups is the minimizer of ((ref)) \[ \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}. \]
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)), 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.
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 Greene2005. However, the latent class model is only a small part of 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)) 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 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.
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:
where $\tilde{f}$ is the likelihood function defined in ((ref)), and we incorporate estimates from Step 3, as detailed in Section (ref).
To determine the optimal number of mixtures, we introduce a new information criterion for this task:
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)) \[ \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). \]
The outline of the estimation procedure is as follows:
We examine classification consistency in Section (ref). Subsequently, we discuss the large sample properties of the post-classification estimators in Section (ref).
Assumption (ref) imposes a condition of weak dependence across $t$, noting that independence across $i$ is not required for classifications. Assumption (ref) requires that $x$ and $v$ have finite $q$-th moment. Assumption (ref) is the classic full rank condition. Assumption (ref) stipulates that the coefficients are $\kappa$-th order differentiable, a standard condition for nonparametric or semiparametric estimation. Assumption (ref) requires that at least one of the coefficients, including the variance of $v$, must differ across groups.
Assumption (ref) 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))); a similar condition was imposed in 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) is satisfied. We do not have the usual bias and variance tradeoff here, as explained below. The results of Theorem (ref) 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) (ii) in the subsequent section. With the aforementioned technical conditions, we demonstrate the consistency of the classification.
Theorem (ref) (i) establishes the uniform convergence of $\hat{\vartheta}_{i}$ for $i=1,2,\ldots,N$ provided $m\to\infty$. Building on this, Theorem (ref) (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}$.
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.
Assumption (ref) further imposes independence across $i.$ Assumption (ref) says the number of members in each group is proportional to $N$. This condition is not necessary, but it facilitates expositions. Assumption (ref) specifies the distributional conditions on the error terms. As explained in Section (ref), these conditions are essential for the identification of $\sigma_{u}^{2}$, and common in the literature, (see, for e.g., YaoZhangKum2019). Assumption (ref) 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). However, condition $\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)), 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). 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(
\right)_{\left(p+1\right)\times\left(m-1+mp\right)}, \]
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)) of Appendix (ref). 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). 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.
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}$).
Proposition (ref) 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.
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) and only present Design 3 here.
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 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.
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). 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)) and ((ref))), which cannot effectively handle mixtures with more than two components. In Appendix (ref), 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.
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).
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) for the benchmark specification with $c_{\lambda}=\tilde{c}_{\lambda}=1$. Additionally, Table (ref) 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) in Appendix (ref) 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) - (ref), 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), obtained by averaging over 500 MC iterations, are reasonable.
Illustrated in Figures (ref), (ref) and (ref) in Appendix (ref) 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) 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) and (ref) 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.
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 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 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 Fengetal2017 and 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) in Appendix (ref).
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 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 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 Greene2005:
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 (ref) using the method described in the previous section, setting the tuning parameters as in Section (ref). 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). Different values of $c_{\lambda}$ and $\tilde{c}_{\lambda}$ deliver the same classification results.
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) 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).
Figure (ref) 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) 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 Greene2005 and suggest some considerable time-variations for both groups, with group 2 banks enjoying larger economies of scale.
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). 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), 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
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)). 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) in Appendix (ref). 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.
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 Mugnier2025.