Estimating the Number of Components in Panel Data Finite Mixture Regression Models with an Application to Production Function Heterogeneity

Estimating the Number of Components in Panel Data Finite Mixture Regression Models with an Application to Production Function Heterogeneity Yu Hao Faculty of Business and Economics The University of Hong Kong [email removed] Hiroyuki Kasahara Vancouver School of Economics The University of British Columbia [email removed]

Introduction

Finite mixture models provide a flexible and natural framework for representing unobserved heterogeneity across a finite number of latent classes. Due to their adaptability, these models have found extensive applications in empirical studies across various disciplines since the seminal introduction of a two-component normal mixture model by Pearson1894. In economics, often utilizing panel data, finite mixture models have been particularly influential in capturing unobserved individual-specific effects in fields such as labor economics, health economics, and industrial organization.\footnote{For instance, Heckman1984 utilize finite mixture models as an alternative method for handling unobserved heterogeneity in unemployment duration analysis. Similarly, Keane1997 and Cameron1998 employ these models to analyze dynamic decision-making regarding schooling and occupational choices in the presence of unobserved heterogeneity in human capital. In health economics, Deb1997 propose a finite mixture negative binomial model to capture the unobserved dispersion in healthcare utilization among elderly populations. Additionally, finite mixture models are extensively used in consumer segmentation within industrial organization, as demonstrated by Kamakura1989 and Andrews2003. Kasahara2009 and Hu2012 establish the conditions under which finite mixture models are non-parametrically identified using panel data.} The theoretical underpinnings and practical applications of finite mixture models have been extensively explored by Titterington1985, Lindsay1995, and McLachlan2000.

A critical aspect in the application of finite mixture models is determining the appropriate number of components, which often corresponds to the number of latent individual types or abilities in economic models. Selecting an incorrect number of components can lead to estimation biases, excessive computational costs, or identification difficulties. Therefore, developing reliable statistical procedures for accurately selecting the number of components in finite mixture models is essential for empirical analysis.

This paper investigates the statistical properties and practical performance of methods for determining the number of components in panel data finite mixture regression models, including the Likelihood Ratio Test (LRT) and the Akaike and Bayesian Information Criteria (AIC and BIC). Panel data finite mixture models offer significant advantages over their cross-sectional counterparts by utilizing repeated observations for each unit, thereby enhancing parameter identification and mitigating common identification issues associated with cross-sectional analyses. However, the existing literature on asymptotic properties of panel data finite mixture regression models is limited, particularly when the component-specific regression errors are normally distributed—a commonly assumed baseline model in regression analyses. This study addresses this gap by rigorously examining and comparing these estimation methods under normally distributed regression errors.

Specifically, this paper develops the LRT and BIC methods for consistently selecting the number of components in various panel data finite mixture regression models and applies these methods to examine production function heterogeneity. We first consider a baseline scenario where regression errors are conditionally independent and normally distributed across periods, conditional on the latent type. Acknowledging the limitations of the normality assumption, we also analyze cases where regression errors follow more flexible, component-specific normal mixtures. Moreover, we extend the analysis to dynamic panel data finite mixture models, incorporating lagged outcomes as covariates under a Markovian assumption, which accommodates regression errors modeled by an autoregressive process. We also consider a nonparametric procedure to estimate a lower bound on the number of components for finite mixture models under conditional independence using the rank test of Kleibergen06 and Kasahara2014.

Testing for the number of components in normal mixture regression models remains challenging, as standard asymptotic regularity conditions fail due to parameter non-identifiability, singular Fisher information, and boundary parameter issues. While numerous studies have addressed the likelihood ratio test (LRT) for finite mixture models ghoshsen85book, chernofflander95jspi, lemdanipons97spl, chenchen01cjstat, chenchen03sinica, cck04jrssb, garel01jspi, garel05jspi, Chen14joe, deriving its asymptotic distribution through Gaussian processes dacunha99as, liushao03as, zhuzhang04jrssb, azais09esaim, these approaches fail in cross-sectional normal regression models due to: (i) infinite Fisher information for testing, (ii) unbounded log-likelihood function, and (iii) linear dependence between density derivatives for mean and variance parameters.\footnote{chenli09as, Chen2012jasa, and kasaharashimotsu15jasa analyze the LRT asymptotics for univariate models, with Kasahara2019 extending to multivariate cases. Additionally, Amengual2022 develop a score-type test for cross-sectional normal mixtures.}

To the best of our knowledge, whether the aforementioned problems (i)--(iii) of the cross-sectional normal mixture still arise in panel data normal finite mixture models or their extensions remains unknown in the literature. Furthermore, no likelihood-based test has yet been developed for testing the null hypothesis of an $M_0$-component model against an alternative $(M_0+1)$-component model for $M_0\geq 1$ in panel normal regression mixture models with conditionally independent errors.

We show that the higher-order degeneracy of problem (iii) disappears in panel normal mixture models with conditional independence and dynamic panel normal mixture models with a Markov assumption, but problems related to (i) and (ii) arise. Consequently, the existing approaches in dacunha99as, liushao03as, zhuzhang04jrssb, and azais09esaim do not directly apply to this class of panel normal mixture models. We impose bounds on the component-specific variance parameters and mixing proportions to address the unboundedness and infinite Fisher information, respectively, and then analyze the asymptotic distribution of the LRT using reparameterization orthogonal to the direction in which the Fisher information matrix is singular. The likelihood ratio of an $(M_0+1)$-component model against the $M_0$-component model is approximated with local quadratic-form expansion with squares and cross-products of the reparameterized parameters. We demonstrate that the asymptotic null distributions of the LRT statistic are characterized by the maximum of $M_0$ random variables. Building on the LRT tests, we propose a sequential hypothesis testing approach for consistently estimating the number of components.

This paper contributes to the literature in several ways. First, we analyze the likelihood ratio test (LRT) for determining the number of components in panel data normal regression finite mixture models, covering both conditional independence and dynamic panel cases with lagged dependent variables. While kasaharashimotsu15jasa and Kasahara2019 study the LRT for cross-sectional univariate and multivariate normal mixture regression models, respectively, we show that the asymptotic distribution of the LRT statistic in panel data differs due to the absence of higher-order dependencies when repeated outcome measurements are available under conditional independence or Markov assumptions.\footnote{From a technical perspective, our panel data finite mixture models constitute a special case of the multivariate normal framework studied by Kasahara2019. The primary distinction lies in our assumption of either conditional independence or a Markov process within each mixture component, which facilitates improved identification of model parameters. We show that this structure yields asymptotic null distributions for the likelihood ratio test that differ from those in Kasahara2019, mainly because our panel data setting does not involve higher-order singularities.}

Second, while it is known that the log-likelihood function of normal mixture models is unbounded Hartigan1985, whether this issue occurs in panel data remains unexplored. We show that the likelihood ratio test (LRT) statistic is also unbounded in panel normal mixture models with conditionally independent errors when the time dimension is finite, though this issue diminishes as the time dimension grows. This unboundedness may lead to excessive rejection rates in the LRT, which we address by imposing explicit bounds on the component-specific variance parameters.

Third, we establish the consistency of the Bayesian Information Criterion (BIC) and the inconsistency of the Akaike Information Criterion (AIC) for selecting the true number of components in panel data finite mixture models. Our consistency result generalizes keribin00sankhya by relaxing its restrictive condition (P2) through higher-order rank conditions and utilizing a generalized form of Le Cam's differentiability in quadratic mean (DQM) framework liushao03as,kasahara2018arXiv.

Fourth, we empirically analyze production technology heterogeneity using panel data from Chilean manufacturing plants. Our findings reveal significant variation in output elasticities of material inputs and the stochastic processes governing factor-augmented technological changes within narrowly defined industries. This contrasts sharply with standard production function estimation methods, which typically impose homogeneous coefficients across plants olley1996dynamics, levinsohn2003estimating, Ackerberg2015. Our results highlight the necessity of explicitly accounting for unobserved plant heterogeneity beyond Hicks-neutral technological differences in empirical analyses of production functions LiSasaki17arxiv, doraszelski2018measuring, Balat19mimeo, Kasahara2022esri.

The identification and estimation of latent group structures in panel data has received attention in recent studies Kasahara2009, LinNg12jem, Bonhomme15ecma, AndoBai16jae, SuShiPhillips16ecma, LuSu17qe. Finite mixture modeling provides a practical, model-based approach to determining unobserved group structures. Choosing the number of groups is often a prerequisite for classifying each individual's group membership. We can estimate the number of groups in panel data regression models by applying our proposed sequential hypothesis testing approach.

This paper closely relates to Kasahara2022esri, which analyzes the nonparametric identification and estimation of finite mixture production models with unobserved heterogeneity, explicitly assuming a known number of mixture components. In contrast, our study emphasizes testing and estimating the number of technology types, focusing specifically on heterogeneity in output elasticities of material inputs. By leveraging first-order condition expressions, our approach enables more flexible finite mixture model specifications than those studied by Kasahara2022esri.\footnote{Additionally, we analyze a production function specification with factor-augmented technological changes, explicitly testing for plant-level heterogeneity in labor-augmented technological processes, a specification not explored by Kasahara2022esri.} However, as our analysis does not recover the entire production function, our framework cannot address heterogeneity in output elasticities of predetermined inputs such as capital, nor heterogeneity in the stochastic processes governing Hicks-neutral technological changes. Thus, our paper complements Kasahara2022esri by addressing distinct yet related aspects of mixture modeling in production function analysis.

The rest of this paper is organized as follows. In Section (ref), we introduce several classes of panel data finite mixture regression models studied in this paper, along with empirical examples. In Section (ref), we discuss the failure of regularity conditions when analyzing panel data finite mixture normal regression models. Section (ref) analyzes the consistency of Maximum Likelihood Estimation (MLE). Section (ref) analyzes the likelihood ratio test (LRT) for testing $H_0: M=1$ against $H_1: M=2$, while Section (ref) considers the LRT for testing $H_0: M=m$ against $H_1: M=m+1$ for $m\geq 2$. Section (ref) develops a sequential hypothesis testing procedure, and Section (ref) analyzes the consistency of the Bayesian Information Criterion (BIC) for selecting the number of components. Section (ref) discusses the estimation of a lower bound for the number of components using the rank test. Section (ref) presents simulation results, while Section (ref) presents empirical analysis.

In what follows, all limits are taken as $n \rightarrow \infty$ unless otherwise stated. Let $:=$ denote "equals by definition." For a $k\times 1$ vector $\boldsymbol{a}$ and a function $f(\boldsymbol{a})$, let $\nabla_{\boldsymbol{a}}f(\boldsymbol{a})$ denote the $k\times 1$ vector of partial derivatives $(\partial/\partial \boldsymbol{a})f(\boldsymbol{a})$, and let $\nabla_{\boldsymbol{a}\boldsymbol{a}^{\top}}f(\boldsymbol{a})$ denote the $k\times k$ matrix of second partial derivatives $(\partial/\partial \boldsymbol{a}\partial \boldsymbol{a}^{\top})f(\boldsymbol{a})$. Let $||\cdot||$ denote the Euclidean norm. We adopt the convention that capitalized letters, such as $\boldsymbol{W}$, represent random variables, whereas their lowercase counterparts, such as $\boldsymbol{w}$, denote evaluation points.

Heteroskedastic finite mixture panel regression model

We consider finite mixture regression models with panel data, where the panel length \( T \geq 2 \) is fixed, and the number of cross-sectional observations \( n \) tends to infinity. Given \( M \geq 2 \) components, we assume that \( y_t \) is conditionally independent over time. The conditional probability density function of \(\{y_t\}_{t=1}^T\) given \(\{\boldsymbol{x}_t\}_{t=1}^T\), where \( y_t \in \mathbb{R} \) and \(\boldsymbol{x}_t \in \mathbb{R}^q\), is expressed as:

where \( f(y_t \mid \boldsymbol{x}_t; \boldsymbol{\theta}_j) \) denotes the density function of the \( j \)-th component for \( y_t \) given \(\boldsymbol{x}_t\), belonging to a parametric class, and \(\alpha_j\) represents the population proportion of the \( j \)-th component. Let \(\boldsymbol{\vartheta}_M := (\boldsymbol{\alpha}^\top, \boldsymbol{\theta}_1^\top, \ldots, \boldsymbol{\theta}_M^\top)^\top \in \Theta_{\boldsymbol{\vartheta}_M}\), where \(\boldsymbol{\alpha}^\top := (\alpha_1, \ldots, \alpha_{M-1})^\top\) and \(\alpha_M = 1 - \sum_{j=1}^{M-1} \alpha_j\). The vector $\boldsymbol x_t$ may include both time-varying and time-invariant regressors.

We are interested in estimating the number of components $M$ when we have correctly specified parametric mixture models given in ((ref)). As a baseline model, we specify $f(\{y_t\}_{t=1}^T \mid \{\boldsymbol{x}_t\}_{t=1}^T;\boldsymbol{\theta}_j)$ using the normal density function as

where $\boldsymbol{\theta}_j := (\mu_j, \sigma_j^2, \boldsymbol{\beta}_j^\top)^\top \in \Theta_{\boldsymbol\theta}$ with $\sigma_j^2 \in \Theta_{\sigma^2} :=\mathbb{R}_{++}$, and $\boldsymbol\beta_j \in \Theta_{\boldsymbol\beta}$ while $\phi(t) := (2\pi)^{-1/2} \exp(-\frac{t^2}{2})$ is the standard normal probability density function. Given that the assumption of a normal density function is restrictive, we also consider more flexible parametric models using a mixture of normal density functions:

where we assume $\mu_{j1}<\mu_{j2}< \cdots<\mu_{j{\cal{K}}}$ for $j=1,...,M$ while the component-specific parameter is given by $\boldsymbol{\theta}_j := (\boldsymbol\tau_j,\boldsymbol\mu_j, \sigma_j^2, \boldsymbol{\beta}_j^\top)^\top \in \Theta_{\boldsymbol\theta}$, $\boldsymbol\tau_j=(\tau_{j1},...,\tau_{j{\cal{K}}-1}) \in\Theta_{\boldsymbol\tau}$, $\tau_{j{\cal{K}}}:=1-\sum_{k=1}^{{\cal{K}}-1}\tau_{jk} $, $\boldsymbol\mu_j=(\mu_{j1},...,\mu_{{\cal{K}}})^\top \in \Theta_{\mu}^{{\cal{K}}}$ , $\sigma_j^2 \in \Theta_{\sigma^2}$, $\boldsymbol\beta_j \in \Theta_{\boldsymbol\beta}$.\footnote{Testing procedures under other parametric classes of mixture models can be developed based on the results of this paper or existing results in the literature, provided the regularity conditions are satisfied zhuzhang04jrssb.} In ((ref)), the mean parameter $\mu_{jk}$ varies across ${\cal{K}}$ components, but neither the variance parameter $\sigma_j$ nor the coefficient $\boldsymbol\beta_j$ varies across ${\cal{K}}$ components. The constant variance parameter assumption prevents another source of unboundedness in this context, while the constant coefficient $\boldsymbol\beta_j$ assumption ensures that the mixture structure arises solely from the mixture distribution of regression error terms.

To simplify our analysis, we assume that ${{\cal{K}}}\geq 2$ is known to the researcher and that all true mixing proportions $\tau_{jk}$ are non-zero.\footnote{Testing the number of components $K$ simultaneously with the number of components $M$ is certainly an important issue but beyond the scope of this paper, and left for future research.}

Here, $M$ and ${\cal{K}}$ serve distinct roles in our model specification. Specifically, $M$ captures the number of latent types representing permanent unobserved heterogeneity across units, whereas increasing ${\cal{K}}$ enhances the flexibility of the i.i.d. error term distributions. Our primary interest in this paper lies in analyzing the permanent unobserved heterogeneity represented by $M$.

In our empirical application, we will analyze plant heterogeneity beyond Hicks-neutral technology differences by estimating the number of latent technology types. The following example provides a finite mixture specification for analyzing plant heterogeneity in input elasticities.

This specification can be extended to a finite mixture variant of the correlated random effects framework Mundlak1978, Chamberlain1984, where the unobserved individual-specific effects have finite discrete support. In this extension, the covariate vector \( \boldsymbol{x}_t \) may include either the individual-specific averages of the time-varying regressors or even the complete sequence of regressors \(\{\boldsymbol x_t\}_{t=1}^T\) over the entire observation period.

The strict exogeneity assumption on \(\boldsymbol{x}_t\) can be relaxed by adopting a sequential exogeneity framework, where \(\boldsymbol{x}_t\) includes covariates observed up to period \(t-1\). Under sequential exogeneity, coefficients become period-specific, as the dimension of \(\boldsymbol{x}_t\) may vary over time. Although this approach has the advantage of relaxing strict exogeneity, it also substantially increases the number of parameters to estimate, particularly as the number of mixture components grows. To maintain analytical tractability, we focus below on the simpler scenario where the dimension of \(\boldsymbol{x}_t\) is constant and the associated coefficients are time-invariant but extending our theoretical results to accommodate sequential exogeneity is conceptually straightforward.

The assumption of conditional independence of unobserved idiosyncratic shocks over time in the model ((ref)) may be seen as restrictive. To address this limitation, we extend our analysis to finite mixture dynamic panel data models that incorporate lagged outcomes as covariates by considering models where \( \{y_t\}_{t=1}^T \) follows a first-order Markov process, conditional on both the latent component type and the sequence of covariates \( \{\boldsymbol{x}_t\}_{t=1}^T \). Specifically, for the baseline dynamic panel data model, the conditional probability density of \( \{y_t\}_{t=1}^T \) given \( \{\boldsymbol{x}_t\}_{t=1}^T \) is given by ((ref)) with

with

where $\boldsymbol\theta_j=(\mu_{1,j},\mu_j,\sigma_{1,j},\sigma_{j},\boldsymbol\beta_{1,j}^{\top},\boldsymbol\beta_j^{\top},\rho_j)^{\top}$. In an extended case, the densities within each component are flexibly modeled as mixtures of normal densities:

with $\boldsymbol\theta_j=(\boldsymbol\tau_j,\boldsymbol\mu_{1,j}^{\top},\boldsymbol\mu_j^{\top},\sigma_{1,j}^2,\sigma_{j}^2,\boldsymbol\beta_{1,j}^{\top},\boldsymbol\beta_j^{\top},\rho_j)^{\top}$, where $\boldsymbol\tau_j=(\tau_{j1},...,\tau_{{j\cal{K}}-1})^{\top}$ and $\tau_{j{\cal{K}}} = 1 - \sum_{k=1}^{{\cal{K}}-1} \tau_{jk}$, $\boldsymbol\mu_{1,j}=(\mu_{1,j1},...,\mu_{1,j\cal{K}})^{\top}$, and $\boldsymbol\mu_j=(\mu_{j1},...,\mu_{{j\cal{K}}})^{\top}$. The following example illustrates that the finite mixture dynamic panel data model incorporates a scenario in which the regression errors follow a first-order autoregressive (AR(1)) process.

Factor augmented technological changes are also a popular way to incorporate heterogeneity in production functions beyond the Hicks-neutral technological change doraszelski2018measuring,Zhang2019,Raval2019. The heterogeneity in the stochastic process of factor augmented technological changes can be modelled using finite mixture dynamic panel data model as the following example illustrates.

The number of components, denoted by $M_0$, is defined as the smallest integer $M$ such that the data density of $\boldsymbol{w}$ admits the representation ((ref)). Consider a random sample of $n$ with a panel length of $T$ independent observations $\{\boldsymbol{W}_{i}\}_{i=1}^n$, where $\boldsymbol{W}_i = \{ (Y_{it},\boldsymbol{X}^{\top}_{it})\}_{t=1}^T$ is drawn from a true $M_0$-component density $g_M(\{y_{t}\}_{t=1}^T|\{\boldsymbol{x}_{t}\}_{t=1}^T;\boldsymbol\vartheta_{M_0}^*)$ in Equation ((ref)), where $\boldsymbol\vartheta_{M_0}^*$ represents the true parameter value. Because component distributions can be identified only up to permutation, we assume that $\mu_1^*<\mu_2^*<\cdots <\mu_{M_0}^*$ for identification.\footnote{More generally, we may consider a lexicographical order: $\boldsymbol\theta_1^*<\boldsymbol\theta_2^*<\cdots<\boldsymbol\theta_{M_0}^*$.}

In Examples (ref)–(ref), the number of components corresponds to the number of latent technology types, reflecting plant-level heterogeneity in output elasticities of material input or factor-augmented technology stochastic processes. We propose a procedure based on the Likelihood Ratio Test (LRT) for the finite mixture model specified in equation ((ref)) to evaluate the null hypothesis of $H_0: M=1$ against $H_1: M=2$. This provides a systematic method to assess whether significant heterogeneity exists across plants in terms of input elasticities or factor-augmented technologies. Moreover, we further develop the LRT for testing \[ \text{$H_0: M=m$ \ against \ $H_1: M=m+1$} \] where $m$ is some known integer. By sequentially testing $H_0: M=m$ against $H_1: M=m+1$ for $m=1,2,...$, this procedure consistently estimates the number of distinct plant types.

Additionally, we consider determining the number of components using the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), with BIC consistency analyzed in Section (ref). We further propose estimating a lower bound on the number of components based on the nonparametric rank test approach of Kasahara2014, thus avoiding parametric assumptions. These methods are empirically applied to analyze the plant heterogeneity in production function using Examples (ref)–(ref) in Section (ref).

For notational brevity, let $\boldsymbol w = \{y_{t},\boldsymbol{x}_{t}\}_{t=1}^T$ and write $g_M(\boldsymbol w;\boldsymbol\vartheta_M):=g_M(\{y_{t}\}_{t=1}^T|\{\boldsymbol{x}_{t}\}_{t=1}^T;\boldsymbol\vartheta_{M})$ in ((ref)), and let $f(\boldsymbol w;\boldsymbol\theta_j):=f( \{y_t\}_{t=1}^T| \{\boldsymbol{x}_t\}_{t=1}^T;\boldsymbol\theta_j)$ be component-specific density function so that ((ref)) is written as $g_M(\boldsymbol w;\boldsymbol\vartheta_M) = \sum_{j=1}^M \alpha_j f(\boldsymbol w;\boldsymbol\theta_j).$

The failure of the regularity conditions

We consider a random sample of $n$ observations $\{\boldsymbol{W}_{i}\}_{i=1}^n$, where $\boldsymbol{W}_i = \{ (Y_{it},\boldsymbol{X}^\top_{it})^\top \}_{t=1}^T$ from an $M_0$-component density $g_{M_0}(\boldsymbol{w}; \boldsymbol{\vartheta}_{M_0})$ defined in equation ((ref)):

where $\boldsymbol{\vartheta}_{M_0}^* = (\boldsymbol{\theta}_1^*,\boldsymbol{\theta}_2^*,\ldots,\boldsymbol{\theta}_{M_0}^*,\alpha_1^*,\ldots,\alpha_{M_0 -1}^*) \in \Theta_{\boldsymbol{\vartheta}_{M_0}}$ and $\alpha_{M_0}^*=1-\sum_{j=1}^{M-1}\alpha_j^*$. We assume $\mu_{1}^* < \mu_{2}^*, \ldots, < \mu_{M_0}^*$ in the true parameters for identification.

To examine the failure of the regularity conditions of the LRTS in panel finite mixture models, consider testing the null hypothesis $H_0: M=1$ against the alternative hypothesis $H_1: M=2$, where $\boldsymbol{W}_i = \{ (Y_{it}, \boldsymbol{X}_{it}^{\top})^{\top} \}_{t=1}^T$, drawn from a true one-component density $f(\boldsymbol{w}; \boldsymbol{\theta}^*)=\prod_{t=1}^Tf(y_t|\boldsymbol x_t;\boldsymbol\theta^*)$ or $f_1(y_1|\boldsymbol x_1;\boldsymbol\theta^*)\prod_{t=2}^Tf(y_t|y_{t-1},\boldsymbol x_t;\boldsymbol\theta^*)$. The two-component model $g_2(\boldsymbol{w}; \boldsymbol{\vartheta}_2) = \alpha f(\boldsymbol{w}; \boldsymbol{\theta}_1 ) + (1 - \alpha) f(\boldsymbol{w}; \boldsymbol{\theta}_2)$ can generate the true one-component density in two cases: (1) $\boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^*$ and (2) $\alpha = 0$ or $1$. Consequently, $H_0: M=1$ can be partitioned into two sub-hypotheses: $H_{01}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2$ and $H_{02}: \alpha (1-\alpha) = 0$. The regularity conditions of the LRTS fail in any finite mixture models: under $H_{01}$, $\alpha$ is not identified, and the Fisher information matrix for the other parameters becomes singular; under $H_{02}$, $\alpha$ is on the boundary of the parameter space, and either $\boldsymbol{\theta}_1$ or $\boldsymbol{\theta}_2$ is not identified.

Analyzing the asymptotic distribution of the LRTS for the cross-sectional normal mixture, i.e., with $T=1$, is even more challenging because chenli09as: (i) the Fisher information for testing $H_{02}$ is not finite, (ii) the log-likelihood function is unbounded Hartigan1985, and (iii) the first-order derivative of $g_2(\boldsymbol{w}; \boldsymbol{\vartheta}_2)$ with respect to $\sigma_j^2$ is linearly dependent on its second-order derivative with respect to $\mu_j$. We now examine whether problems (i)--(iii) are present in panel normal mixture models with $T\geq 2$.

Regarding problem (i), the issue of the infinite Fisher information for testing $H_{02}$ arises in the panel normal mixture model. The score for testing $H_{02}: \alpha=0$ takes the form \[ \left.\frac{\partial g_2(\boldsymbol{W} ; \alpha,\boldsymbol\theta_1,\boldsymbol\theta_2)}{\partial \alpha} \right|_{\alpha=0,\boldsymbol\theta_2=\boldsymbol\theta^*}= \frac{f(\boldsymbol{W};\boldsymbol\theta_1)}{f(\boldsymbol{W};\boldsymbol\theta^* )}-1, \] where $f(\boldsymbol{W};\boldsymbol\theta)= \prod_{t=1}^T \phi((Y_t-\mu-\boldsymbol X_t^{\top}\boldsymbol\beta)/\sigma)/\sigma$. Then, $\mathbb{E}[\{f(\boldsymbol{W};\boldsymbol\theta_1)/f(\boldsymbol{W};\boldsymbol\theta^* )-1\}^2]=\infty$ when $\sigma_1^2> 2\sigma^{*2}$. This infinite Fisher information also arises for normal mixture density ((ref)) in place of ((ref)) and/or for dynamic panel models ((ref)). For more details, please refer to Proposition (ref). Because the infinite Fisher information causes difficulty in deriving the asymptotic distribution under $H_{02}$, throughout this paper, we focus on testing $H_{01}$ by restricting the value of $\alpha$ to be away from $0$ and $1$ by assuming that $\alpha\in[c_{1}, 1-c_{1}]$ for some positive constant $c_{1}$. Appendix (ref) discusses the asymptotic distribution under $H_{02}$ under the constraint $\sigma_1^2\leq 2\sigma^{*2}-c$ for some $c>0$ based on the framework of andrews01em.

Because we focus on $H_{01}$, our test may not have power against the local alternatives with $\alpha_n\rightarrow 0$ as discussed in Appendix (ref).

Related to problem (ii), the LRTS in normal mixture models with panel data becomes unbounded as the sample size $n$ goes to $\infty$. Define the likelihood ratio statistic with respect to the true parameter under $H_0$ as $ LR_n^*(\boldsymbol{\vartheta}_{2}) := 2\left\{ \sum_{i=1}^n \log g_2(\boldsymbol{W}_i;\boldsymbol{\vartheta}_{2}) - \sum_{i=1}^n\log f(\boldsymbol{W}_i;\boldsymbol{\theta}^*)\right\}, $ where $g_2$ is the density of the two-component finite mixture distribution in ((ref)) with $M=2$ and $\boldsymbol{\theta}^*$ is the true parameter value under $H_0$.

The assumption of non-compactness for the parameter space of $\sigma_j^2$ is crucial for Proposition (ref). Without allowing parameters to approach boundary points (i.e., if compactness were imposed), one could not demonstrate the divergence of the likelihood ratio statistic, thus highlighting the source of unboundedness.

For problem (iii), consider the two-component panel mixture model without covariates given by

where $\boldsymbol\theta_j=(\mu_j,\sigma_j^2)^{\top}$. Then, when evaluated under the null hypothesis $H_{01}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^*$ with $\alpha(1 - \alpha) \neq 0$, the first-order derivative of $g_2(\{y_t\}_{t=1}^T;\alpha,\boldsymbol\theta_1,\boldsymbol\theta_2)$ with respect to $\sigma_j^2$ is not linearly dependent on its second-order derivative with respect to $\mu_j$.

Consequently, the panel mixture model ((ref)) with component density function ((ref)) is strongly identifiable, and the best rate of convergence for estimating the mixing distribution is $n^{-1/4}$ when the number of components is unknown chen95as; see also our analysis under local alternatives in Appendix (ref). The key assumption underlying this result is that the component density function can be factored into a product of density functions, reflecting limited dependence in regression errors across periods under either the conditional independence or Markov assumption. In contrast, the strong identifiability does not hold for the cross-sectional normal mixture or the multivariate normal mixture with general dependence, and its convergence rate becomes as slow as $n^{-1/8}$ when the number of components is over-specified kasaharashimotsu15jasa,Kasahara2019.

As in any finite mixture models, however, the standard asymptotic analysis breaks down in testing $H_{01}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^*$ because $\alpha$ is not identified under $H_{01}$; in addition, the first-order derivative at the true value $\boldsymbol\vartheta^*_2 = (\alpha, (\boldsymbol\theta^*)^{\top}, (\boldsymbol\theta^*)^{\top})^{\top}$ is linear dependent as

Consequently, we derive the asymptotic distribution of the LRTS using a fourth-order Taylor series approximation for the log-likelihood function.

Consistency of the MLE

To address the issue of unboundedness, following hathaway85as, we restrict the parameter space to ensure that the ratio of each component-specific variance to the unconditional variance is bounded away from zero, i.e., $\min_{j} \left( \frac{\sigma_j^2}{\sum_{k=1}^M \alpha_k \sigma_k^2} \right) \geq c_{2},$ for some constant \( c_{2} > 0 \). This constraint enforces a positive lower bound on the component variances when the unconditional variance is strictly positive with $\sum_{k=1}^M \alpha_k \sigma_k^2\geq c_3>0$. We also assume that $\alpha_j$ is bounded away from $0$ and $1$ by a small constant $c_1>0$ so that $\boldsymbol\alpha\in\Theta_{\boldsymbol\alpha,c_1}:=\{\boldsymbol\alpha: \alpha_j\in [c_1,1-c_1], \sum_{j=1}^{M-1}\alpha_j\leq 1-c_1\}$. This assumption prevents issues related to infinite Fisher information, allowing us to focus on $H_{01}$. Let $\boldsymbol c=(c_{1},c_{2},c_3)^{\top}$, and define a restricted parameter space for $\boldsymbol\vartheta_M$ for a model with the component-specific density function ((ref)) by \[ {\bar \Theta}_{\boldsymbol\vartheta_M}(\boldsymbol c) := \{ \boldsymbol\vartheta_M \in \Theta_{\boldsymbol\vartheta_M} : \boldsymbol \alpha\in \Theta_{\boldsymbol\alpha,c_1}, \ \min_j \sigma_j^2\geq c_2 \sum_{k=1}^M \alpha_k\sigma_k^2,\ \sum_{k=1}^M \alpha_k \sigma_k^2\geq c_3 \}. \] For the dynamic panel model ((ref)), we define a restricted parameter space for $\boldsymbol\vartheta_M$ analogously, but with additional restrictions on $\sigma_{1,j}$: $\min_{j} \sigma_{1,j}^2 \geq c_{2} \sum_{k=1}^M \alpha_k \sigma_{1,k}^2 $ and $\sum_{k=1}^M \alpha_k \sigma_{1,k}^2\geq c_3$.

We consider the MLE of $\boldsymbol\vartheta_M$ over a restricted parameter space ${\bar \Theta}_{\boldsymbol\vartheta_M}(\boldsymbol c)$ for the M-component model:

where

is the the log likelihood function of ${\boldsymbol{\vartheta}}_M$. The MLE $\widehat{\boldsymbol\vartheta}_M$ is potentially one of many possible global maximum points of $\ell_n^M({\boldsymbol{\vartheta}}_M)$ when $M$ is greater than the true number of components.

For $M>M_0$, define a set of parameter values for the $M$-component density that generates the true $M_0$-component density by \[ \Theta_{\boldsymbol\vartheta_M}^*:=\left\{\boldsymbol\vartheta_M\in \Theta_{\boldsymbol\vartheta_M}: g_M(\boldsymbol w;\boldsymbol\vartheta_M)=g_{M_0}(\boldsymbol w;\boldsymbol\vartheta_{M_0}^*)\quad\text{for all $\boldsymbol w\in \boldsymbol{\mathcal{W}}$}\right\}. \] For example, when $M_0=1$ with the true density function $f(\boldsymbol w;\boldsymbol\theta^*)$, we have \[ \Theta_{\boldsymbol\vartheta_2}^*:= \left\{ (\alpha,\boldsymbol{\theta}_1,\boldsymbol{\theta}_2) \in \Theta_{\boldsymbol\vartheta_2}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^*; \alpha=1 \text{ and } \boldsymbol\theta_1=\boldsymbol\theta^*; \alpha=0 \text{ and } \boldsymbol\theta_2=\boldsymbol\theta^* \right\},\] where $\boldsymbol{\theta}_1$ and $\boldsymbol{\theta}_2$ are component-specific parameters.

The following proposition establishes the consistency of the MLE for over-specified mixture models.

Consequently, Proposition (ref) suggests that the MLE $\widehat{\boldsymbol{\vartheta}}_M$ converges almost surely to a set of parameters for which the true density function $g_{M_0}(\boldsymbol{w}; \boldsymbol\theta_{M_0}^*)$ emerges within the space of overspecified $M$-component density functions.

Likelihood ratio test for $H_0: M=1$ against $H_1: M=2$

We first consider testing $H_0: M=1$ against $H_1: M=2$. Let $\widehat{\boldsymbol{\theta}}_0$ be the one-component MLE that maximizes the one-component log likelihood function $\ell_n^1(\boldsymbol\theta):= \sum_{i=1}^n \log f(\boldsymbol{W}_i; \boldsymbol{\theta}) $. Define the LRTS of testing $H_{01}$ as

where $\hat{\boldsymbol\vartheta}_2$ is the MLE under a restricted parameter space as defined by ((ref)).

To derive the asymptotic distribution of the LRTS, we consider the following one-to-one reparameterization of $\boldsymbol{\theta}_1$ and $\boldsymbol{\theta}_2$ given $\alpha$ Kasahara2012:

where $\boldsymbol{\nu}$ and $\boldsymbol{\lambda}$ are both $q \times 1$ reparameterized parameter vectors.

For the model ((ref)) with density ((ref)), we have $\boldsymbol{\nu}=(\nu_\mu,\nu_\sigma,\boldsymbol{\nu}_{\boldsymbol\beta}^{\top})^{\top}$ and $\boldsymbol\lambda = (\lambda_{\mu},\lambda_{\sigma}, \boldsymbol\lambda_{\beta}^{\top})^{\top}= ( \mu_1 - \mu_2, \sigma_1^2 - \sigma_2^2, (\boldsymbol \beta_1 - \boldsymbol\beta_2)^{\top})^{\top}$ while for the model ((ref)) with ((ref)), $\boldsymbol{\nu}$ and $\boldsymbol\lambda$ are given as $\boldsymbol{\nu}=(\boldsymbol\nu_{\boldsymbol\tau}^{\top},\boldsymbol\mu^{\top},\nu_\sigma,\boldsymbol{\nu}_{\boldsymbol\beta}^{\top})^{\top}$ and $\boldsymbol\lambda = (\boldsymbol\lambda_{\boldsymbol\tau}^{\top},\boldsymbol\lambda_{\boldsymbol\mu}^{\top},\lambda_{\sigma}, \boldsymbol\lambda_{\boldsymbol\beta}^{\top})^{\top}= ( \boldsymbol\tau_1^{\top}-\boldsymbol\tau_2^{\top},\boldsymbol\mu_1^{\top} - \boldsymbol\mu_2^{\top}, \sigma_1^2 - \sigma_2^2, \boldsymbol \beta_1^{\top} - \boldsymbol\beta_2^{\top})^{\top}$. We may also define $\boldsymbol\nu$ and $\boldsymbol\lambda$ appropriately fort the dynamic models ((ref)) with ((ref)) or ((ref)). Here, the parameter $\boldsymbol{\lambda}$ captures a deviation from the one-component model, such that $\boldsymbol \nu=\boldsymbol\theta^*$ and $\boldsymbol\lambda=\boldsymbol 0$ under $H_0: M=1$.

Define the elements of $\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$ as $\boldsymbol{\theta} =(\theta_1,\theta_2,...,\theta_{q})^{\top}$ and $\boldsymbol{\lambda}=(\lambda_1,\lambda_2,...,\lambda_{q})^{\top}$.

Taking the derivative with respect to $\boldsymbol\lambda$ reveals that $\nabla_{\boldsymbol{\lambda}} \log g_2(\boldsymbol{w};\alpha,\boldsymbol{\nu} + ( 1- \alpha) \boldsymbol{\lambda},\boldsymbol{\nu} - \alpha \boldsymbol{\lambda})|_{\boldsymbol\lambda=\boldsymbol 0} = \boldsymbol{0}$. Therefore, the Fisher information matrix is singular, and the standard quadratic approximation fails. Consequently, we use the second-order derivative with respect to $\boldsymbol{\lambda}$ to identify $\boldsymbol{\lambda}$.

Let $f^*$ and $\nabla f^*$ denote $f(\boldsymbol{W}; \boldsymbol{\theta}^*)$ and $\nabla f(\boldsymbol{W}; \boldsymbol{\theta}^*)$, respectively. Define the vector $\boldsymbol{s}(\boldsymbol{W})$ as

where $\widetilde{\text{vech}}_2\left({ {\nabla}_{\boldsymbol\theta\boldsymbol\theta^{\top}} f^*}/{f^*}\right)$ is a $q(q+1)/2\times 1$ vector that contains all unique elements of $q$-dimensional symmetric array ${ {\nabla}_{\boldsymbol\theta\boldsymbol\theta^{\top}} f^*}/{f^*}$, multiply by its frequency, and stacks them into a vector. The function $\boldsymbol{s}_{\boldsymbol{\lambda} \boldsymbol\lambda}(\boldsymbol w)$ essentially comprises the second-order derivatives of the log-likelihood function with respect to $\boldsymbol\lambda$, serving as a score function for identifying $\boldsymbol{\lambda}$. We refer to $\boldsymbol{s}(\boldsymbol{w})$ as a score function. An explicit expression for the score function $\boldsymbol{s}(\boldsymbol{w})$ can be derived using Hermite polynomials, as shown in Appendix (ref), from which the non-singularity of $\boldsymbol{\mathcal{I}}:=\mathbb{E}[\boldsymbol s(\boldsymbol W)\boldsymbol s(\boldsymbol W)^{\top}]$ follows.

Partition $\boldsymbol{\mathcal{I}}$ and define

where $\boldsymbol{G}_{\boldsymbol\lambda,\boldsymbol\nu} \sim \mathcal{N}(\boldsymbol 0,(\boldsymbol{\mathcal{I}}_{\boldsymbol\lambda,\boldsymbol\nu})^{-1})$.

Define $\widehat{\boldsymbol{t}}_{\lambda}$ as

such that $\widehat{\boldsymbol{t}}_{\boldsymbol\lambda}$ is a projection of a Gaussian random variable $\boldsymbol{G}_{\boldsymbol{\lambda}}$ on a cone $\Lambda_{\boldsymbol\lambda}: = \Big\{ v(\boldsymbol{\lambda}) : \boldsymbol{\lambda} \in \mathbb{R}^q\Big\}$, where $v(\boldsymbol{\lambda})$ is a vector of unique elements of $\boldsymbol\lambda\boldsymbol\lambda^{\top}$, where the diagonal elements are divided by 2 :

The following proposition establishes the asymptotic distribution of the LRTS under the null hypothesis $H_0: M=1$.

For Assumption (ref)(c), the non-singularity of $\boldsymbol{\mathcal{I}}$ is proved in Lemma (ref) for the model with component-specific density function ((ref)) with ((ref)). For other models, while analytically proving the non-singularity of $\boldsymbol{\mathcal{I}}$ is not straightforward, each element of the score function $\boldsymbol s(\boldsymbol W)$ is a linear transformation of products of posterior probabilities and Hermite polynomials of different degrees. Therefore, we expect $\boldsymbol{\mathcal{I}}=\mathbb{E}[\boldsymbol s(\boldsymbol W) \boldsymbol s(\boldsymbol W)^\top]$ to be finite and non-singular.

Likelihood ratio test for $H_0: M=m$ against $H_1: M=m+1$

If the data is generated from the $M_0$-component model ((ref)), the $(M_0 + 1)$-component model

gives rise to the true density ((ref)) in two cases: (i) two components have identical parameters or (ii) one component has zero mixing proportion. Accordingly, we partition the null hypothesis $H_0: M=M_0$ into two as $H_0 = H_{01}\cup H_{02}$, with $H_{01}: \boldsymbol{\theta}_h = \boldsymbol{\theta}_{h+1} = \boldsymbol{\theta}_{h}^*$ for some $h=1,\ldots,M_0$ and $\alpha_j>0$ for all $j$, and $H_{02}: \alpha_h = 0$ for some $h=1,\ldots,M_0+1$.

Partition $H_{02}$ as $H_{02}=\cup_{h=1}^{M_0}H_{0,2h}$, where $H_{0,2h}: \alpha_h=0$. Define the subset of $\Theta_{\vartheta_{M_0+1}}^*$ corresponding to $H_{0,2h}$ as $\Theta_{\boldsymbol \vartheta_{M_0+1},2h}^* = \{ \boldsymbol \vartheta_{M_0+1} \in \Theta_{\boldsymbol \vartheta_{M_0+1}} : \alpha_h = 0; (\alpha_j, \boldsymbol\theta_j) = (\alpha_j^*, \boldsymbol\theta_j^*)\text{ for } j < h; (\alpha_j, \boldsymbol\theta_j)= (\alpha_{j-1}^*, \boldsymbol\theta_{j-1}^*)\text{ for } j > h\}$. Because $\boldsymbol\theta_h$ is not identified when $\alpha_h = 0$, we take the supremum of the variance of $\nabla_{\alpha_h} \log g_{M_0 + 1}(\boldsymbol W_i; \boldsymbol \vartheta_{M_0+1})$ over $\boldsymbol \vartheta_{M_0+1} \in \Theta_{\boldsymbol \vartheta_{M_0+1},2h}^*$ to examine the finiteness of the Fisher information matrix for testing $H_{0,2h}: \alpha_h = 0$. The Fisher information is infinite unless there is an a priori restriction on the values of component-specific variance, $\sigma_j^2$.

Because the restriction on the values of $\sigma_j^2$ in Proposition (ref) is difficult to justify and not easy to enforce in practice, we focus on testing $H_{01}$.

Partition $H_{01}$ as $H_{01}=\cup_{h=1}^{M_0}H_{0,1h}$, where $H_{0,1h}: \boldsymbol\theta_{h}=\boldsymbol\theta_{h+1}$ with $\mu_1<\cdots < \mu_h=\mu_{h+1}<\cdots < \mu_{M_0+1}$. We impose these inequality constraints on $\mu_j$ for component identification. There are $M_0$ ways to describe the $M_0$ component null model in the space of $(M_{0}+1)$ component models, and each way corresponds to the null hypothesis of $H_{0,1h}: \boldsymbol\theta_{h}=\boldsymbol\theta_{h+1}$ for $h=1,2,..., M_0$. Testing $H_{0,1h}: \boldsymbol\theta_{h}=\boldsymbol\theta_{h+1}$ in the $M_0$-component null models is similar to testing $H_{01}: \boldsymbol\theta_1=\boldsymbol\theta_2$ in the one-component null model in Section (ref).

Define the subset of $\Theta_{\vartheta_{M_0+1}}$ corresponding to $H_{0,1h}$ as

for $h = 1,\ldots,M_0$. The set $\Theta_{\boldsymbol \vartheta_{M_0+1},1}^* := \cup_{h=1}^{M_0} \Theta_{\boldsymbol \vartheta_{M_0+1},1h}^*$ corresponds to $H_{01}=\cup_{h=1}^{M_0}H_{0,1h}$.

We first develop a test statistic for testing test $H_{0,1h}$ by restricting the estimators under the $(M_0 + 1)$-component model to be in a neighbourhood of $\Theta_{\boldsymbol \vartheta_{M_0+1},1}^*$.

Recall that $\mu_1^{*} < \mu_2^* \ldots < \mu_{M_0}^*$. Let $ \underline{\Theta}_{\mu}$ and $\overline{\Theta}_{\mu}$ denote the lower and upper bounds of $\Theta_{\mu}$, respectively. Define

Then, $\{\Theta_{\boldsymbol\theta,h}^*\}_{h=1}^{M_0}$ is a partition of $\Theta_{\boldsymbol\theta}$ such that $\Theta_{\boldsymbol\theta,h}^*$ is a neighbourhood containing $\boldsymbol\theta_h^{*}$ but not $\boldsymbol\theta_j^{*}$ for $j \neq h$. For $h=1,\ldots M_0$, define a restricted parameter space $\boldsymbol{\Psi}_h^* \subset \bar\Theta_{\boldsymbol\vartheta_{M_0 + 1}}(\boldsymbol c)$ as

Note that $\boldsymbol{\Psi}_h^* \cap \Theta_{\boldsymbol \vartheta_{M_0+1},1h}^* \neq \emptyset$ and $\boldsymbol{\Psi}_h^* \cap \Theta_{\boldsymbol \vartheta_{M_0+1},1l}^* = \emptyset$ if $h \neq l$, and $\cup_{h=1}^{M_0}\boldsymbol{\Psi}_h^*= \bar\Theta_{\boldsymbol\vartheta_{M_0 + 1}}(\boldsymbol c)$.

Let $\widehat{\boldsymbol{\Psi}}_h^* $ and $\widehat{\Theta}_{\boldsymbol \theta,h}^* $ be consistent estimators of ${\boldsymbol{\Psi}}_h^* $ and ${\Theta}^*_{\boldsymbol\theta,h}$, which can be constructed from a consistent estimator of $\boldsymbol{\vartheta}_{M_0}^*$ in the $M_0$-component model. We test $H_{0,1h}: \boldsymbol\theta_{h}=\boldsymbol\theta_{h+1}$ by estimating the $(M_0+1)$-component model under the restriction that $\boldsymbol\vartheta^{M_0+1}\in \widehat{\boldsymbol{\Psi}}_h^*$.

For $h=1,2,..., M_0$, define the local MLE that maximizes the log-likelihood function of the $(M_0+1)$-component model under the constraint that $\boldsymbol\vartheta_{M_0+1}\in \widehat{\boldsymbol\Psi}_h^*$ in ((ref)) by \[ \ell^{M_0+1}_n (\widehat{\boldsymbol{\vartheta}}_{M_0 + 1}^h) = \arg\sup_{\boldsymbol{\vartheta}_{M_0 + 1} \in \widehat{\boldsymbol\Psi}_h^* } \ell^{M_0+1}_n (\boldsymbol{\vartheta}_{M_0+1}). \]

Under $H_{0}: M=M_0$, $\boldsymbol\Psi_h^*$ contains a set of parameters $ \Theta_{\boldsymbol \vartheta_{M_0+1},1h}^*$ defined in ((ref)) such that $g_{M_0 + 1}(\boldsymbol{w}; \boldsymbol{\vartheta}_{M_0 + 1})$ is equal to the true density $g_{M_0}(\boldsymbol{w}; \boldsymbol{\vartheta}_{M_0 }^*)$ for any $\boldsymbol{\vartheta}_{M_0 +1}\in \Theta_{\boldsymbol \vartheta_{M_0+1},1h}^*$. Then, by the analogous argument in the proof of Proposition (ref), the local MLEs are consistent, i.e., under $H_0: M = M_0$, for $h=1,2,...,M_0$, $\inf_{\boldsymbol{\vartheta}_{M_0+1} \in \Theta_{\boldsymbol \vartheta_{M_0+1},1h}^*} |\widehat{\boldsymbol{\vartheta}}_{M_0+1}^h - \boldsymbol{\vartheta}_{M_0+1}| \overset{p}{\to} 0$.

Consider the local LRTS for testing $H_{0,1h}: \boldsymbol\vartheta_h = \boldsymbol\vartheta_{h+1}$ defined by

where $LR^{M_0,h}_n$ converges in distribution to the random variable $(\widehat{\boldsymbol{t}}_{\boldsymbol\lambda}^h )^{\top} \boldsymbol{\mathcal{I}}^h_{\boldsymbol\lambda,\boldsymbol\nu} \widehat{\boldsymbol{t}}^h_{\boldsymbol\lambda}$, as defined in ((ref)), which is analogous to $(\widehat{\boldsymbol{t}}_{\boldsymbol\lambda} )^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol\lambda,\boldsymbol\nu} \widehat{\boldsymbol{t}}_{\boldsymbol\lambda}$ in ((ref)) for testing $H_0: M=1$.

Because $\cup_{h=1}^{M_0}\boldsymbol{\Psi}_h^*= \bar\Theta_{\boldsymbol\vartheta_{M_0 + 1}}(\boldsymbol c)$, the LRTS is identical to the maximum of the local LRTS over $h=1,2,...,M_0$:

where $\hat{\boldsymbol{\vartheta}}_{M_0+1}$ is the MLE defined by ((ref)). Then, because $LR^{M_0,h}_n\overset{d}{\rightarrow}(\widehat{\boldsymbol{t}}_{\boldsymbol\lambda}^h )^{\top} \boldsymbol{\mathcal{I}}^h_{\boldsymbol\lambda,\boldsymbol\nu} \widehat{\boldsymbol{t}}^h_{\boldsymbol\lambda}$, the asymptotic distribution of $LR^{M_0}_n$ is characterized by the maximum of $M_0$ random variables.

In practice, we implement parametric bootstrap to obtain the bootstrap p-value for testing $H_0: M= m$ against $H_1: M= m+1$.

Sequential hypothesis testing

To estimate the number of components, we sequentially test $H_0: M=m$ against $H_1: M=m+1$ starting from $m=1$, and then $m=2,\ldots,\overline M$, where $\overline M$ is the upper bound for the number of components, which is assumed to be known and larger than $M_0$. The first value for $r$ that leads to a nonrejection of $H_0$ gives our estimate for $M_0$. Robin2000 develop a similar sequential hypothesis test for estimating the rank of a matrix.

For $m=1, \ldots, \overline M$, let $c^{m}_{1-q_n}$ denote the $100(1-q_n)$ percentile of the cumulative distribution function of a random variable $ \max \{(\widehat{\boldsymbol{t}}^1_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^1_{\boldsymbol\lambda,\boldsymbol\eta} \widehat{\boldsymbol{t}}^1_{\boldsymbol\lambda} , \ldots, (\widehat{\boldsymbol{t}}^{m}_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^{m}_{\boldsymbol\lambda,\boldsymbol\eta} \widehat{\boldsymbol{t}}^{m}_{\boldsymbol\lambda} \}$ for testing $H_0: M=m$ in Propositions (ref). Let $\widehat c^{m}_{1-q_n}$ be a consistent estimator of $c^{m}_{1-q_n}$. Then, our estimator based on sequential hypothesis testing (SHT, hereafter) is defined as

The estimator $\widehat M_{\text{LRT}}$ depends on the choice of the significance level $q_n$. The following proposition states that $\widehat M_{\text{LRT}}$ converge to $M_0$ in probability as $n \rightarrow \infty$ when $-n^{-1}\log q_n=o(1)$ and $q_n=o(1)$.

Let \[ Q_n^M( \boldsymbol{\vartheta}_M):=n^{-1} \sum_{i=1}^n \log g_M(\boldsymbol{W}_i; \boldsymbol{\vartheta}_M) \text{ and } Q^M( \boldsymbol{\vartheta}_M):= \mathbb{E}[\log g_M(\boldsymbol{W}_i; \boldsymbol{\vartheta}_M)]. \]

Assumptions (ref)(a)--(e) ensure the consistency and asymptotic normality of $\widehat{\boldsymbol{\vartheta}}_M$ for $M\leq M_0$, where (c)--(e) correspond to Assumption A6 of white82em. Per Assumption (ref)(f), the Kullback--Leibler information criterion of the model relative to the true $M_0$-component model strictly decreases as the number of components $M$ increases for $M < M_0$.

Selecting the significance level $q_n$ in finite samples is challenging. Due to the complexity of analyzing the optimal choice of $q_n$, we leave this for future research and recommend presenting results across conventional levels (1%, 5%, 10%).

To compare BIC and LRT performance across different $q_n$ values, we also recommend conducting a simulation exercise similar to a parametric bootstrap; researchers can generate datasets from the estimated model with $M$ components, apply methods with different values of $q_n$, and analyze the frequency distribution of estimated components as shown in Tables (ref)--(ref) in Section (ref).

Consistent estimation of the number of components by the penalized likelihood method

We may also estimate the number of components by the penalized maximum likelihood estimator \[ \widehat{M}_{PL} = \arg\max_{M\in\{1,2,...,\overline{M}\}} p\ell_{n}^M(\hat{\boldsymbol{\vartheta}}_{M}), \] where \[ p\ell_{n}^M(\hat{\boldsymbol{\vartheta}}_{M}):=\ell_n^M(\hat{\boldsymbol{\vartheta}}_{M}) - p_{n,k_M}, \] and $\hat{\boldsymbol{\vartheta}}_{M}$ is the MLE defined by ((ref)) while $ p_{n,k_M}$ is a penalization term with $k_M:= \text{dim}(\boldsymbol{\vartheta}_{M})$ representing the number of estimated parameters for mixture models ((ref)).

The following proposition states that the penalized maximum likelihood estimator of the number of components is consistent under the regularity condition.

Assumption (ref) corresponds to Assumption (C1) of keribin00sankhya, specifying conditions for penalty functions. The BIC penalty function, given by \( p_{n,k} = \frac{k}{2}\log(n) \), satisfies this assumption, whereas the AIC penalty function, \( p_{n,k} = k\), does not satisfy Assumption (ref)(b).

Assumption (ref) is a key high-level condition preventing overestimation. In Appendix (ref), we provide explicit conditions under which Assumption (ref) holds by extending the previous asymptotic analysis, originally developed for testing $H_{0}: M = M_0$ against $H_{1}: M = M_0 + 1$, to a more general setting of testing $H_{0}: M = M_0$ against $H_{1}: M = M_1$, where $M_1 > M_0 + 1$. Under this extension, we derive conditions ensuring that $\ell_{n}^{M_1}(\hat{\boldsymbol\vartheta}_{M_1}) - \ell_n^{M_0}({\boldsymbol{\vartheta}}_{M_0}^*) = O_p(1)$ for $M>M_0+1$. The primary challenge encountered in generalizing our earlier analysis is the increased severity of the singularity in the Fisher Information matrix, which necessitates expansions of order higher than second-order in the asymptotic analysis.

To address the singularity of the Fisher information matrix, we reparameterize the model and expand the density ratio $ {g_{M}(\boldsymbol{w}; \boldsymbol\vartheta_{M})}/{g_{M_0}(\boldsymbol{w}; \boldsymbol\vartheta_{M_0}^*)} - 1$ beyond the second order. This expansion yields a quadratic-form approximation of the log-likelihood difference \( \ell_n^{M}(\hat{\boldsymbol{\vartheta}}_{M}^j) - \ell^{M_0}_n( {\boldsymbol{\vartheta}}_{M_0}^*) \) that captures higher-order terms necessary, building on a generalized version of Le Cam’s differentiable in quadratic mean (DQM) expansion liushao03as,kasahara2018arXiv. We then demonstrate that \( \ell_n^{M}(\hat{\boldsymbol{\vartheta}}_{M}^j)-\ell^{M_0}_n({\boldsymbol{\vartheta}}_{M_0}^*)=O_p(1) \), analogously to arguments used in the proof of Proposition (ref). Appendix (ref) contains the detailed discussion and Assumption (ref) presents the conditions for obtaining $\ell_n^{M}(\hat{\boldsymbol{\vartheta}}_{M}^j)-\ell^{M_0}_n({\boldsymbol{\vartheta}}_{M_0}^*)=O_p(1)$. Assumption (ref) (b) is a key condition requiring that the unique elements of ${\nabla_{\boldsymbol\theta_h} f(\boldsymbol{w}; \boldsymbol\theta_h^*)}/{g_{M_0}(\boldsymbol{w}; \boldsymbol\vartheta_{M_0}^*)}$, ${\nabla_{\boldsymbol\theta_h\otimes \boldsymbol\theta_h} f(\boldsymbol{w}; \boldsymbol\theta_h^*)}/{g_{M_0}(\boldsymbol{w}; \boldsymbol\vartheta_{M_0}^*)}$, ..., ${\nabla_{\boldsymbol\theta_h^{\otimes p_h}} f(\boldsymbol{w}; \boldsymbol\theta_h^*)}/{g_{M_0}(\boldsymbol{w}; \boldsymbol\vartheta_{M_0}^*)}$ are linearly independent and their expectation is finite for $h=1,2,...,M_0$, where the value of $p_h$ indicates a necessary order for identification.

The following proposition formalizes the results derived in Appendix (ref).

leroux92as established that, under conditions similar to Assumption (ref), the maximum penalized likelihood method produces an estimator that asymptotically does not underestimate the number of components. Building upon the work of dacunha99as, keribin00sankhya derived regularity conditions necessary for the consistency of the maximum penalized likelihood estimator. Our Assumption (ref)(b) corresponds to condition (P2) in keribin00sankhya. However, Keribin's condition (P2) is not satisfied when the number of components \( M \) is moderately larger than the true number \( M_0 \), as it relies on the sufficiency of a second-order expansion of the the density ratio for identification. In contrast, our Assumption (ref)(b) accommodates settings where parameter identification requires a higher-order rank condition, thereby extending the applicability of the consistency result.

Nonparametric estimation of the lower bound of the number of components by the rank test

We also estimate the lower bound of the number of components without imposing the parametric assumption on error distributions by extending a method proposed by Kasahara2014 which in turn is based on the rank test of Kleibergen06.

We partition the support of $Y_{t}\in\mathcal{Y}_t$ into $|\Delta_{t}|$ mutually exclusive and exhaustive subsets $\Delta_{t}=\{\delta_{1}^t, \ldots ,\delta_{|\Delta_{t}|}^t\}$ so that $\mathcal{Y}_t=\cup_{i=1}^{|\Delta_{t}|}\delta_{i}$ and $\delta_{i}\cap \delta_{j}=\emptyset$ for $i\neq j$.

For each $k\in \mathcal{T}:=\{1,2,...,T\}$, define $\mathcal{T}_{-k}:=\{s \in \mathcal{T}: s\neq t\}=\{1,..,t-1,t+1,...,T\}$ and let $\mathbf{Y}_{\mathcal{T}_{-k}}=(Y_1,..,Y_{k-1},Y_{k_1},...,Y_T)^\top$ be a vector of $Y_t$s in the group $\mathcal{T}_{-k}$. Partition the support of $\mathbf{Y}_{\mathcal{T}_{-k}}$ into $|\Delta_{\mathcal{T}_{-k}}|$ mutually exclusive and exhaustive subsets $\Delta_{\mathcal{T}_{-k}}=\{\boldsymbol\delta_{1}^{\mathcal{T}_{-k}}, \ldots ,\boldsymbol\delta_{|\Delta_{\mathcal{T}_{-k}}|}^{\mathcal{T}_{-k}}\}$. In particular, we construct this partition by unfolding the tensor product $\otimes_{t\neq k} \Delta_{t}$ using the the Khatri-Rao product denoted by $\odot$ as $\Delta_{\mathcal{T}_{-k}}=\odot_{t\neq k} \Delta_{t}= \Delta_{1}\odot \cdots \odot \Delta_{k-1}\odot \Delta_{k+1}\odot \cdots\odot \Delta_T$ so that $|\Delta_{\mathcal{T}_{-k}}|=|\Delta_t|^{T-1}$.

For each $k\in \mathcal{T}$, we construct a $|\Delta_{t}| \times |\Delta_{\mathcal{T}_{-k}}|$ bivariate probability matrix $\boldsymbol P_{k}$ by arranging $\Pr(Y_k \in \delta_{a},\boldsymbol Y_{-k} \in \boldsymbol\delta_{b}^{\mathcal{T}_{-k}})$ for partition level $(a,b)=(1,1),\ldots,(|\Delta_{t}|,| \Delta_{{\mathcal{T}_{-k}}}|)$ as

Collect the marginal probability distribution of ${Y}_{k}$ and $\mathbf{Y}_{-k}$ over $\Delta_{t}$ and $\Delta_{\mathcal{T}_{-k}}$ conditional on being from the $j$-th component into a vector as

Then, under the conditional independence assumption as in the mixture model ((ref)) but without imposing parametric restrictions, we may represent $\boldsymbol P_k$ as \[ \boldsymbol P_k = \sum_{j=1}^M \alpha^j \boldsymbol p_k^j (\boldsymbol q_k^j)^\top. \] Kasahara2014 shows the rank of $\boldsymbol P_k$ identifies the lower bound of the number of components and develop a sequential hypothesis testing procedure for estimating the rank of $\boldsymbol P_k$ when the empirical quantile of the $Y_t$'s are used to construct the partition. Because there are $T$ possible ways to pick different $k$'s out of $\{1,...,T\}$, we test the maximum of the ranks of $\boldsymbol P_k$ across $k=1,...,T$.

We first develop a rk-statistic of Kleibergen06 for testing the null hypothesis of $\text{rank}(\boldsymbol P_k)=r$. Write the eigenvalue decomposition of ${\boldsymbol P}_k$ as \[ {\boldsymbol P}_k = \mathbf{U}^k\mathbf{S}^k(\mathbf{U}^k)^\top =

^\top, \] where $\mathbf{U}^k$ is a $|\Delta_t| \times |\Delta_t|$ orthonormal matrix and $\mathbf{S}^k$ is a diagonal matrix containing the eigenvalues of $\mathbf{P}_k$ in decreasing order. In the partition of $\mathbf{U}^k$ and $\mathbf{S}^k$ on the right-hand side, $\mathbf{U}^k_{11}$ and $\mathbf{S}^k_1$ are $r \times r$, and the dimensions of the other submatrices are defined conformably. Then, the null hypothesis $\mathcal{H}_0: \mathrm{rank}(\mathbf{Q}_k) = r$ is equivalent to $\mathcal{H}_0: \mathbf{S}^k_2 = 0$. The statistic of Kleibergen06 is based on an orthogonal transformation of $\mathbf{S}^k_2$ given by $\mathbf{\Lambda}^k_r = (\mathbf{A}^k_{r })^\top\mathbf{Q}_k\mathbf{A}^k_{r }$, where \[ \mathbf{A}^k_{r } =

(\mathbf{U}^k_{22} )^{-1}(\mathbf{U}^k_{22}(\mathbf{U}^k_{22})^\top)^{1/2}. \]

Let $\widehat{\boldsymbol P}_k$ be a sample analogue estimator for $\boldsymbol P_k$ for which we have $\sqrt{n}\text{vec}(\widehat{\boldsymbol P}_k-\boldsymbol P_k)\overset{d}{\rightarrow} N(0,\boldsymbol\Sigma_k)$. The following proposition follows from Theorem 1 of Kleibergen06.

Kleibergen06 proposed the statistic called the rk-statistic: \[ \text{rk}^k(r) = n (\widehat{\boldsymbol{\lambda}}_r^k)^\top (\widehat{\boldsymbol{\Omega}}^k_r) ^{-1}\widehat{\boldsymbol{\lambda}}_r^k, \] where $\widehat{\boldsymbol{\Omega}}^k_r$ is a consistent estimator for $\boldsymbol{\Omega}^k_r$. If the assumptions of proposition (ref) hold, $\mathrm{rk}(r)$ converges in distribution to a $\chi^2(|\Delta_t|-r)$ random variable as $n \rightarrow \infty$.

When $T \geq 3$, we test the null hypothesis that $\text{rank}(\boldsymbol P_k) \leq r$ for each $k = 1,\dots,T$. By selecting partitions such that $|\Delta_t| = r + 1$, we define the following ave- and max-rk test statistics:

See Section 3.4 of Kasahara2014 for the asymptotic distribution of these statistics.

In practice, we use the Bayesian bootstrap to obtain the bootstrap p-value for testing the rank, rejecting the null hypothesis $\text{rank}(\boldsymbol P_k)\leq r$ for all $k=1,\dots,T$ at significance level $\alpha$ if the bootstrap p-value is strictly less than $\alpha$; see Appendix (ref). We estimate the lower bound for the number of components by applying the sequential hypothesis testing procedure described in Section 3.2 of Kasahara2014.

Simulation

This section evaluates the relative performance of several approaches—rank tests based on ave-rk and max-rk statistics, LRT , AIC, and BIC—in estimating the number of components across different models through simulation studies. For our simulation and empirical analyses, we set the constraint that $\alpha_j, \tau_{jk} \geq 0.05$ and $\sigma_j \geq 0.05 \hat{\sigma}_0$, effectively setting $c_1=0.05$ and $c_2=0.05$. We estimate the parameter by the MLE under these constraints using the EM algorithm as described in Appendix (ref).

Tables (ref) and (ref) summarize size and power performance of statistical methods for testing $H_0: M=2$ against $H_1: M=3$ in models with conditionally independent errors without covariates, as defined by equations ((ref)) and ((ref)), based on $500$ simulations with $T=3$. Results compare ave-rk, max-rk, LR, AIC, and BIC methods, where we report the selection frequencies between two- and three-components models for AIC and BIC.

Table (ref) presents rejection frequencies at the 5% significance level, using data from two-component models without covariates. The ave-rk and max-rk tests achieve near-nominal levels for larger samples ($N=400$) and distinct means ($\mu=(-1,1)$) but underperform with smaller samples or closer means ($\mu=(-0.5,0.5)$). The LRT consistently maintains near-nominal performance. AIC often overestimates component counts, while BIC remains conservative, rarely overestimating.

Table (ref) evaluates power using data from three-component models. The LRT exhibits excellent power (100%) with clearly separated means ($\mu=(-1.5,0,1.5)$) for both sample sizes ($N=200,400$), outperforming ave-rk and max-rk tests even when mean separations narrow ($\mu=(-0.5,0,1.5)$). While the BIC closely matches LRT performance under clear mean separation, it becomes overly conservative as the means become closer, selecting the true three-component model less frequently than the LRT. The AIC chooses the three-component models more often than other methods, reflecting its tendency to overestimate the number of components.

We also evaluate sequential testing procedures under realistic conditions by simulating datasets from three-component models estimated from real data. Specifically, we consider different sample sizes and various specifications for density functions and stochastic processes governing the error terms. This approach enables us to assess how model misspecification impacts the accuracy of estimating the number of components.

Table (ref) summarizes the selection frequencies (in percentages) for identifying the number of components (from 1 to 6) using the AIC, BIC, LRT, and ave-rk and max-rk tests in simulations for the sample sizes of $n=50$ and $225$ with $T=3$. Data are generated from a three-component model without covariates under normal error density, as defined by equations ((ref)) and ((ref)), estimated using Chilean fabricated metal products industry data. The model features mixing probabilities $\boldsymbol\alpha = [0.352, 0.402, 0.245]$, means $\boldsymbol\mu = [-1.01, -0.557, -0.242]$, and standard deviations $\boldsymbol\sigma = [0.464, 0.187, 0.195]$. Results are reported for correctly specified normal errors ("Normal") and over-specified two-component normal mixtures ("Mixture"), with LRTs evaluated at significance levels $q_n=0.01$, $0.05$, and $0.10$ while the ave- and max-rk tests implemented at 5% significance level. The parameters are estimated by the MLE under the constraint $\alpha_j, \tau_{jk} \geq 0.05$ and $\sigma_j \geq 0.05 \hat{\sigma}_0$ but the simulation results are not so sensitive to changing these tuning parameters as shown in Tables (ref)--(ref) in Appendix (ref).

Panel A of Table (ref) reports the results for a small sample size of $n=50$, showing that both the BIC and LR tests effectively identify the correct number of components when the model is correctly specified with normal error density. Specifically, the LR tests select $M = 3$ most frequently, with accuracy improving as the significance level increases: 74% at $q_n = 0.01$, 82% at $q_n = 0.05$, and 83% at $q_n = 0.10$, whereas BIC selects the correct model ($M = 3$) in 73% of cases. Thus, BIC tends to underestimate the number of components relative to the LR tests under correct specification. On the other hand, AIC selects the correct number of components in 84% of cases, although it overestimates the number in 12% of cases.

Thus, for the correctly specified models, the LR test outperforms BIC with similar component means and small samples, as BIC underestimates component numbers. While BIC is simpler to implement, bootstrapping makes the LR test equally easy. Using all three methods provides a useful range: BIC as the lower bound, AIC as the upper bound, and LR test as a middle ground. Since rank tests typically underestimate component numbers, findings of $M\geq 2$ provide strong evidence for heterogeneity.

Under the over-specified scenario with two-component normal mixture error density specification, both BIC and the LR test select $M=2$ over $M=3$ more frequently compared to the correctly specified normal-density scenario. The LR test is particularly conservative at lower significance levels, choosing $M=2$ in 74% of cases at $q_n=0.01$, but its accuracy for selecting $M=3$ improves to 64% at $q_n=0.10$. The BIC selects $M=2$ in 58% of cases and $M=3$ in 40%. AIC consistently favours more complex models, selecting $M=4$ or higher in 33% of cases.

The ave-rk and max-rk tests are highly conservative, selecting $M=2$ in nearly all cases (90% and 89%, respectively), and rarely identifying the true number of components.

Overall, the LR test and BIC perform comparably, though BIC sometimes underestimates component numbers when mixture components have similar means and sample sizes are small. While BIC is simpler to implement, bootstrapping makes the LR test equally easy. Using all three methods provides a useful range, especially with small samples: BIC as the lower bound, AIC as the upper bound, and the LR test as a middle ground. Furthermore, since rank tests typically underestimate component numbers, findings of $M\geq 2$ also provide strong evidence for heterogeneity.

Panel B of Table (ref) presents results for the sample size matches the actual data at $n=225$ instead of $n=50$, where we observe a substantial improvement in the performance of BIC and LRT as the sample size increases from $n=50$ to $n=225$. Both methods correctly select $M=3$ in approximately 90% or more of cases, even under the over-specified mixture density scenario. In contrast, the AIC continues to favor more complex models than BIC or LRT, while the ave-rk and max-rk tests tend to select $M=2$ over $M=3$ in the majority of cases.

Table (ref) examines the case in which the true data-generating process is a three-component model with a normal mixture error density, as described in ((ref)) with ((ref)), rather than a simple normal error density. Data are simulated using parameter values estimated from the Chilean fabricated metal products industry, as detailed in the table notes, for sample sizes of $n=50$ and $225$.

Similar to the results obtained with the normal error density DGP, both BIC and the LR test tend to frequently select $M=2$ over $M=3$ under the normal mixture error density specification, compared to a normal error density scenario at the small sample size of $n=50$. As the sample size increases from $n=50$ to $n=225$, the accuracy of BIC and LR tests significantly improves, particularly for BIC and the LR test at $q_n=0.01$. At $n=225$, however, the LR test under the incorrect normal error density specification tends to select $M=4$ over $M=3$ more frequently than the correctly-specified normal mixture scenario, suggesting that incorrectly imposing a normal error density may lead to overestimating the number of components with the sufficiently large sample size. The AIC has a higher tendency to select $M=4$ over $M=3$ than the BIC or the LRT while the rank-based tests continue to under-estimate the number of components.

In the Appendix, Tables (ref)–(ref) present simulation results based on estimated parameters for the food products and textiles industries, revealing patterns similar to those in Tables (ref)–(ref).

Appendix (ref) also presents the simulation results for dynamic panel mixture models when the error terms follow AR(1) processes whose innovations are drawn from either a normal distribution or a two-component normal mixture, as defined by equations ((ref)) with ((ref)) or ((ref)), respectively, where both BIC and LRT performs well under the realistic DGP based on an estimated model from the Chilean fabricated metal products industry. See Tables (ref)-(ref).\footnote{Furthermore, simulation results using data generated from the estimated models with labor-augmented technological changes are presented in Tables (ref) and (ref).}

Empirical analysis: plant heterogeneity in production function

In this section, we investigate heterogeneity beyond Hicks-neutral technological differences using finite mixture specifications ((ref)), ((ref)), and ((ref)), as illustrated in Examples (ref), (ref), and (ref). Specifically, we estimate the number of latent technology types characterized by variations in input elasticities and labor-augmented technological changes through the application of the Likelihood Ratio Test (LRT), Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and the nonparametric rank test based on ave-rk statistics.

For our analysis, we use plant-level panel data from Chilean manufacturing plants covering the period 1992--1996.\footnote{Please refer to KASAHARA2008 for details on the dataset.} We focus on the three largest industries---Food Products, Fabricated Metal Products, and Textiles---and include observations of plants with continuous data entry for either $T=3$ or $T=5$ consecutive years within each industry. We use data from 1994--1996 for estimating the model with conditionally independent errors, while data from 1992--1996 are used for dynamic panel mixture models. The selected panel lengths correspond to the minimum length required for non-parametric identification of the respective panel mixture models.\footnote{Non-parametric identification requires a panel length of $T=3$ for models with conditionally independent errors Kasahara2009 and $T=5$ for dynamic panel mixture models Hu2012,Kasahara2022esri.} Table (ref) presents summary statistics for the revenue share of materials and the log of gross output among other variables in these industries. The note section of Table (ref) further explains the sample selection criteria.

Plant heterogeneity in output elasticity of material input

Empirical studies often assume a Cobb-Douglas production technology $\log O_{it} = \gamma_0 + \gamma_v V_{it} + \gamma_\ell \log L_{it} + \gamma_k \log K_{it} +\gamma_z \log Z_{it} + \omega_{it}$, where $O_{it}$, $V_{it}$, $L_{it}$, $K_{it}$, and $Z_{it}$ represent output, material input, labor, capital, and other observed characteristics, respectively. Because input elasticities are constant under Cobb-Douglass specification, the first-order condition for profit maximization ((ref)) in Example (ref) implies

i.e., the ratios of material input to output is constant, where $P_{O,t}$ and $P_{V,t}$ are output and material input prices, respectively.

The implication of Equation ((ref)) that material input shares are identical across plants can be empirically tested. Figure (ref) displays histograms illustrating the distribution of plant-level material input shares in the Food Products, Fabricated Metal Products, and Textiles industries, clearly revealing substantial variability across plants within each industry. In Column 4 of Table (ref), standard deviations of the revenue shares attributable to material costs are considerable, ranging from $0.145$ to $0.216$, even within narrowly defined industries. Consequently, the hypothesis that the ratio of material input to output is constant across plants is overwhelmingly rejected by the data.

While the substantial variability in material input shares across plants suggests potential plant heterogeneity in the output elasticity of material input, the implication derived from equation ((ref)) is strictly valid only under the Cobb-Douglas assumption. Under more general production function specifications, such as CES or translog forms, output elasticities with respect to inputs naturally depend on the levels of materials, labor, capital, and other observable plant characteristics, even in the absence of technological heterogeneity. Additionally, the relationship between material input-to-output ratios and the coefficient $\gamma_v$ may be influenced by idiosyncratic shocks.\footnote{As Gandhi2020 demonstrate, if the production function is specified as $\log O_{it} = \gamma_0 + \gamma_v V_{it} + \gamma_\ell \log L_{it} + \gamma_k \log K_{it} +\gamma_z \log Z_{it} + \omega_{it}+\epsilon_{it}$, where $\epsilon_{it}$ is an i.i.d. mean zero random variable whose realization is unknown when the intermediate input $V$ is selected, then the first-order condition for profit maximization yields $\frac{P_{V,t} V_{it}}{P_{O,t} O_{it}}=\gamma_v \mathbb{E}[e^{\epsilon_{it}}]e^{-\epsilon_{it}}$. Consequently, the relationship between material input-to-output ratios and the coefficient $\gamma_v$ is subject to idiosyncratic shocks.}

To investigate persistent plant heterogeneity in output elasticities of material input within a more general production function framework while accounting for transitory idiosyncratic shocks, we specify that the logarithm of the output elasticity of material input is linearly related to observed plant characteristics ($\mathbf{x}_{it}$), such as the log of capital and other plant-specific attributes, conditional on the latent technological type of the plant, denoted by $D_i = j$, as

where $\boldsymbol{\beta}_j$ represents parameters specific to the $j$-th latent technology type, and $\epsilon_{it}$ denotes an idiosyncratic shock to plant-level input elasticities. The distribution of $\epsilon_{it}$ is type-specific and can be represented as either a single-component normal distribution or a mixture of normals: \[ \epsilon_{it} \overset{iid}{\sim} \mathcal{N}(\mu_{j},\sigma_j^2) \quad \text{or} \quad \epsilon_{it} \overset{iid}{\sim} \sum_{k=1}^{{\cal{K}}} \tau_{jk}\, \mathcal{N}(\mu_{jk},\sigma_j^2)\quad\text{given $D_i = j$}, \] where ${\cal{K}}$ represents the number of components in normal mixture error distributions.

This formulation leads directly to the finite mixture model given in equation ((ref)) with component density functions as outlined in equations ((ref)) or ((ref)). By estimating the number of latent technology types ($M$), we can quantify the degree of heterogeneity in output elasticities of material input among plants, even after controlling for observable plant characteristics.

To determine the number of latent technology types, we conduct the sequential hypothesis test using the LRT, as outlined in Section (ref), setting $q_n=0.01$ guided by the simulation results presented in Section (ref) although the results are not sensitive to an alternative choice of $q_n=0.05$. Additionally, we apply the AIC, the BIC, and the nonparametric rank test in determining the number of technology types.

We conduct analyses both without covariates and with covariates including the log of capital ($\log K$), imported material shares (import), and 4-digit industry classification dummies (CIIU4). To allow for flexible error structures and relax the normality assumption, we model the error terms using normal distributions as well as two-component and three-component normal mixtures.

Panel A of Table (ref) provides a comparative analysis of the estimated number of technology types across the three industries, utilizing the AIC, BIC, LRT, and nonparametric rank test for models with conditionally independent errors as specified in equation ((ref)). As the flexibility of the error term structure increases from normal to multi-component mixtures, the number of components estimated generally decreases for the AIC, BIC, and LRT. For example, without covariates and assuming normally distributed error terms, the LRT estimates a high number of components—8 for Metal, 10+ for Food, and 9 for Textiles. However, when employing a more flexible three-component normal mixture error density specification, the number of components reduces to 4, 8, and 6 for Metal, Food, and Textiles, respectively. This is consistent with our simulation results reported in Table (ref), indicating that a less flexible parametric assumption on the error term distribution may lead to over-estimation of the number of components. Including covariates such as $\log K$, imported material shares, and controls for 4-digit industry classifications (CIIU4) further decreases the estimated number of components for the LRT to 4, 5, and 5 for Metal, Food, and Textiles, respectively.

Comparing the results of the LRT with those from the AIC and the BIC, we find that the BIC estimates a number of types similar to the LRT, whereas the AIC typically estimates a higher number. The number of technology types suggested by both the LRT and the BIC remains high for models with covariates and normal mixture error densities. Furthermore, as reported in the first row of Table (ref), the lower bound of the number of components estimated by the nonparametric rank test is 3 for all industries. Overall, these results provide strong evidence of substantial and persistent plant heterogeneity in material input elasticities across years.

Figure (ref) presents histograms of estimated error terms alongside component-specific error distributions (solid lines) for models incorporating covariates that control for differences in $\log K$, import shares, and 4-digit industry dummies. The models employ iid error distributions specified as three-component normal mixtures ($\mathcal{K}=3$). Plant observations are classified into distinct technology types based on posterior probabilities calculated using Bayes' theorem, with each type indicated by a unique color. Following the LRT model selection criterion, we present a 4-component model for Fabricated Metal, a 5-component model for Food Products, and a 5-component model for Textiles. The results demonstrate clear evidence of non-normality in component-specific error distributions, while the flexible normal mixture specification effectively captures this non-normality within each component, illustrating the advantages of our modeling approach.

To better understand the nature of heterogeneity, Figures (ref) and (ref) present the estimated parameters and the 95% confidence intervals (CI) for two-component models with covariates, where the error densities are modeled using normal and two-component normal mixture distributions, respectively. The CI is computed based on the robust standard errors computed via the sandwich variance estimator to account for potential misspecification. We specifically focus on the two-component models because confidence intervals for some parameters become excessively wide in models with three or more components.\footnote{Parameter estimates for three-component models are provided in Figures (ref) and (ref) in Appendix (ref).}

In Figure (ref), under the normality assumption, the two latent types differ significantly in their estimated mixing proportions, means, and variances, whereas no statistically significant differences are found in the estimated coefficients for $\log K_{it}$ and import shares. Specifically, compared to the first type, the second type exhibits higher mixing proportions, higher means, and lower variances consistently across all three industries. The absence of statistical evidence for differences in $\beta_{\log K}$ and $\beta_{\text{import}}$ suggests that the Cobb-Douglas specification with two-component mixtures adequately captures the observed persistence in material shares.

Figure (ref) presents parameter estimates for two-component models under the normal mixture density, reporting averaged mean parameters and average standard deviation parameters as $\hat\mu_j = \sum_{k=1}^2 \hat\tau_{jk}\hat\mu_{jk}$ and $\widehat{Var}(\epsilon_{it}| D_i=j)=\hat\sigma_j^2+\hat\tau_{j1}(1-\hat\tau_{j1}) (\hat\mu_{j1}-\hat\mu_{j2})^2$, respectively, where bootstrap is used to construct the confidence intervals for these two parameters. Mixing proportions and mean differences between latent types resemble those under normal density, while $\beta_{\log K}$ and $\beta_{\text{import}}$ are insignificant for both types. Wider confidence intervals reflect model complexity.

As discussed in Example (ref), we also consider the case where the error term $\epsilon_{it}$ in equation ((ref)) follows an AR(1) process with $\epsilon_{it}=\rho_j \epsilon_{it-1}+\xi_{it}$ to capture potential persistence in shocks to output elasticities of material input. This leads to the following specification:

where $\xi_{it}$ is drawn from either a normal distribution or multi-component normal mixtures as specified in equation ((ref)). By further specifying the initial distribution of $\log\left(\frac{P_{V,t} V_{it}}{P_{O,t} O_{it}}\right)$ at $t=1$ as:

where $\epsilon_{i1} \overset{iid}{\sim} \mathcal{N}(\mu_{1,j},\sigma_{1,j}^2)$ or $\sum_{k=1}^{{\cal{K}}} \tau_{1,jk}\, \mathcal{N}(\mu_{1,jk},\sigma_{1,j}^2)$ given latent technology type $D_i = j$, this specification yields the conditional density of $\{y_t\}_{t=1}^T$ and $\{\boldsymbol{x}_{t}\}_{t=1}^T$ in the form of equation ((ref)) with ((ref)), where the component-specific densities are given by either equation ((ref)) or ((ref)). In this framework, each component-specific technology type represents a distinct stochastic process governing the output elasticities of material input.

Panel B of Table (ref) reports the estimated number of components selected by the AIC, BIC, and LRT for models specified in equations ((ref))–((ref)).\footnote{Parameters are estimated using the EM algorithm detailed in Appendix (ref), explicitly incorporating the constraint implied by specification ((ref)).} Due to the incorporation of persistence in the error term $\epsilon_{it}$ within the AR(1) framework, the results indicate fewer estimated components compared to models assuming conditionally independent errors. Nonetheless, there remains clear evidence of heterogeneity in the stochastic processes governing the error term $\epsilon_{it}$ especially for Food and Textiles industries; under the specification with covariates and two-component normal mixture innovation densities, the LRT identifies 2, 2, and 2 components for the Metal, Food, and Textiles industries, respectively, whereas the BIC selects 2, 5, and 2 components.

Figure (ref) reports parameter estimates for the AR(1) specification with covariates and two-component normal mixture innovation densities.\footnote{For estimating models with an AR(1) specification, we imposed a restriction that the parameter $\mu_j$ lies within the interval $[-2, 2]$, as the estimated value of $\mu_1$ for the food industry was implausibly low when the model was estimated without this bound.} For Fabricated Metal Products and Textile industries, relative to the first type, the second type exhibits higher mixing proportions, higher means, lower variances, and higher persistence (as captured by the AR(1) coefficient $\rho_j$). Reflecting model complexity, the confidence intervals are generally wide, and differences between latent types are often not statistically significant.

Plant heterogeneity in the stochastic process of factor-augmented technological changes

Factor-augmented technological changes are widely used to capture heterogeneity in production functions beyond the standard Hicks-neutral approach doraszelski2018measuring,Zhang2019,Raval2019. The heterogeneity in the stochastic process of factor-augmented technological changes can be effectively modeled with finite mixture dynamic panel data models.

As discussed in Example (ref), we extend the production function ((ref)) by incorporating labor-augmented technological change \(\delta_{it}\) Demirer2022: \[ O_{it}=e^{\omega_{it}}\bar F_{jt}(\chi_{jt}(V_{it}, e^{\delta_{it}} L_{it}), K_{it}, Z_{it},\epsilon_{it})\ \text{ with }\ \chi_{jt}(V_{it}, e^{\delta_{it}} L_{it}) := \left[\alpha_{V,jt} V_{it}^{1/\varsigma_j} + \alpha_{L,jt} (e^{\delta_{it}} L_{it})^{1/\varsigma_j}\right]^{\varsigma_j}, \] and \(\delta_{it}\) evolves according to an AR(1) process: \[ \delta_{it} = \rho^{\delta}_j \delta_{it-1} + \eta^{\delta}_{it}, \] where \(\eta^{\delta}_{it}\) follows either a normal or a multi-component normal mixture distribution.

Assuming both material and labor inputs are flexibly chosen after observing relevant prices and shocks, profit maximization implies a dynamic relationship for log input ratios:

where $\tilde\eta_{it}^\delta = \eta^{\delta}_{it}/(1-\varsigma_j)$ and $\mu^{\delta}_{jt}$ depends on $\alpha_{V,jt}$, $\alpha_{L,jt}$, $\varsigma_j$, and the input prices. By specifying the density functions of the initial values \( \log(V_{it}/L_{it}) \) at $t=1$, this specification leads to conditional densities as described in equations ((ref)) and ((ref)), with component-specific densities detailed in either ((ref)) or ((ref)).

Table (ref) presents the results of model selection for the dependent variable $\log(V/L)$ (the share of material inputs over labor) across three industries—Metal, Food, and Textiles—using AIC, BIC, and Likelihood Ratio (LR) metrics. The table compares various model specifications, including those without covariates, those with conditioning variables ($\log K$, import shares, and 4-digit industry classification dummies). Both static models and models with autocorrelated error terms (AR1) are analyzed, and the error term structures range from normal distributions to more flexible two-component, three-component, and 4-component normal mixtures.

For the model without mixture components or covariates, AIC overestimates the number of components, predicting $10$ components for all industries. In contrast, BIC is more conservative, predicting $3$ components for Metal and Textiles, and $6$ for Food. More flexible mixture error structures generally lead to fewer predicted components, particularly under BIC and LR metrics. Introducing mixture error structures reduces the number of predicted components across all metrics. For example, under the Plain Mixture 2-Component model, AIC predicts $6$ components for Metal, $2$ for Food, and $2$ for Textiles, while BIC predicts $3$, $2$, and $2$ components, respectively. LR aligns more closely with BIC, predicting $2$, $4$, and $2$ components, respectively. Including covariates ($\log K$) and CIIU4 controls in the models reduces the predicted number of components in most cases, particularly when using BIC and LR. For example, under the $lnK$, Import, 2-Component Control CIIU4 model, LR predicts $1$ components for Metal, $3$ for Food, and $1$ for Textiles, compared to $2$, $4$, and $2$ under the corresponding Plain Mixture model.

Conclusion

This paper develops statistical methods for determining the number of components in panel data finite mixture regression models with normally distributed or more flexible normal mixture errors. We establish that panel data structures eliminate higher-order degeneracy problems present in cross-sectional normal mixture models while retaining issues of unbounded likelihood and infinite Fisher information. We address these challenges and derive the asymptotic null distribution of the LRT statistic, demonstrate the consistency of BIC and inconsistency of AIC, and propose a sequential hypothesis testing approach for consistent component selection.

The empirical application to Chilean manufacturing data provides compelling evidence of substantial plant heterogeneity in production technology. Using finite mixture models with conditionally independent normal mixture errors, we identify 4 or 5 distinct technology types across three major industries: Fabricated Metal Products, Food Products, and Textiles. Our analysis reveals significant variation in output elasticities of material inputs across plants within narrowly defined industries, with nonparametric rank tests consistently identifying at least three components across all industries. Dynamic panel analysis incorporating lagged dependent variables as covariates suggests that while some heterogeneity reflects transitory variations, substantial permanent technological differences persist across plants.

These findings contrast sharply with standard production function estimation methods that impose homogeneous coefficients across plants. Our results highlight the importance of explicitly accounting for unobserved plant heterogeneity beyond Hicks-neutral technological differences, with significant implications for policy analysis where homogeneous assumptions may mischaracterize productivity distributions.

The methodological framework extends beyond production function analysis to other contexts with latent group structures in panel data. Future research could extend the analysis to more general production function specifications, develop methods for simultaneous component and distributional selection, and investigate relationships between estimated technological types and observable plant characteristics. Our work contributes both methodologically and empirically to understanding plant heterogeneity, providing useful tools for uncovering latent structures while demonstrating substantial technological diversity that challenges common assumptions in empirical production analysis.