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.
130,721 characters · 14 sections · 55 citation commands
Testing the Number of Components in Finite Mixture Normal Regression Models with Panel Data 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]
Finite mixture models offer a natural representation of heterogeneity across a finite number of classes. Because of their flexibility, they have been used in empirical applications in various fields since the proposal of a two-component normal mixture model by Pearson1894. In economics, finite mixtures are frequently used to model unobserved individual-specific effects in labor economics, health economics, and industrial organization, as well as in other fields.\footnote{For example, Heckman1984 use the finite mixture model to provide an alternative method of accounting for unobserved heterogeneity in the analysis of the single-spell duration times of unemployed workers. Keane1997 and Cameron1998 analyze a dynamic model of schooling and occupational choices with unobserved heterogeneous human capital. Likewise, finite mixture models have been applied in health economics. Deb1997 develop a finite mixture negative binomial count model that accounts for the unobserved dispersion of medical care utilization by the elderly. Kamakura1989 and Andrews2003 model consumer segmentation in marketing in industrial organizations.} The theoretical properties of finite mixture models and examples of their applications have been discussed by numerous authors, such as Titterington1985, Lindsay1995, and McLachlan2000.
The number of components is a crucial parameter in finite mixture models. In economic applications, the number of components often represents the number of unobservable types or abilities of individuals. Choosing an arbitrary number of parameters may lead to overestimation or underestimation of the level of heterogeneity. Using too few components may result in biased estimation, while using too many components can be computationally costly and the model becomes singular because of identification problems. Thus, developing a statistical procedure to determine the number of components is essential.
Testing for the number of components in normal mixture regression models has long been an unsolved problem. The regularity conditions of the likelihood ratio test (LRT) for standard asymptotic analysis fail in finite mixture models because of issues such as non-identifiable parameters, the singularity of the Fisher information matrix, and the true parameter being on the boundary of the parameter space. Numerous papers have been written on the subject of the LRT for the number of components ghoshsen85book, chernofflander95jspi, lemdanipons97spl, chenchen01cjstat, chenchen03sinica, cck04jrssb, garel01jspi, garel05jspi, Chen14joe, and the asymptotic distribution of the LRT statistic for general finite mixture models has been derived as a function of the Gaussian process dacunha99as, liushao03as, zhuzhang04jrssb, azais09esaim. However, the key assumptions in these works are violated in cross-sectional normal regression models because normal mixtures possess additional undesirable mathematical properties: (i) the Fisher information for testing is not finite, (ii) the log-likelihood function is unbounded, and (iii) the second derivative of the density function for the mean parameter is linearly dependent on its first derivative for the variance parameter. Kasahara2015a analyze the asymptotic distribution of the LRT statistics of a cross-sectional univariate finite mixture normal regression model, and Kasahara2019 develop a multivariate extension. Amengual2022 develop a score-type test for a cross-sectional normal mixture model.
This paper develops a likelihood ratio-based test for determining the number of components in finite mixture normal regression models with panel data, where outcome variables are conditionally independent across periods given the latent type within each unit. To the best of our knowledge, it is not known in the literature whether the aforementioned problems (i)--(iii) of the cross-sectional normal mixture still arise in the panel normal mixture. 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 the panel normal regression mixture models with conditional independent errors.\footnote{Kasahara2014 develop a procedure to estimate a lower bound on the number of components consistently in finite mixture models in which each component distribution has independent marginals, which includes panel normal regression mixture models with conditionally independent errors as a special case.}
We show that problems related to (i) and (ii) arise, but the higher-order degeneracy of problem (iii) disappears in panel normal mixture models with conditional independence. Following chenli09as and Kasahara2015a, we consider a penalized likelihood ratio test (PLRT) and an expectation-maximization (EM) test to deal with the unboundedness and analyze the asymptotic distribution of the PLRT 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 penalized likelihood ratio test statistic (PLRTS) and the EM test statistic are characterized by the maximum of $M_0$ random variables, which we can easily simulate. Building on the PLRT and EM tests, we propose a sequential hypothesis testing approach for consistently estimating the number of components. In simulations, our proposed PLRT and EM tests demonstrate favorable finite sample properties. Moreover, a sequential hypothesis testing approach accurately selects the correct number of components with high frequency, surpassing selection procedures based on the Akaike information criterion (AIC) and the Bayesian information criterion (BIC).
This paper makes several contributions to the literature. First, it analyzes the likelihood ratio-based test for the number of components in panel normal regression mixture models with conditional independence. Kasahara2015a and Kasahara2019 analyze likelihood ratio-based tests for the number of components in cross-sectional univariate normal mixture regression models and multivariate normal mixture models, respectively. We demonstrate that the asymptotic distribution of the PLRT and EM tests for panel normal regression mixture models with conditionally independent errors differs from that of the univariate/multivariate normal mixture models in the aforementioned two papers because higher-order dependency does not occur when the repeated measurement of the outcome variables is available in the panel data. Furthermore, we develop a sequential hypothesis testing approach for consistently estimating the number of components.
Second, while it is well known that the log-likelihood function of normal mixture models is unbounded Hartigan1985, it is unknown whether the related unboundedness problem arises in panel data. We show that the likelihood ratio test statistic is unbounded in panel normal mixture models with conditionally independent errors when the time dimension of panel data is finite. This unboundedness causes over-rejection of the LRT. We introduce a penalty function to prevent the likelihood ratio test statistics from being unbounded, using computational experiments to determine the data-driven penalty function. We also develop an R package {NormalRegPanelMixture} Hao2017 that contains the EM test module and asymptotic distribution simulation module.
Third, we conduct an empirical analysis of the number of production technology types using panel data from Japanese and Chilean manufacturing firms and provide strong evidence of substantial heterogeneity in production function coefficients across firms within narrowly defined industries. This is an important contribution to the literature on production function estimation, as most empirical applications assume the homogeneity of production function coefficients across firms using the standard production function estimation methods developed by olley1996dynamics, levinsohn2003estimating, and Ackerberg2015. Our empirical finding suggests that it is essential to incorporate unobserved heterogeneity into the production function coefficients across firms in applications LiSasaki17arxiv, doraszelski2018measuring, Balat19mimeo, Kasahara2022esri.
The EM test approach was introduced by lcm09bm and chenli09as to test homogeneity in finite mixture models. Li2010 develop an EM test for the null hypothesis of $M_0$ components applicable to general $M_0\geq 2$, and Kasahara2015a propose an EM test for normal regression mixture models to test the null of $M_0\geq 2$. The EM approach is also applied to test homogeneity in multivariate mixtures Niu2011 and subgroup analyses Shen2015. More recently, Liu2018 extend the EM test to mixtures of the general location-scale family distribution, and Kasahara2019 develop an EM test for multivariate normal mixture models. Building upon the literature, this paper develops an EM test for panel normal regression mixture models with conditionally independent errors.
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.
The rest of this paper is organized as follows. In Section (ref), we define the finite normal mixture panel regression model. In Section (ref), we demonstrate the PLRT for testing the homogeneity of a normal mixture panel regression against a two-component model as a precursor to obtaining the general test of $M_0$ components. Section (ref) generalizes the result to testing $M_0$ components against $M_0 + 1$ components. Section (ref) introduces the EM test for testing $M_0$ components against $M_0 + 1$ components. Section (ref) derives the asymptotic distribution of the PLRT and EM tests under local alternatives. Section (ref) develops a consistent estimator for the number of components based on sequential hypothesis testing. Section (ref) presents the simulated results of the tests. Section (ref) provides an empirical application. In the following, $:=$ denotes “equals by definition,” and boldface letters denote vectors or matrices.
We consider finite mixture normal regression models with panel data, where the panel length $T$ is fixed and the number of cross-sectional observations $n$ goes to infinity. Define $\boldsymbol{w} := \{y_{t},\boldsymbol{x}_{t},\boldsymbol{z}_{t}\}_{t=1}^T$ with $y_t \in \mathbb{R}, \boldsymbol{x}_t \in \mathbb{R}^q, \boldsymbol{z}_t \in \mathbb{R}^p$. Given $M \ge 2$, denote the density of an $M$-component model that represents the conditional density function of $\{y_t\}_{t=1}^T$ given $\{\boldsymbol{x}_{t},\boldsymbol{z}_{t}\}_{t=1}^T$ as
where $\boldsymbol{\vartheta}_M=(\boldsymbol{\alpha}^\top,\boldsymbol\theta_1^\top,...,\boldsymbol\theta_M^\top,\gamma^\top)^\top\in\Theta_{\boldsymbol{\vartheta}_M}$, $\boldsymbol\alpha^\top :=(\alpha_1,...,\alpha_{M-1})$, $\alpha_M=1-\sum_{j=1}^{M-1}\alpha_j$ and
is the $j$-th component density function with $\mu_j \in \Theta_{\mu} \subset \mathbb{R}$ , $\sigma_j^2 \in \Theta_{\sigma} \subset \mathbb{R}_{++}$, $\boldsymbol\beta_j \in \Theta_{\boldsymbol\beta} \subset \mathbb{R}^{q}$, $\boldsymbol\gamma \in \Theta_{\boldsymbol\gamma} \subset \mathbb{R}^p$, and $\phi(t) = (2\pi)^{-1/2} \exp(-\frac{t^2}{2})$ is the standard normal probability density function. We collect the component-specific parameters into $\boldsymbol{\theta}_j := (\mu_j, \sigma_j^2, \boldsymbol{\beta}_j^\top)^\top \in \Theta_{\boldsymbol\theta}$, and the regression coefficient $\boldsymbol\gamma$ for a vector $\boldsymbol z$ is assumed to be common across components.
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},\boldsymbol{Z}^{\top}_{it} )^{\top} \}_{t=1}^T$ from a true $M_0$-component density $f_M(\boldsymbol{w};\boldsymbol\vartheta_{M_0}^*)$ defined in equation ((ref)) with $\boldsymbol\vartheta_{M_0}^*=((\boldsymbol{\alpha}^*)^\top,(\boldsymbol\theta_1^*)^\top,...,(\boldsymbol\theta_{M_0}^*)^\top,(\gamma^*)^\top)^\top$. The superscript $*$ signifies 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}^*$.}
Our goal is to test $$ H_0: M = M_0\ \text{ against }\ H_A: M = M_0 +1. $$
We begin by developing the PLRT to test the null hypothesis $H_0: M=1$ against the alternative hypothesis $H_1: M=2$. 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}_{it}^{\top}, \boldsymbol{Z}_{it}^{\top})^{\top} \}_{t=1}^T$, drawn from a true one-component density $f(\boldsymbol{w}; \boldsymbol{\gamma}, \boldsymbol{\theta})$ defined in equation ((ref)). Now consider a two-component mixture density function
where $\boldsymbol{\vartheta}_2 = (\alpha, \boldsymbol{\theta}_1^{\top}, \boldsymbol{\theta}_2^{\top}, \boldsymbol{\gamma}^{\top})^{\top} \in \Theta_{\boldsymbol{\vartheta}_2}$ and $\alpha$ is the mixing probability of the first component. The two-component model 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, the null hypothesis $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 for a standard asymptotic analysis fail in 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.
As discussed in the introduction, analyzing the asymptotic distribution of the LRTS for the cross-sectional normal mixture is challenging because of its undesirable mathematical properties 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 $f_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$. The presence of problems (i)--(iii) in panel normal mixture models with $T\geq 2$ is not well understood in the literature because, to the best of our knowledge, no studies have examined them so far.
Regarding problem (i), we note that the issue of the infinite Fisher information for testing $H_{02}$ also arises in the panel normal mixture model. For brevity, let us consider the case without $(\boldsymbol{X}, \boldsymbol{Z})$. The score for testing $H_{02}: \alpha=0$ takes the form \[ \left.\frac{\partial f_2(\boldsymbol{W} ; \mu_1,\sigma_1^2,\mu_2,\sigma_2^2)}{\partial \alpha} \right|_{\alpha=0,\mu_2=\mu^*,\sigma_2^2=\sigma^{*2}}= \frac{f(\boldsymbol{W};\mu_1,\sigma_1^2)}{f(\boldsymbol{W};\mu^*,\sigma^{*2} )}-1, \] where $f(\boldsymbol{W};\mu,\sigma^2)= \prod_{t=1}^T \phi((y_t-\mu)/\sigma)/\sigma$ and $\phi(\cdot)$ is the standard normal density function. When $\sigma_1^2> 2\sigma^{*2}$, $\mathbb{E}[\{f(\boldsymbol{W};\mu_1,\sigma_1^2)/f(\boldsymbol{W};\mu^*,\sigma^{*2} )-1\}^2]=\infty$. For more details, please refer to Proposition (ref). Because the infinite Fisher information causes difficulty in deriving the asymptotic distribution under $H_{02}$, this paper focuses on testing $H_{01}$. We define $ \Upsilon^*_1 := \{ (\alpha, \boldsymbol{\gamma}, \boldsymbol{\theta}_1, \boldsymbol{\theta}_2) \in \Theta_{\vartheta_2} : \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^* \text{ and } \boldsymbol{\gamma} = \boldsymbol{\gamma}^* \}$, which is the subspace of $\Theta_{\vartheta_2}$ that corresponds to $H_{01}$. Note that because we focus on $H_{01}$, our test may not have power against the local alternatives with $\alpha_n\rightarrow 0$. We analyze the asymptotic distribution of the PLRTS under the contiguous local alternatives in Section (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 f_2(\boldsymbol{W}_i;\boldsymbol{\vartheta}_{2}) - \sum_{i=1}^n\log f(\boldsymbol{W}_i;\boldsymbol{\gamma}^*,\boldsymbol{\theta}^*)\right\}, $$ where $f_2$ is the density of the two-component finite mixture distribution in ((ref)) with $M=2$ and $((\boldsymbol{\gamma}^*)^{\top},(\boldsymbol{\theta}^*)^{\top})^{\top}$ is the true parameter value under $H_0$. Let $\tilde{\boldsymbol{\vartheta}}_{2,n}$ be the maximum likelihood estimator for the two-component model, i.e., $\tilde{\boldsymbol{\vartheta}}_{2,n}=\arg\max_{\boldsymbol{\vartheta}_{2}\in \Theta_{\boldsymbol{\vartheta}_{2}}} LR_n^*(\boldsymbol{\vartheta}_{2})$.
To deal with unboundedness, we consider a penalized maximum likelihood estimator (PMLE) as in Chen2009a using the following penalty function:
This penalty function circumvents the problem of unbounded log likelihood by preventing a variance parameter estimate from nearing zero. The parameter $a_n$ is selected such that the penalty's impact becomes asymptotically negligible for the distribution of the PMLE. Refer to conditions C1--C3 in the proof of Proposition (ref).
Let \[ \hat{\boldsymbol{\vartheta}}_2 = \arg \max_{{\boldsymbol{\vartheta}}_2 \in \Theta_{\boldsymbol{\vartheta}_2} } \sum_{i=1}^n \log f_2(\boldsymbol{W}_i ; \boldsymbol{\vartheta}_2 )+\tilde p_n({\boldsymbol{\vartheta}}_2) \] denote the PMLE under the two-component model. Define a set of parameter values for the two-component density that generates the true one-component density by $\Theta^*_2:= \{ (\alpha,\boldsymbol{\gamma},\boldsymbol{\theta}_1,\boldsymbol{\theta}_2) \in \Theta_{\vartheta_2}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^* \text{ and } \boldsymbol{\gamma} = \boldsymbol{\gamma}^*; \alpha=1 \text{ and } \theta_1=\theta^*; \alpha=0 \text{ and } \theta_2=\theta^* \}$. $\boldsymbol{\theta}_1$ and $\boldsymbol{\theta}_2$ are component-specific parameters and $\boldsymbol{\gamma}$ is a parameter vector common across components. The following proposition establishes the consistency of the PMLE.
It should be noted that $f_2(\boldsymbol{w}; \boldsymbol{\vartheta}_2^*) = f(\boldsymbol{w}; \boldsymbol\gamma^*,\boldsymbol\theta^*)$ for any $\boldsymbol{\vartheta}_2^* \in \Theta_{2}^*$. Consequently, Proposition (ref) suggests that the PMLE $\hat{\boldsymbol{\vartheta}}_2$ converges in probability to a set of parameters for which the true density function $f(\boldsymbol{w}; \boldsymbol\gamma^*,\boldsymbol\theta^*)$ emerges within the space of two-component density functions.
For problem (iii), we show that in normal mixture models with panel data, the first-order derivative of $f_2(\boldsymbol{w}; \boldsymbol{\vartheta}_2)$ with respect to $\sigma_j^2$ is not linearly dependent with its second-order derivative with respect to $\mu_j$ (See Proposition (ref)(c)). Consequently, the panel mixture model ((ref)) with the 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 Proposition (ref)(a). In contrast, the strong identifiability does not hold for the cross-sectional normal mixture, and its convergence rate becomes as slow as $n^{-1/8}$ when the number of components is over-specified Kasahara2015a.
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}, (\boldsymbol\gamma^*)^{\top})^{\top}$ is linear dependent as
To deal with this linear dependency, we analyze the asymptotic distribution of the LRTS by developing a higher-order approximation for the log-likelihood function.
To extract the direction of the Fisher information matrix singularity, we adapt the reparameterization approach by Kasahara2012 and consider the following one-to-one reparameterization of $\boldsymbol{\theta}_1$ and $\boldsymbol{\theta}_2$ given $\alpha$:
where $\boldsymbol{\nu}$ and $\boldsymbol{\lambda}$ are both $(q+2) \times 1$ reparameterized parameter vectors with $\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}$. We also write $\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$ as $\boldsymbol{\theta} =(\theta_1,\theta_2,\theta_3,...,\theta_{q+2})^{\top}:= (\mu,\sigma^2,\beta_1,...\beta_q)^{\top}$ and $\boldsymbol{\lambda}=(\lambda_1,\lambda_2,\lambda_3,...,\lambda_{q+2})^{\top} := (\lambda_\mu,\lambda_\sigma,\lambda_{\beta_1},...,\lambda_{\beta_q})^{\top}$.
This reparameterization is essential for analyzing the asymptotic distribution of the PLRTS in light of the linear dependency in ((ref)). The reparameterized parameter $\boldsymbol{\lambda}$ captures a deviation from the one-component model, where its first-order derivatives of the log density are identically equal to zero under $H_0: M=1$. Consequently, this reparameterization facilitates the derivation of an approximate quadratic-form criterion function, which is based on the fourth-order Taylor series approximation of the log-likelihood function, to characterize the asymptotic distribution of the LRTS.
Define the space for reparameterized parameters as $$ \boldsymbol{\psi} := (\boldsymbol{\gamma}^{\top},\boldsymbol{\nu}^{\top},\boldsymbol{\lambda}^{\top})^{\top} \in \Theta_{\boldsymbol\psi}, $$ where $\Theta_{\boldsymbol{\psi}} = \{ \boldsymbol{\psi}: \boldsymbol{\gamma} \in \Theta_{\boldsymbol\gamma}, \boldsymbol{\nu} + ( 1 - \alpha) \boldsymbol{\lambda} \in \Theta_{\boldsymbol\theta}, \boldsymbol{\nu} - \alpha \boldsymbol{\lambda} \in \Theta_{\boldsymbol\theta}\}.$ Under the null hypothesis $H_{01}: \boldsymbol{\theta}_1 = \boldsymbol{\theta}_2 = \boldsymbol{\theta}^*$, we have $\boldsymbol{\lambda} = (0,\ldots, 0)^{\top}$ and $\boldsymbol{\nu} = \boldsymbol{\theta}^*$. We rewrite the reparameterized parameters under the null hypothesis as $(\boldsymbol{\psi}^{*})^{\top} = ((\boldsymbol{\gamma}^*)^{\top}, (\boldsymbol{\theta}^*)^{\top},0,\ldots,0)^{\top}$. Under the reparameterized parameter space, the density function and its logarithm are expressed as
Write $\boldsymbol{\psi}$ as $\boldsymbol{\psi} = (\boldsymbol{\eta}^{\top},\boldsymbol{\lambda}^{\top})^{\top}$ with $\boldsymbol{\eta} = (\boldsymbol{\gamma}^{\top},\boldsymbol{\nu}^{\top})^{\top}$, where $\boldsymbol{\eta}^* = ((\boldsymbol{\gamma}^*)^{\top},(\boldsymbol{\nu}^* )^{\top})^{\top}$ and $\boldsymbol{\lambda}^* = \boldsymbol{0}$. Denote the parameter spaces of $\boldsymbol\eta$ and $\boldsymbol\lambda$ by $\Theta_{\boldsymbol\eta}\subset \mathbb{R}^{p+q+2}$ and $\Theta_{\boldsymbol\lambda}\subset \mathbb{R}^{q+2}$, respectively.
Under this reparameterization, the first-order derivatives of the reparameterized log density with respect to the reparameterized parameters $\boldsymbol{\eta}$ are identical to those under the one-component model, and the first-order derivative with respect to $\boldsymbol{\lambda}$ is a zero vector:
With $\nabla_{\boldsymbol\lambda} l(\boldsymbol{w};\boldsymbol{\psi}^*,\alpha) = 0$, the Fisher information matrix is singular, and the standard quadratic approximation fails. Consequently, the information on $\boldsymbol{\lambda}$ is provided by the second-order derivative of $l(\boldsymbol{w};\boldsymbol{\psi},\alpha)$ with respect to $\boldsymbol{\lambda}$. We use the second-order derivative with respect to $\boldsymbol{\lambda}$ to identify $\boldsymbol{\lambda}$:
When $\alpha$ is bounded away from $0$ and $1$, the elements of $\nabla_{\boldsymbol{\lambda} \boldsymbol{\lambda}^{\top}} l(\boldsymbol{W};\boldsymbol{\psi}^*,\alpha)$ are mean-zero random variables.
Note that unlike the cross-sectional models analyzed by Kasahara2015a, there exists no collinearity between these first- and second-order derivatives for the panel models. This distinction is indeed important, as it highlights the differences in the asymptotic distribution of the LRTS for the panel models compared with the cross-sectional models. The absence of collinearity between the first- and second-order derivatives in the panel models leads to different convergence rates and asymptotic properties.
Let $f^*$ and $\nabla f^*$ denote $f(\boldsymbol{W};\boldsymbol{\gamma}^*,\boldsymbol{\theta}^*)$ and $\nabla f(\boldsymbol{W};\boldsymbol{\gamma}^*,\boldsymbol{\theta}^*)$, respectively. Define the vector $\boldsymbol{s}(\boldsymbol{W})$ as
The term $\widetilde{\nabla}_{\boldsymbol\theta\boldsymbol\theta^{\top}} f^*$ denotes the second-order derivatives of the density function $f^*$ with respect to the parameters $\boldsymbol{\theta}$. The coefficients $c_{jk}$ are used to adjust the scaling of these second-order derivatives. The function $\boldsymbol{s}(\boldsymbol{w})$ comprises the second-order derivatives of the log-likelihood function with respect to the reparameterized parameter $\boldsymbol\lambda$. This function, $\boldsymbol{s}_{\boldsymbol{\lambda} \boldsymbol\lambda}(\boldsymbol w)$, serves as a score function for identifying $\boldsymbol{\lambda}$. Consequently, $\boldsymbol{s}(\boldsymbol{w})$ is referred to as a score function. An explicit expression for the score function $\boldsymbol{s}(\boldsymbol{w})$ can be derived using Hermite polynomials, as elaborated in Appendix (ref).
Collect the relevant normalized reparameterized parameters and define $\boldsymbol{t}(\boldsymbol{\psi},\alpha)$ as
where $v(\boldsymbol{\lambda})$ is a vector of unique elements of $\boldsymbol\lambda\boldsymbol\lambda^{\top}$ given by
the length of which is $q_\lambda:=(q+2)(q+3)/2$.
Let $L_n(\boldsymbol{\psi},\alpha):= \sum_{i=1}^n l(\boldsymbol{W}_i;\boldsymbol{\psi}^*,\alpha)$ be the reparameterized log-likelihood function and define the normalized score vector \[ \boldsymbol S_n := n^{-1/2} \sum_{i=1}^n {\boldsymbol s}(\boldsymbol W_i). \] Then, taking the fourth-order Taylor expansion of $L_n(\boldsymbol{\psi},\alpha)$ around $(\boldsymbol\psi^*,\alpha)$, we may write $2\{L_n(\boldsymbol{\psi},\alpha)- L_n(\boldsymbol{\psi}^*,\alpha)\}$ as a quadratic function of $\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha)$ as
where $\boldsymbol {\mathcal{I}}_n$ is the negative of the sample Hessian defined in the proof of Proposition (ref) and $\boldsymbol G_n:=\boldsymbol{\mathcal{I}}_n^{-1}\boldsymbol S_n$. Let $\boldsymbol{\mathcal{I}}=\mathbb{E}[\boldsymbol s(\boldsymbol W)\boldsymbol s(\boldsymbol W)^{\top}]$.
The non-singularity of $\boldsymbol{\mathcal{I}}$ in Proposition (ref)(c) highlights the difference between the cross-sectional normal mixture and the panel data normal mixture models. In particular, as shown in equation ((ref)) in Appendix (ref), the first-order derivative of $f_2(\boldsymbol w;\boldsymbol{\vartheta}_2)$ with respect to $\sigma_j^2$ is linearly independent of its second-order derivative with respect to $\mu_j$ when $T\geq 2$, which ensures that the higher-order degeneracy of problem (iii) does not arise. Intuitively, the availability of repeated observations within each individual unit provides better identification, even for over-parameterized models, and reduces the degree of higher-order degeneracy.
The set of feasible values of $\sqrt{n} \boldsymbol{t}(\boldsymbol{\psi},\alpha)$ is given by the shifted and rescaled parameter space for $(\boldsymbol\eta,v(\boldsymbol\lambda))$ defined as $\Lambda_{n} := \sqrt{n} (\Theta_{\boldsymbol\eta}-\eta^*) \times \sqrt{n} \alpha(1-\alpha)v( \Theta_{\boldsymbol\lambda})$, where $v(A):= \{ t \in \mathbb{R}^{q_\lambda}: t = v(\lambda) \text{ for some } \lambda \in A \subset \mathbb{R}^{q+2}\}$. Because $\Lambda_n/\sqrt{n} $ is locally approximated by a cone $\Lambda:= \mathbb{R}^{p+q+2} \times v(\mathbb{R}^{q+2})$, we can apply Lemma 2 of Andrews1999 to approximate the distribution of the supremum of the right-hand side of ((ref)) as \[ \max_{\boldsymbol\psi\in \Theta_{\boldsymbol\psi}}2\{L_n(\boldsymbol{\psi},\alpha)- L_n(\boldsymbol{\psi}^*,\alpha)\} \overset{d}{\rightarrow} \boldsymbol G ^{\top} \boldsymbol{\mathcal{I}} \boldsymbol G-\inf_{\boldsymbol t\in \Lambda} (\boldsymbol t-\boldsymbol G)'\boldsymbol{\mathcal{I}} (\boldsymbol t-\boldsymbol G), \] where $\boldsymbol G = \boldsymbol{\mathcal{I}}^{-1} \boldsymbol S\sim N(0,\boldsymbol{\mathcal{I}}^{-1})$. This allows us to characterize the asymptotic distribution of the LRTS.
For each $\alpha\in (0,1)$, define the reparameterized PMLE as
with $\hat{\boldsymbol{\psi}} := (\hat{\boldsymbol{\gamma}}^{\top},\hat{\boldsymbol{\nu}}^{\top},\hat{\boldsymbol{\lambda}}^{\top})^{\top}$, where $\Theta_{\boldsymbol\psi}$ is defined as the space of $\boldsymbol{\psi}$ such that the $\boldsymbol{\vartheta}_2$ implied is in $\Theta_{\boldsymbol{\vartheta}}$ and $\sigma_j^2(\boldsymbol{\psi},\alpha)$ is the value of $\sigma_j$ implied by the value of $\boldsymbol{\psi}$ and $\alpha$ (e.g., $\sigma_1^2(\boldsymbol\psi,\alpha)=\nu_\sigma +(1-\alpha)\lambda_\sigma$).
Let $(\hat{\boldsymbol{\gamma}}_0,\hat{\boldsymbol{\theta}}_0)$ be the one-component MLE that maximizes the one-component likelihood function $L_{0,n}(\boldsymbol{\gamma},\boldsymbol{\theta}) := \sum_{i=1}^n \log f(\boldsymbol{W}_i;\boldsymbol{\gamma},\boldsymbol{\theta}) $. Define the LRTS and the PLRTS of testing $H_{01}$ with a small positivity constant $\epsilon$ on $\alpha$ as, respectively,
A hard bound is imposed on the values of $\alpha$ to avoid an issue of the infinite Fisher information for testing $H_{02}$. However, the LRTS may have reduced power if the true value of $\alpha$ does not satisfy the constraint $[\epsilon, 1- \epsilon]$ given an ad hoc constant $\epsilon>0$. For this reason, in Section (ref), we also develop the EM test, which does not impose a direct constraint on the value of $\alpha$.
With $\boldsymbol{s}(\boldsymbol{W})$ in ((ref)), partition $\boldsymbol{\mathcal{I}}=\mathbb{E}[\boldsymbol s(\boldsymbol W)\boldsymbol s(\boldsymbol W)^{\top}]$ and define
where $\boldsymbol{S}_{\boldsymbol\lambda,\boldsymbol\eta} \sim N(0,\boldsymbol{\mathcal{I}}_{\boldsymbol\lambda,\boldsymbol\eta} )$. Define a set that characterizes the feasible values of $\sqrt{n}\boldsymbol{t}_\lambda(\boldsymbol{\lambda},\alpha)$ when $n \to \infty$ by the cone
Define $\hat{\boldsymbol{t}}_{\lambda}$ as
where $\hat{\boldsymbol{t}}_{\lambda}$ is a projection of a random Gaussian random variable $\boldsymbol{G}_{\boldsymbol{\lambda}}$ on a cone $\Lambda_{\boldsymbol\lambda}$.
The following proposition establishes the asymptotic distribution of the LRTS or PLRTS under the null hypothesis $H_0: M=1$.
Proposition (ref)(a) implies that $\hat{\boldsymbol{\theta}}_j - \boldsymbol{\theta}^*=O_p(n^{-1/4})$ for $j=1,2$. The $n^{1/4}$ convergence rate is a consequence of the linear dependency in ((ref)), where the identification of the parameter $\boldsymbol{\theta}$ relies on the fourth-order Taylor approximation of the log-likelihood function. This rate is also the best convergence rate for an over-parameterized mixture under the strong identifiability condition chen95as. When we choose the penalty function so that $\sum_{j=1}^2 p_n(\sigma_j^2(\hat{\boldsymbol{\psi}},\alpha))=o_p(1)$ under the null hypothesis of $M=1$, $PLR_n$ has the same asymptotic null distribution as $LR_n$.
In this section, we build upon the analysis from the previous section and derive the asymptotic distribution of the PLRTS for testing the null hypothesis of $M_0$ components against an alternative of $(M_0+1)$ components, where $M_0 \ge 2$.
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},\boldsymbol{Z}^\top_{it} )^\top \}_{t=1}^T$ from an $M_0$-component density $f_{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}^*,\boldsymbol{\gamma}^*) \in \Theta_{\boldsymbol{\vartheta}_{M_0}}$ and $\alpha_{M_0}^*=1-\sum_{j=1}^{M-1}\alpha_j^*$.
Let the density of the $(M_0 + 1)$-component model be defined by
where $\boldsymbol{\vartheta}_{M_0+1} = (\boldsymbol{\theta}_1,\boldsymbol{\theta}_2,\ldots,\boldsymbol{\theta}_{M_0+1}, \alpha_1,\ldots,\alpha_{M_0},\boldsymbol{\gamma}) \in \Theta_{\vartheta_{M_0+1}}$ as defined in ((ref)). We assume $\mu_{1}^* < \mu_{2}^*, \ldots, < \mu_{M_0}^*$ in the true parameters for identification.
The $(M_0 + 1)$-component model ((ref)) gives rise to the true density ((ref)) in two cases: (i) two components have the same mixing parameter and (ii) one component has zero mixing proportion. Accordingly, we partition the null hypothesis of $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 $H_{02}: \alpha_h = 0$ for some $h=1,\ldots,M_0+1$.
We first analyze the infinite Fisher information problem for testing $H_{02}$. 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 \[
\] The score for testing $H_{0,2h}: \alpha_h = 0$ takes the form $\nabla_{\alpha_h} \log f_{M_0+1}(\boldsymbol W_i, \boldsymbol \vartheta_{M_0+1}) = [f(\boldsymbol W_i; \mu_h, \sigma_h^2 ) - f(\boldsymbol W_i; \mu_{M_0}^*, \sigma_{M_0}^{2*} )] / f_{M_0}(\boldsymbol W_i, \boldsymbol \vartheta_{M_0}^*)$. Because $(\mu_h,\sigma_h^2)$ is not identified when $\alpha_h = 0$, the Fisher information matrix of the LRTS for testing $H_{0,2h}: \alpha_h = 0$ depends on the supremum of the variance of $\nabla_{\alpha_h} \log f_{M_0 + 1}(\boldsymbol W_i; \boldsymbol \vartheta_{M_0+1})$ over $\boldsymbol \vartheta_{M_0+1} \in \Upsilon_{2h}^*$. The Fisher information is infinite unless there is an a priori restriction on the values of $\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 $\Upsilon^*_1 := \cup_{h=1}^{M_0} \Upsilon^*_{1h}$ corresponds to $H_{01}=\cup_{h=1}^{M_0}H_{0,1h}$.
Suppose that the null hypothesis of $M=M_0$ holds with the true density ((ref)). Because any parameter in $\Upsilon^{*}_{1} = \cup_{h=1}^{M_0} \Upsilon^*_{1h}$ can generate the true density $f_{M_0}(\boldsymbol{w}; \boldsymbol{\vartheta}^*_{M_0}) =\sum_{j=1}^{M_0} \alpha_0^{j*} f(\boldsymbol{w};\boldsymbol{\gamma}^*,\boldsymbol{\theta}_j^{*})$, we need to restrict the estimators under the $(M_0 + 1)$-component model to be in a neighborhood of $\Upsilon^*_{1h}$ to test $H_{0,1h}$.
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 $D_1^* = [\underline{\Theta}_{\mu}, \frac{\mu_1^* + \mu_2^* }{2}] \times \Theta_{\beta} \times \Theta_{\sigma^2}$, $D_h^* = [\frac{\mu_{h-1}^* + \mu_{h}^* }{2}, \frac{\mu_{h}^* + \mu_{h+1}^* }{2}] \times \Theta_{\beta} \times \Theta_{\sigma^2}$ for $h = 2,\ldots,M_0 - 1$, $D_{M_0 }^* = [\frac{\mu_{M_0-1}^* + \mu_{M_0}^* }{2}, \overline{\Theta}_{\mu}] \times \Theta_{\beta} \times \Theta_{\sigma^2}$. Then, $D_h^* \subset \Theta_{\theta}$ is a neighborhood containing $\theta_h^{*}$ but not $\theta_j^{*}$ for $j \neq h$. For $h=1,\ldots M_0$, given a small positive constant $\epsilon>0$, define a restricted parameter space $\boldsymbol{\Psi}_h^* \subset \Theta_{\vartheta_{M_0 + 1}}(\epsilon)$ as
Note that $\boldsymbol{\Psi}_h^* \cap \Upsilon_{1h}^* \neq \emptyset$ and $\boldsymbol{\Psi}_h^* \cap \Upsilon_{1l}^* = \emptyset$ if $h \neq l$, and $\cup_{h=1}^{M_0}\boldsymbol{\Psi}_h^*= \Theta_{\vartheta_{M_0 + 1}}(\epsilon)$.
Let $\hat{\boldsymbol{\Psi}}_h^* $ and $\hat{D}^*_h$ be consistent estimators of ${\boldsymbol{\Psi}}_h^* $ and ${D}^*_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 $\vartheta^{M_0+1}\in \hat{\boldsymbol{\Psi}}_h^*$.
For $h=1,2,..., M_0$, define the local PMLE that maximizes the log-likelihood function of the $(M_0+1)$-component model under the constraint that $\boldsymbol\vartheta_{M_0+1}\in \hat{\boldsymbol\Psi}_h^*$ in ((ref)) by \[ \hat{\boldsymbol{\vartheta}}_{M_0 + 1}^h = \operatorname*{arg\,max}_{\boldsymbol{\vartheta}_{M_0 + 1} \in \hat{\boldsymbol\Psi}_h^* } L_{M_0+1,n}(\boldsymbol{\vartheta}_{M_0 + 1}) + \tilde p_n(\boldsymbol{\vartheta}_{M_0+1 }), \] where $$ L_{M,n}(\boldsymbol{\vartheta}_{M}) : = \sum_{i=1}^n \log f_{M}(\boldsymbol{W}_i;\boldsymbol{\vartheta}_{M})\quad\text{and}\quad \tilde p_n(\boldsymbol\vartheta_{M}) :=\sum_{j=1}^{M} p_n (\sigma_j^2; \hat{\sigma}_{0,j}^2). $$ with
where $\hat{\sigma}_{0,j}^2$ is a root-$n$ consistent estimator of $\sigma_{0,j}^2$ from the $M_0$-component model under the null hypothesis. Because $\hat{\sigma}_j^2 - \sigma_{0,j}^2 = O_p(n^{-1/4})$ under the null hypothesis (cf. Proposition (ref)(a)), $p_n(\hat\sigma_j^2; \hat{\sigma}_{0,j}^2) =o_p(1)$ when $a_n$ is chosen to be $o(n^{1/4})$.
Under $H_{0}: M=M_0$, $\boldsymbol\Psi_h^*$ contains a set of parameters $\Upsilon_{1h}^*$ defined in ((ref)) such that $f_{M_0 + 1}(\boldsymbol{w}; \boldsymbol{\vartheta}_{M_0 + 1})$ is equal to $f_{M_0}(\boldsymbol{w}; \boldsymbol{\vartheta}_{M_0 }^*)$ for any $\boldsymbol{\vartheta}_{M_0 +1}\in\Upsilon_{1h}^*$ and is therefore the density function from which the data are generated. These penalized likelihood estimators are consistent.
Consider the local PLRTS for testing $H_{0,1h}: \boldsymbol\vartheta_h = \boldsymbol\vartheta_{h+1}$ defined by
The test utilizing the local PLRTS, denoted by $PLR^{M_0,h}_n$, possesses power solely against local alternatives within the restricted parameter space of $\boldsymbol{\Psi}_h^*$. To guarantee power against local alternatives over a wide range of directions, we consider the PLRTS characterized by the maximum of the local PLRTS for $h=1,...,M_0$, as defined by
Because $\Theta_{\vartheta_{M_0 + 1}}(\epsilon)= \cup_{h=1}^{M_0} \hat{\boldsymbol\Psi}_h^*$, $PLR_n(M_0)$ is identical to $\max_{\boldsymbol{\vartheta}_{M_0 + 1} \in \Theta_{\boldsymbol{\vartheta}_{M_0 + 1} }(\epsilon) } \{L_{M_0+1,n} (\boldsymbol{\vartheta}_{M_0 + 1})+ \tilde p_n(\boldsymbol{\vartheta}_{M_0+1 })\} - L_{M_0,n}(\hat{\boldsymbol{\vartheta}}_{M_0 })$.
To derive the asymptotic null distribution of $PLR_n(M_0)$, collect the score vector for testing $H_{0,1h}$ for $h = 1,\ldots, M_0$ into one vector as
where
with $\widetilde{ \nabla}_{\boldsymbol\theta_h \boldsymbol\theta^{\top}_h} f(\boldsymbol{W};\boldsymbol{\gamma}^*,\boldsymbol{\theta}^*_h) := (c_{11} \nabla_{\theta_{h1}\theta_{h1}} f^*,...,c_{(q+2)(q+2)}\nabla_{\theta_{h,q+2}\theta_{h,q+2}} f^*,c_{12}\nabla_{\theta_{h1}\theta_{h2}} f^*,...,c_{(q+1)(q+2)}\nabla_{\theta_{h,q+1}\theta_{h,q+2}} f^*)^{\top}$ for $\boldsymbol{\theta}_h:=(\theta_{h1},\theta_{h2},\theta_{h3},...,\theta_{h,q+2})^{\top}:=(\mu_h,\sigma_h^2,\beta_{h1},...,\beta_{hq})^{\top}$ and $c_{jk}=1/2$ for $j\neq k$ and $c_{jk}=1$ for $j=k$. Define
Then, the asymptotic distribution of the normalized score function is given by \[ \tilde{\boldsymbol S}_n: = \frac{1}{\sqrt{n}} \sum_{i=1}^n \tilde{\boldsymbol s}(\boldsymbol W_i) \overset{d}{\to} \tilde {\boldsymbol S} \sim N(\boldsymbol 0, \tilde{\boldsymbol{\mathcal{I}}}), \] where, in view of ((ref)), $ \tilde {\boldsymbol S} $ may be partitioned as $ \tilde {\boldsymbol S} =( \tilde {\boldsymbol S}_{\boldsymbol\eta}^{\top}, \tilde {\boldsymbol S}_{\boldsymbol\lambda\boldsymbol\lambda}^{\top})^{\top}$ with $n^{-1/2} \sum_{i=1}^n \tilde{\boldsymbol s}_{\boldsymbol\eta}(\boldsymbol W_i)\overset{d}{\to} \tilde {\boldsymbol S}_{\boldsymbol\eta}$ and $n^{-1/2} \sum_{i=1}^n \tilde{\boldsymbol s}_{\boldsymbol\lambda\boldsymbol\lambda}(\boldsymbol W_i)\overset{d}{\to} \tilde {\boldsymbol S}_{\boldsymbol\lambda\boldsymbol\lambda}$ .
Let $\boldsymbol{\tilde{S}}_{\boldsymbol\lambda,\boldsymbol\eta} := (\boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^1, \ldots, \boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^{M_0} )^\top:= \boldsymbol{\tilde{S}}_{\boldsymbol\lambda\boldsymbol\lambda} - \tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol\lambda\boldsymbol\eta} \tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol\eta} ^{-1} \boldsymbol{\tilde{S}}_{\boldsymbol\eta} \sim N(0, \tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol\lambda,\boldsymbol\eta})$ be a $\mathbb{R}^{M_0 (q+2)(q+1)/2}$-valued random vector. For $h=1,2,...,M_0$, define $\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol\lambda,\boldsymbol\eta}^h := \mathbb{E}[\boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^h ( \boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^h)^\top]$ and $\boldsymbol{{G}}_{\boldsymbol\lambda,\boldsymbol\eta}^h := ({\boldsymbol{\mathcal{I}}}_{\boldsymbol\lambda,\boldsymbol\eta}^h)^{-1} \boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^h$.
Define $\hat{\boldsymbol t}^h_{\boldsymbol\lambda} $ analogously to $\hat{\boldsymbol t}_{\boldsymbol\lambda} $ as
The local quadratic-form approximation of the log-likelihood function $LR^{M_0,h}_n$ around $\Upsilon^*_{1h} \subset \Theta_{\vartheta_{M_0 + 1}}$ has an identical structure to the approximation that we derive in Section (ref) in testing $H_{01}$ in the test of homogeneity. Consequently, we can show that $PLR^{M_0,h}_n\overset{d}{\to} (\hat{\boldsymbol{t}}^h_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^h_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^h_{\boldsymbol\lambda}$. Then, given ((ref)), the asymptotic null distribution of the PLRTS for testing $H_{01}$ is given by the maximum over $(\hat{\boldsymbol{t}}^h_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^h_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^h_{\boldsymbol\lambda}$s for $h=1,2,..., M_0$.
The asymptotic null distribution of $PLR_n({M_0})$ is non-standard, but it is straightforward to simulate the random variable from the asymptotic null distribution using the estimates. Specifically, we simulate a draw of $\boldsymbol{\tilde{S}}_{\boldsymbol\lambda,\boldsymbol\eta}= (\boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^1, \ldots, \boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^{M_0} )^\top$ from $N(0, \hat{\tilde{\boldsymbol{\mathcal{I}}}}_{\boldsymbol\lambda,\boldsymbol\eta})$, where $\hat{\tilde{\boldsymbol{\mathcal{I}}}}_{\boldsymbol\lambda,\boldsymbol\eta}$ is a sample analogue estimator of ${\tilde{\boldsymbol{\mathcal{I}}}}_{\boldsymbol\lambda,\boldsymbol\eta}$. Then, compute $\boldsymbol{{G}}_{\boldsymbol\lambda,\boldsymbol\eta}^h = (\hat{\boldsymbol{\mathcal{I}}}_{\boldsymbol\lambda,\boldsymbol\eta}^h)^{-1} \boldsymbol{{S}}_{\boldsymbol\lambda,\boldsymbol\eta}^h$ and obtain $\hat{\boldsymbol t}^h_{\boldsymbol\lambda}$ analogously to ((ref)) using an estimator of $\boldsymbol{\mathcal{I}}_{\boldsymbol\lambda,\boldsymbol\eta}^h$ for $h=1,...,M_0$, and a simulated random draw is computed as $ \max \{(\hat{\boldsymbol{t}}^1_{\boldsymbol\lambda})^\top \hat{\boldsymbol{\mathcal{I}}}^1_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^1_{\boldsymbol\lambda} , \ldots, (\hat{\boldsymbol{t}}^{M_0}_{\boldsymbol\lambda})^\top \hat{\boldsymbol{\mathcal{I}}}^{M_0}_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^{M_0}_{\boldsymbol\lambda} \}$. Appendixes (ref) and (ref) present an expression for the score functions using Hermit polynomials.
This section develops an EM test used for testing the hypothesis $H_0: M=M_0$ against the alternative hypothesis $H_A: M = M_0 + 1$. A key limitation of the PLRT, as discussed in the previous section, is that the computation of mixing probabilities, denoted as $\alpha_j$, is subject to a hard constraint, which is dictated by an arbitrary choice of bounds. The EM test, in contrast, circumvents the need to impose an explicit constraint on the $\alpha_j$ values. It achieves this by performing a limited number of EM steps, starting from a predetermined set of $\alpha_j$ values. The EM test approach offers certain advantages, including computational simplicity and less stringent assumptions.
Let $\mathcal{T}$ be a finite set of numbers in $(0,0.5]$ with $0.5\in \mathcal{T}$, let $p(\tau)\leq 0$ be a penalty term that is continuous in $\tau$, $p(0.5)=0$, and let $p(\tau)\rightarrow -\infty$ as $ \tau$ goes to 0. Specifically, we choose $$ p(\tau) := \log(2\min\{\tau,1-\tau\}).$$
For each $\tau_0 \in \mathcal{T}$, let $\tau^{(1)} (\tau_0) = \tau_0$, and define the restricted PMLE by \[ \boldsymbol{\vartheta}_{M_0 + 1}^{h(1)}(\tau_0) = \underset{\boldsymbol{\vartheta}_{M_0 + 1} \in \Theta_{\boldsymbol{\vartheta}_{M_0 +1}}^h(\tau) }{\arg\max} {PL}_n(\boldsymbol{\vartheta}_{M_0 + 1},\tau_0), \] where $\Theta_{\boldsymbol{\vartheta}_{M_0 +1}}^h(\tau_0) := \{ \boldsymbol\theta \in \hat{\boldsymbol{\Psi}}_h : \alpha_h / (\alpha_h + \alpha_{h+1} ) = \tau_0 \}$ and $$ {PL}_n(\boldsymbol{\vartheta}_{M_0 + 1},\tau):= L_{M_0+1,n}(\boldsymbol{\vartheta}_{M_0 + 1}) + \tilde p_n(\boldsymbol{\vartheta}_{M_0+1}) +p(\tau). $$
Starting from $(\boldsymbol{\vartheta}_{M_0 + 1}^{h(1)}(\tau_0),\tau^{h(1)}(\tau_0 ) )$ with $\tau^{h(1)}(\tau_0 )=\tau_0$, update $\boldsymbol{\vartheta}_{M_0 + 1}^{h(k)}(\tau_0) $ and $\tau^{h(k)}(\tau_0)$ by the following generalized EM algorithm. Denote the estimators after the $k$-th round of EM algorithm iteration by $\vartheta_{M_0 + 1}^{h(k)}$ and $\tau^{h(k)}$. In the E-step, for $i=1,\ldots,N$ and $j = 1,\ldots, M_0 + 1 $, compute the weight for observation $i$ and type $j$ as
where, for brevity, we drop the superscript $h$ and its dependency on $\tau_0$ from the notations, such as in $w_{ij}^{h(k)}(\tau_0)$.
In the M-step, we update $\boldsymbol{\alpha}$ and $\tau$ by
We also update $\boldsymbol{\theta}_j$ and $\boldsymbol\gamma$ as
where $\tilde{\boldsymbol{x}}_{it} = (1,\boldsymbol{x}_{it}^\top)^\top$. In the updating procedure, $\boldsymbol{\vartheta}_{M_0 + 1}^{h(k+1)}(\tau_0)$ is not restricted to be in $\hat{\boldsymbol{\Psi}}_h^*$.
For each $\tau_0 \in \mathcal{T}$ and each step $k$, define
With a predetermined finite number $K$, define the local EM test statistic by taking the maximum of $M_n^{h(k)}(\tau_0)$ across different values of $\tau_0$ as
The test statistic $EM_n^{h}$ tests $H_{0,1h}: \boldsymbol\theta_h=\boldsymbol\theta_{h+1}$ and has a power against the local alternative that splits the $h$-th component of the null $M_0$-component model into two different components. To achieve power against a wide range of local alternatives, we consider the EM test statistic that takes the maximum of $M_0$ local EM test statistics:
Therefore, the asymptotic null distribution of the EM test statistic $EM_n(M_0)$ is the same as that of the PLRTS.
We derive the asymptotic distribution of the PLRTS and EM test statistic under local alternatives. For brevity, we focus on testing $H_0: M=1$ against $H_A: M=2$. Consider the following local alternative to the homogeneous model $f(\boldsymbol w;\boldsymbol\gamma^*,\boldsymbol\theta^*)$ with $\boldsymbol\theta^*=(\mu^*,\sigma^{*2},(\boldsymbol{\beta}^*)^{\top})^{\top}$. For brevity, we omit the common parameter $\boldsymbol \gamma$ in this section. In a reparameterized parameter, $\boldsymbol\psi^* = ((\boldsymbol\nu^*)^{\top}, (\boldsymbol{\lambda}^*)^{\top})^{\top}$. For $\alpha^*\in (0,1)$ and a local parameter $\boldsymbol h=(\boldsymbol {h}_{\boldsymbol\nu}^{\top},\boldsymbol {h}_{\boldsymbol\lambda}^{\top})^{\top}$ with $\boldsymbol h_{\boldsymbol \lambda}\in v(\boldsymbol\Theta_{\boldsymbol\lambda})$, we consider a sequence of contiguous local alternatives $(\alpha_n,\boldsymbol{\psi}_n^{\top})^{\top} = (\alpha_n,\boldsymbol\nu_n^{\top},\boldsymbol\lambda_n^{\top})\in \boldsymbol{\Theta}_\alpha\times \boldsymbol{\Theta}_{\boldsymbol\nu}\times\boldsymbol{\Theta}_{\boldsymbol\lambda}$ such that, with $\boldsymbol{t}_{\boldsymbol\lambda}(\boldsymbol\lambda,\alpha)$ given by ((ref)),
Equivalently, the non-reparameterized contiguous local alternatives are given by
for $\boldsymbol\nu_n = \boldsymbol \nu^* + n^{-1/2} \boldsymbol{h}_{\nu}$ and $\boldsymbol\lambda_n= (\lambda_{1,n},\lambda_{2,n},....,\lambda_{q+2,n})^{\top}$ with
where $\boldsymbol{h}_{\boldsymbol\lambda}=(h_{\lambda,1}^2,...,h_{\lambda,q+2}^2,h_{\lambda,1}h_{\lambda,2},....,h_{\lambda,q+1}h_{\lambda,q+2})^{\top}$. The local alternatives are of order $n^{1/4}$ rather than $n^{1/2}$. See the discussion following Proposition (ref).
The following proposition provides the asymptotic distribution of the PLRT and EM test statistics under contiguous local alternatives.
Importantly, a set of contiguous local alternatives considered in ((ref)) excludes a sequence such that $\alpha_n\rightarrow 0$ or $1$.
To estimate the number of components, we sequentially test $H_0: M=r$ against $H_1: M=r+1$ starting from $r=1$, and then $r=2,\ldots,\bar M$, where $\bar M$ is the upper bound for the number of components, which is assumed to be 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, \bar M$, let $c^{M}_{1-q_n}$ denote the $100(1-q_n)$ percentile of the cumulative distribution function of a random variable $ \max \{(\hat{\boldsymbol{t}}^1_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^1_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^1_{\boldsymbol\lambda} , \ldots, (\hat{\boldsymbol{t}}^{M}_{\boldsymbol\lambda})^\top \boldsymbol{\mathcal{I}}^{M}_{\boldsymbol\lambda,\boldsymbol\eta} \hat{\boldsymbol{t}}^{M}_{\boldsymbol\lambda} \}$ for $M=M_0$ in Propositions (ref) and (ref). Let $\hat 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 estimators $\hat M_{\text{PLR}}$ and $\hat M_{\text{EM}}$ depend on the choice of the significance level $q_n$. The following proposition states that $\hat M_{\text{PLR}}$ and $\hat M_{\text{EM}}$ converge to $M_0$ in probability as $n \rightarrow \infty$ when $-n^{-1}\ln q_n=o(1)$ and $q_n=o(1)$.
Let $Q_n^M( \boldsymbol{\vartheta}_M):=n^{-1} \sum_{i=1}^n \ln f_M(\boldsymbol{w}_i; \boldsymbol{\vartheta}_M)$ and $Q^M( \boldsymbol{\vartheta}_M):= \mathbb{E}[\ln f_M(\boldsymbol{w}_i; \boldsymbol{\vartheta}_M)]$, where $f_M(\boldsymbol{w}_i; \boldsymbol{\vartheta}_M)$ is defined in ((ref)) for $M=1,...,\bar M$.
Assumptions (ref)(a)--(e) ensure the consistency and asymptotic normality of $\hat{\boldsymbol{\vartheta}}_M$, 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$.
In this section, we examine the finite sample performance of the EM test and the PLRT by simulation. We test $H_0 : M = M_0 $ against $H_1 : M = M_0 + 1$ for the model with $M_0=2$ and $3$.
We develop a data-dependent empirical formula for $a_n$ by selecting a formula that ensures that the empirical rejection probabilities match the nominal size (5%) across various null models and sample sizes, as reported in Table (ref) in Appendix D. Specifically, for the model without conditioning variables, we derive the following data-dependent empirical formula for testing the null hypotheses of $M_0=1, 2, 3, 4$:
where $\omega(\boldsymbol{\vartheta}_{M_0};M_0)$ is the misclassification probability as defined in Melnykov2010 for each of the null models. The parameters $\hat{\rho}_1^{M_0}$, $\hat{\rho}_2^{M_0}$, $\hat{\rho}_3^{M_0}$, $\hat{\rho}_4^{M_0}$, and $\hat{\rho}_5^{M_0}$ are chosen as follows. Across different null models, sample sizes, and various candidate values of $a_n$, we estimate the empirical rejection probabilities at the $5\%$ significance level by simulations and denote them by $\hat s$. For example, when testing $H_0: M_0=2$, we repeatedly simulate the 500 datasets under each of the $48$ null model parameters and sample sizes $(N,T,\alpha,\mu,\sigma) \in \{100,500\} \times \{ 2,5,10\} \times \{(0.5,0.5),(0.2,0.8)\} \times \{ (-1,1),(-0.5,0.5), (-0.5,0.8)\} \times \{ (1,1), (1.5,0.75), (0.8,1.2)\}$ and test the null hypothesis of $H_0: M_0=2$ by the EM test using one of the six values of $a_n\in \{0.01, 0.05,0.1,0.2, 0.3, 0.4\}$. For each of the $108 \times 6= 648$ combinations of the parameter values, sample sizes, and $a_n$ values, let $\hat s$ denote the fraction of simulated datasets that lead to the rejection of the null hypothesis at the $5\%$ significance level. Using these $648$ observations of $\{\hat{s},N,T,\omega(\boldsymbol{\vartheta}_{2};2), a_n\}$, we run the following regression: \[
\] where $\hat{\rho}_1^{M_0}$, $\hat{\rho}_2^{M_0}$, $\hat{\rho}_3^{M_0}$, $\hat{\rho}_4^{M_0}$, and $\hat{\rho}_5^{M_0}$ in ((ref)) denote the corresponding estimates. Table (ref) in the Appendix reports the estimates. Note that the data-dependent formula ((ref)) is obtained by setting $\hat s=0.05$ and solving for $a_n$ in the above equation.
For the model with conditioning variables, we find that the value of $a_n$ that gives accurate Type I errors is sensitive to the dimension of covariates, and developing a data-dependent empirical formula for $a_n$ is difficult. Consequently, we choose a constant value of $a_n$ that depends only on the number of components $M_0=1, 2, 3,$ and $4$ as follows: $a_n= 0.1617 \text{ if } M_0 =1 ; a_n= 0.0025 \text{ if } M_0 = 2; a_n= 0.0567 \text{ if } M_0 = 3; a_n= 0.4858 \text{ if } M_0 = 4; \text{ and } a_n= 0.5 \text{ if } M_0 \geq 5.$ These penalty terms for the regression with covariates are chosen by averaging the predictions of the penalty function for the null parameters used in the simulations. For example, the penalty term for $M_0 = 2$ is chosen by generating $a_n$ using the formula for all of the combinations of $(N, T,\alpha, \mu,\sigma)$ in Table (ref) for $M_0 = 2$ and taking the average across the predicted values of $\hat a_n$. For $M_0 \ge 5$, we use the parametric bootstrap method to obtain the critical values for our empirical application, where we set $a_n = 0.5$.
Table (ref) displays the simulated Type I error rates for the EM test when we examine the null hypothesis $H_0: M=2$ against the alternative hypothesis $H_1: M=3$. A total of 2,000 repetitions are used for the asymptotic distribution, and 1,000 repetitions are used for the bootstrap distribution. Moreover, the PLRT with simulated critical values is considered.
The table presents the results for four distinct null models, as explained in the table's footnote. Utilizing the asymptotic distribution, the EM test sizes generally approximate the nominal 5% level. Nonetheless, the test may be undersized in instances where $T \ge 5$. Furthermore, the test size is larger when the mixing proportions are equal ($\boldsymbol{\alpha} = (0.5,0.5)$) than when they are unequal ($\boldsymbol{\alpha} = (0.2, 0.8)$). The bootstrapped EM test demonstrates satisfactory performance.
For the PLRT, 2,000 repetitions are conducted, and the results are reported for cases where a constraint is applied to $\alpha_j \in [\epsilon, 1 - \epsilon]$ with $\epsilon=0.1$. The value of $a_n$ for the PLRT is chosen to be 10 times larger than its value for the EM test. The findings suggest that the PLRT is slightly oversized.
Table (ref) reports the rejection frequency of testing $H_0: M_0 = 2$ under 12 alternative three-component mixture models, as elaborated in the table's footnote. For both the EM test and the PLRT, the test power is greater when the distances between $\mu_j$s are larger and equal, such as $(\mu_1,\mu_2,\mu_3)=(-1,0,1)$ or $(-1.5,0,1.5)$, as opposed to unbalanced distances such as $(-1,0,2)$ or $(-0.5,0,1.5)$. The power is also improved when the mixture probabilities are equal ($\boldsymbol{\alpha} = (1/3,1/3,1/3)$) relative to when they are unequal ($\boldsymbol{\alpha} = (1/4,1/2,1/4)$). The power increases with both the time dimension $T$ and the cross-sectional sample size $N$. Reflecting a larger actual rejection frequency of the PLRT under $H_0: M_0 = 2$ in Table (ref), the power of the PLRT is often higher than that of the EM test, although the EM test sometimes has higher power, especially when the mixing probabilities are unequal.
Table (ref) displays the simulated Type I error rates of the EM test using the asymptotic distribution for testing $H_0: M_0=3$ against $H_1: M_0=4$. Six null models are considered with varying $(\alpha_1,\alpha_2, \alpha_3)$ and $(\mu_1,\mu_2, \mu_3)$ values. The EM test generally yields accurate Type I errors.
The Type I error rates of the EM test with conditioning variables under the null $M_0 = 2$ are examined using 500 repetitions. The results presented in Table (ref) indicate a slightly oversized test for small samples with $(N, T) = (200, 2)$, but overall, the finite sample properties are satisfactory.
In our empirical application examining production function heterogeneity in Japan and Chile, we find evidence that the number of components is frequently greater than 5 when we sequentially apply our EM test to estimate the number of components. We also investigate the performance of the SHT using the EM test in comparison with the AIC and the BIC when the data are generated from a five-component model in a realistic setting. Specifically, we simulate 100 datasets from the estimated five-component model of the Chilean textile industry in our empirical application and apply these three methods to select the number of components in each of the 100 datasets. Here, we apply the EM test at the 5% significance level to sequentially test the null hypothesis $H_0: M=M_0$ for $M_0=1,2,...,7$, and we determine the number of components to be $M_0$ when we fail to reject $H_0: M=M_0$, as in ((ref)).
Table (ref) presents the frequencies at which the three methods select the number of components in this simulation. The table demonstrates that the proposed SHT selects the correct number of components 72% of the time, while it underestimates the true number of components 25% of the time. Conversely, the AIC overestimates the number of components 86% of the time, and the BIC underestimates the number of components by selecting a four-component model 41% of the time and accurately estimates the number of components 58% of the time. Overall, in this simulation, our proposed SHT approach outperforms both the AIC and the BIC.
In this section, we conduct an empirical application of our proposed test for the number of components in a finite mixture production function model, the identification of which is analyzed in Kasahara2022esri. Specifically, we estimate the number of types of input elasticities in production functions using panel data from Japanese publicly traded firms in the machinery industry and data from Chilean manufacturing firms.
Consider the input and output panel data of $n$ firms over $T$ years, $\{\{Y_{it}, V_{it}, L_{it},$ and $K_{it}\}_{t=1}^T\}_{i=1}^N$, where $Y_{it}$, $V_{it}$, $L_{it}$, and $K_{it}$ represent the output, intermediate input, labor, and capital of firm $i$ in year $t$, respectively. We denote the logarithms of the corresponding variables by lowercase letters as $(y_{it}, v_{it}, l_{it}, k_{it})$, with, for example, $y_{it} = \log(Y_{it})$.
We use a finite mixture specification to capture the unobserved heterogeneity in a firm's input elasticities. We are interested in testing the number of production technology types. Assume that there are $M$ discrete types of production technologies and define the latent random variable $D_i \in \{ 1,2,\ldots, M\}$ to represent the production technology type of firm $i$. If $D_i = j$, then firm $i$ is of type $j$. The population proportion of type $j$ is denoted by $\alpha_j=\Pr(D=j)$. The production function for type $j$ is Cobb--Douglas, and the output is related to inputs as
with
where $\gamma_t^j$ represents the aggregate productivity shock of type $j$ in year $t$, $\omega_{it}$ is the serially correlated productivity shock, and $\epsilon_{it}$ is the idiosyncratic productivity shock.
We assume that an intermediate input $V_{it}$ is flexibly chosen by firm $i$ after observing the aggregate shock $\gamma_t^j$ and the serially correlated productivity shock $\omega_{it}$. The variable $\epsilon_{it}$ represents a mean-zero i.i.d. random variable, the realization of which is unknown when the intermediate input $V$ is selected. Denote the information available to a firm for making decisions on $V_{it}$ by $\mathcal{I}_{it}$. Denote the information available to a firm for making decisions on $V_{it}$ by $\mathcal{I}_{it}$.
To identify the intermediate input elasticity of the production function, we introduce the following assumptions (cf. Kasahara2022esri).
In Assumption (ref)(a), each firm's production function belongs to one of the $M$ types. Assumption (ref)(b) assumes that the idiosyncratic productivity shock follows a normal distribution. Assumption (ref)(c) assumes that both the aggregate shock $\gamma_t^j$ and the serially correlated productivity shock $\omega_{it}$ are observed when intermediate inputs are chosen but idiosyncratic productivity shocks are unknown. Assumption (ref) states that firms observe input and output prices when deciding on $V_{it}$. Assumption (ref) assumes that $V_{it}$ is chosen to maximize the current expected period profit conditional on the value of $(K_{it},L_{it})$.\footnote{We are agnostic about the timing of choosing $K_{it}$ and $L_{it}$ as long as they are either determined before $V_{it}$ or simultaneously chosen with $V_{it}$. It is reasonable to assume that capital input $K_{it}$ is determined before the value of $V_{it}$ is chosen. However, labor input $L_{it}$ may be flexibly chosen simultaneously with $V_{it}$ after $\gamma_t^j$ and $\omega_{it}$ are observed. Even when labor input is simultaneously chosen with intermediate input, equation ((ref)) and the corresponding first-order condition characterize the intermediate input choice once we interpret $L_{it}$ in ((ref)) as the optimal value chosen by firm $i$, as discussed in Ackerberg2015.}
Given the above Assumptions (ref), (ref), and (ref), we derive an empirical specification based on the first-order condition of the profit maximization problem ((ref)), following the idea developed by Gandhi2020 and extending it to a finite mixture production function modeled by Kasahara2022esri. Note that $E[\exp(\epsilon_{it})|D_i=j]=\exp( \sigma_j^2/2)$ for $\epsilon_{it}\sim N(0,\sigma_j^2)$. Then, because $\delta_{v,j}= \frac{\partial F_i^{j}(V_{it},K_{it},L_{it})/\partial V_{it}} {{F_i^{j}(V_{it},K_{it},L_{it})}/{V_{it}}}$ for the Cobb--Douglas production function, the first-order condition with respect to $V_{it}$ in ((ref)) together with the production function ((ref)) implies that
where \[ s_{it}:=\log\left(\frac{P_{V,t} V_{it}}{P_{Y,t} Y_{it}}\right) \] is the logarithm of the ratio of the intermediate input cost to revenue.
Collect the observed data as $\boldsymbol{W}_i = \{ s_{it},\log K_{it}\}_{t=1}^T$. Let $\mu_j = \log \delta_{v,j} + \frac{1}{2} \sigma_j^2$ and define a type-specific parameter to be $\boldsymbol{\theta}_j = (\mu_j,\sigma_j)$, where $\delta_{v,j}$ can be identified from $\boldsymbol{\theta}_j$ as $\delta_{v,j}=\exp(\mu_j-\sigma_j^2/2)$. Collect the parameters of each type and the mixing probability as $\boldsymbol{\vartheta}_M = (\alpha_1,\ldots,\alpha_{M-1}, \boldsymbol{\theta}_1^{\top},\ldots,\boldsymbol{\theta}_M^{\top})^{\top}$. Recall that $\epsilon_{it} \overset{iid}{\sim} N(0, \sigma_j^2)$ over $i$ and $t$ conditional on the technology type $D_i=j$. Then, from ((ref)), we can write the density function of $s_{i1},..., s_{iT}$ as a mixture of type-specific likelihood density similar to the density function in equation ((ref)):
The PMLE is defined as \[ \hat{\boldsymbol{\vartheta}}_{M} = \arg\max_{\boldsymbol{\vartheta}_M} \sum_{i=1}^{n} \log f_{M}(\boldsymbol{W}_i;\boldsymbol{\vartheta}_{M}) +\tilde p_n(\boldsymbol{\vartheta}_{M}). \]
As an alternative specification, we allow the elasticity of output for intermediate input to be a function of $\log K_{it}$ as $\log\delta_{v,j}=\beta_{0,j} + \beta_{k,j} \log K_{it}$. This results in the logarithm of the ratio of intermediate input cost to revenue being linearly related to $\log K_{it}$ as $ s_{it} = \mu_j +\beta_{k,j} \log K_{it} - \epsilon_{it} $ for $D_i=j$ with $\mu_j = \beta_{0,j}+ \frac{1}{2}\sigma_{j}^2$. In this case, the conditional density function of $\{ s_{it}\}_{t=1}^T$ given $\{\log K_{it}\}_{t=1}^T$ is
In addition, we consider a specification in which we include not only $\log K_{it}$ but also $\log L_{it}$ as a regressor:
We apply the EM test to two producer-level datasets to determine the number of production technology types. We use the production data from Japanese publicly traded firms from 2003 to 2007 and Chilean manufacturing plants from 1992 to 1996.\footnote{Please refer to kasahara21 and KASAHARA2008 for the details of the datasets of the Japanese publicly traded firms and the Chilean manufacturing plants, respectively. } We clean the data and use the firms/plants with continuous data entry for five years to ensure that we have balanced panel data. We focus on the three largest industries in terms of the number of firms and plants for each country (chemical, machine, and electronics for Japan and food products, fabricated metal products, and textiles for Chile). Table (ref) presents the summary statistics for the revenue share of intermediate materials and the log of gross output in these industries. The within-industry standard deviations of the revenue share of intermediate materials are substantial across all industries, suggesting that the intermediate input elasticities differ across firms within the narrowly defined industries.
To determine the number of components, we test the null hypothesis $H_0: M=M_0$ against $H_1: M=M_0+1$ by applying the EM test at the 5% significance level sequentially for $M_0 = 1,\ldots,5$. If we fail to reject the null hypothesis at a certain $M_0 = M$, then we conclude that there are $M$ types of intermediate input elasticities. We consider both the models without conditioning variables ((ref)) and the models with conditioning variables ((ref))--((ref)).
Tables (ref) and (ref) report the results of the EM test for the model without conditioning variables ((ref)) for the Japanese and the Chilean industries with a panel length of $T=3,4,5$ and a null model of $M=1,...,5$. For all industries in both countries and all panel lengths, we reject the null hypothesis of $H_0: M=M_0$ for all $M_0=1,2,3,4,$ and $5$ at the 5% significance level, which indicates that there are at least five types of intermediate input elasticities. This result reflects the considerable and persistent heterogeneity in the revenue share of intermediate materials across firms or plants, providing strong evidence for substantial heterogeneity in intermediate input elasticities across firms' production functions among Japanese and Chilean producers. Our findings serve as a caution against the conventional empirical practice of estimating the Cobb--Douglas production function, which assumes that elasticity parameters are common across firms. Given the strong evidence of heterogeneity in the production function coefficients, incorporating heterogeneity in production function coefficients in empirical applications is warranted and should be encouraged.
One possible reason for the estimated number of technology types being greater than 5 is that the assumption of the Cobb--Douglas production function may be too restrictive. When the production function is not Cobb--Douglas, the revenue share of intermediate materials generally depends on the value of production inputs Gandhi2020. For this reason, we test the number of technology types when the revenue share of intermediate materials depends on the values of capital input and labor input by estimating models ((ref))--((ref)).
Table (ref) presents the results of the SHT and the BIC when we estimate the mixture regression model with $\log K_{it}$ in ((ref)) using data with a panel length of $T=3$. For the Japanese chemical, electronics, and machinery industries, the SHT suggests that the data are generated from seven- to nine-component models; concurrently, the BIC selects models with at least 10 components. For the Chilean food industry, the SHT indicates a 10-component model, while the BIC chooses an eight-component model. In contrast, the SHT and the BIC respectively select models with seven and six components for the Chilean fabricated metal products industry and the Chilean textile industry.
Table (ref) reports the results for the model that includes both $\log K_{it}$ and $\log L_{it}$ as regressors. Across six industries, the SHT and the BIC in Table (ref) both select models with at least five components, providing evidence for substantial heterogeneity in production technology across firms and plants. Comparing the results of Table (ref) with those of Table (ref), the selected number of components for the model with $\log K_{it}$ and $\log L_{it}$ is smaller than that for the model with only $\log K$. This suggests that the number of components may be overestimated if we do not consider a sufficiently flexible production function specification by excluding some regressors.
The selection of the number of components in a finite normal mixture panel regression model is a crucial practical issue that must be addressed with care. Arbitrary choice of the number of components can result in biased estimates and invalid inferences, and can reduce the credibility of the final outcomes. To tackle this issue, this study proposes the PLRT and an EM test and derives their asymptotic distribution for the null hypothesis of a model with $M_0$ components against the alternative hypothesis with $(M_0 + 1)$ components. We also develop a procedure to consistently select the number of components by sequentially applying the PLRT and EM tests. Through a simulation exercise, we demonstrate that the proposed SHT procedure exhibits good performance in finite samples.
As an empirical application, we estimate the number of production technology types using producer-level panel data from Japan and Chile. We find that most industries in our dataset exhibit a level of heterogeneity that requires a five-or-more-component mixture model when using the Cobb--Douglas production specification or a specification in which the elasticity of inputs depends on capital and labor input linearly. This provides strong evidence of the presence of unobserved heterogeneity in technology types. One important caveat of our empirical exercise is that the class of production functions that we investigate may be restrictive. Investigating production function heterogeneity with more flexible function forms is an important future research topic.