#  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_ADDRESS>Address for correspondence: Yu (Jasmine) Hao, Faculty of Business
and Economics, The University of Hong Kong. We are very grateful for the
comments from Chun Pang Chow, Vadim Marmer, and Kevin Song. We are also
grateful for the IAAE grant at the 2019 IAAE Conference. This research is
supported by the Natural Science and Engineering Research Council of Canada.
Hiroyuki Kasahara
Vancouver School of Economics
The University of British Columbia
<EMAIL_ADDRESS>
This paper introduces a likelihood ratio-based test for examining the null
hypothesis of an $M_{0}$-component model versus an alternative
$(M_{0}+1)$-component model in the context of normal mixture panel regression.
We demonstrate that the first-order derivative of the density function for the
variance parameter in the normal mixture panel regression model is linearly
independent from its second-order derivative for the mean parameter, unlike
the cross-sectional normal mixture. However, similar to the cross-sectional
normal mixture, the likelihood ratio test statistic for the panel normal
mixture remains unbounded. To manage this unboundedness, we use a penalized
maximum likelihood estimator and derive the asymptotic distribution of the
penalized likelihood ratio test and expectation-maximization test statistics
using a fourth-order Taylor expansion of the log-likelihood function for
reparameterized parameters. A sequential hypothesis testing approach is
developed for consistently estimating the number of components. Simulation
experiments reveal that the proposed tests have good finite sample
performance. We apply these tests to estimate the number of production
technology types for a finite mixture Cobb–Douglas production function model.
Our findings suggest heterogeneity in output elasticities for intermediate
goods, indicating firm-specific variation in production functions beyond
Hicks-neutral productivity terms.

##  1 Introduction

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.111For example, Heckman and Singer, (1984) 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. Keane
and Wolpin, (1997) and Cameron and Heckman, (1998) analyze a dynamic model of
schooling and occupational choices with unobserved heterogeneous human
capital. Likewise, finite mixture models have been applied in health
economics. Deb and Trivedi, (1997) develop a finite mixture negative binomial
count model that accounts for the unobserved dispersion of medical care
utilization by the elderly. Kamakura and Russell, (1989) and Andrews and
Currim, (2003) 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 Titterington et al., (1985), Lindsay, (1995), and McLachlan and Peel,
(2004).

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 (see, e.g., Ghosh and Sen,, 1985;
Chernoff and Lander,, 1995; Lemdani and Pons,, 1997; Chen and Chen,, 2001,
2003; Chen et al.,, 2004; Garel,, 2001, 2005; Chen et al.,, 2014), and the
asymptotic distribution of the LRT statistic for general finite mixture models
has been derived as a function of the Gaussian process (Dacunha-Castelle and
Gassiat,, 1999; Liu and Shao,, 2003; Zhu and Zhang,, 2004; Azaïs et al.,,
2009). 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. Kasahara and
Shimotsu, (2015) analyze the asymptotic distribution of the LRT statistics of
a cross-sectional univariate finite mixture normal regression model, and
Kasahara and Shimotsu, (2019) 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.222Kasahara and Shimotsu,
(2014) 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 Chen and Li, (2009) and Kasahara and
Shimotsu, (2015), 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. Kasahara and
Shimotsu, (2015) and Kasahara and Shimotsu, (2019) 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 (Hartigan,, 1985), 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 (Hao,,
2017) 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 Olley and Pakes, (1996), Levinsohn and Petrin, (2003),
and Ackerberg et al., (2015). Our empirical finding suggests that it is
essential to incorporate unobserved heterogeneity into the production function
coefficients across firms in applications (Li and Sasaki,, 2017; Doraszelski
and Jaumandreu,, 2018; Balat et al.,, 2019; Kasahara et al.,, 2022).

The EM test approach was introduced by Li et al., (2009) and Chen and Li,
(2009) to test homogeneity in finite mixture models. Li and Chen, (2010)
develop an EM test for the null hypothesis of $M_{0}$ components applicable to
general $M_{0}\geq 2$, and Kasahara and Shimotsu, (2015) 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 (Niu et
al.,, 2011) and subgroup analyses (Shen and He,, 2015). More recently, Liu et
al., (2018) extend the EM test to mixtures of the general location-scale
family distribution, and Kasahara and Shimotsu, (2019) 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 (Kasahara and Shimotsu,, 2009; Lin and
Ng,, 2012; Bonhomme and Manresa,, 2015; Ando and Bai,, 2016; Su et al.,, 2016;
Lu and Su,, 2017). 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 2, we define the
finite normal mixture panel regression model. In Section 3, 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 4 generalizes the result to testing $M_{0}$ components
against $M_{0}+1$ components. Section 5 introduces the EM test for testing
$M_{0}$ components against $M_{0}+1$ components. Section 6 derives the
asymptotic distribution of the PLRT and EM tests under local alternatives.
Section 7 develops a consistent estimator for the number of components based
on sequential hypothesis testing. Section 8 presents the simulated results of
the tests. Section 9 provides an empirical application. In the following, $:=$
denotes “equals by definition,” and boldface letters denote vectors or
matrices.

##  2 Heteroskedastic finite mixture panel normal regression model

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\geq 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

$\begin{split}f_{M}(\boldsymbol{w};\boldsymbol{\vartheta}_{M})&=\sum_{j=1}^{M}\alpha_{j}f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\theta}_{j}),\end{split}$
(1)

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

$f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\theta}_{j})=\prod_{t=1}^{T}\frac{1}{\sigma_{j}}\phi\left(\frac{y_{t}-\mu_{j}-\boldsymbol{x}^{\top}_{t}\boldsymbol{\beta}_{j}-\boldsymbol{z}^{\top}_{t}\boldsymbol{\gamma}}{\sigma_{j}}\right)$
(2)

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 (1). 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
(1) 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.333More
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.$

##  3 Likelihood ratio test for $H_{0}:M=1$ against $H_{A}:M=2$

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 (2). Now consider a two-component mixture density function

$f_{2}(\boldsymbol{w};\boldsymbol{\vartheta}_{2})=\alpha
f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\theta}_{1})+(1-\alpha)f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\theta}_{2}),$

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 (cf. Chen and Li,, 2009): (i) the Fisher
information for testing $H_{02}$ is not finite, (ii) the log-likelihood
function is unbounded (Hartigan,, 1985), 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 5. 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 6.

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 (1) 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})$.

######  Proposition 1.

Suppose that the true model is described by the one-component model with
$(\boldsymbol{\gamma},\boldsymbol{\theta})=(\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*})$.
Then, for any positive constant $M>0$,
$\Pr\Big{(}LR_{n}^{*}(\tilde{\boldsymbol{\vartheta}}_{2,n})\leq M\Big{)}\to 0$
as $n\to\infty$.

To deal with unboundedness, we consider a penalized maximum likelihood
estimator (PMLE) as in Chen and Tan, (2009) using the following penalty
function:

$\tilde{p}_{n}({\boldsymbol{\vartheta}}_{M}):=\sum_{j=1}^{M}p_{n}(\sigma_{j}^{2})\quad\text{with}\quad
p_{n}(\sigma_{j}^{2}):=-a_{n}\\{{{\sigma}_{0}^{2}}/{\sigma_{j}^{2}}+\log({\sigma_{j}^{2}}/{{\sigma}_{0}^{2}})-1\\}.$
(3)

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 6.

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.

######  Assumption 1.

(a) $\boldsymbol{X}$ and $\boldsymbol{Z}$ have finite second moments, and
$\text{Pr}(\boldsymbol{X}^{\top}\boldsymbol{\beta}+\boldsymbol{Z}^{\top}\boldsymbol{\gamma}\neq\boldsymbol{X}^{\top}\boldsymbol{\beta}^{*}+\boldsymbol{Z}^{\top}\boldsymbol{\gamma}^{*})>0$
for
$(\boldsymbol{\beta}^{\top},\boldsymbol{\gamma}^{\top})^{\top}\neq((\boldsymbol{\beta}^{*})^{\top},(\boldsymbol{\gamma}^{*})^{\top})^{\top}$.
(b) $a_{n}>0$ and $a_{n}=o(n^{1/4})$ in the penalty function (3).

######  Proposition 2.

Suppose that Assumption 1 holds. Then, under the null hypothesis
$H_{0}:M_{0}=1$,
$\inf_{\boldsymbol{\vartheta}_{2}\in\Theta_{2}^{*}}|\hat{\boldsymbol{\vartheta}}_{2}-\boldsymbol{\vartheta}_{2}|\to_{p}0$.

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 2 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 3(c)).
Consequently, the panel mixture model (1) with the component density function
(2) 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
(cf. Chen,, 1995). See Proposition 4(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 (cf. Kasahara and Shimotsu,, 2015).

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

$\displaystyle\nabla_{\boldsymbol{\theta}_{1}}\log
f_{2}(\boldsymbol{w};\boldsymbol{\vartheta}^{*}_{2})=\frac{\alpha}{1-\alpha}\nabla_{\boldsymbol{\theta}_{2}}\log
f_{2}(\boldsymbol{w};\boldsymbol{\vartheta}^{*}_{2}).$ (4)

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 Kasahara and Shimotsu, (2012) and
consider the following one-to-one reparameterization of
$\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$ given $\alpha$:

$\begin{pmatrix}\boldsymbol{\lambda}\\\
\boldsymbol{\nu}\end{pmatrix}:=\begin{pmatrix}\boldsymbol{\theta}_{1}-\boldsymbol{\theta}_{2}\\\
\alpha\boldsymbol{\theta}_{1}+(1-\alpha)\boldsymbol{\theta}_{2}\end{pmatrix}\text{
so that }\begin{pmatrix}\boldsymbol{\theta}_{1}\\\
\boldsymbol{\theta}_{2}\end{pmatrix}=\begin{pmatrix}\boldsymbol{\nu}+(1-\alpha)\boldsymbol{\lambda}\\\
\boldsymbol{\nu}-\alpha\boldsymbol{\lambda}\end{pmatrix},$ (5)

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 (4). 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

$\displaystyle g(\boldsymbol{w};\boldsymbol{\psi},\alpha)$
$\displaystyle=\alpha
f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\nu}+(1-\alpha)\boldsymbol{\lambda})+(1-\alpha)f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\nu}-\alpha\boldsymbol{\lambda})\quad\text{and}$
(6) $\displaystyle l(\boldsymbol{w};\boldsymbol{\psi},\alpha)$
$\displaystyle=\log g(\boldsymbol{w};\boldsymbol{\psi},\alpha).$

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:

$\begin{split}\nabla_{\boldsymbol{\eta}^{\top}}l(\boldsymbol{w};\boldsymbol{\psi}^{*},\alpha)=\frac{\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\theta}^{\top})^{\top}}f(\boldsymbol{w};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*})}{f(\boldsymbol{w};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*})}\quad\text{and}\quad\nabla_{\boldsymbol{\lambda}}l(\boldsymbol{w};\boldsymbol{\psi}^{*},\alpha)=\boldsymbol{0}.\end{split}$
(7)

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}$:

$\nabla_{\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top}}l(\boldsymbol{w};\boldsymbol{\psi}^{*},\alpha)=\alpha(1-\alpha)\frac{\nabla_{\boldsymbol{\theta}\boldsymbol{\theta}^{\top}}f(\boldsymbol{w};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*})}{f(\boldsymbol{w};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*})}.$
(8)

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 Kasahara and Shimotsu,
(2015), 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

$\displaystyle\boldsymbol{s}(\boldsymbol{W})=\begin{pmatrix}\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{W})\\\
\boldsymbol{s}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}(\boldsymbol{W})\end{pmatrix},\quad\text{where}\underset{(p+q+2)\times
1}{\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{W})}:=\frac{\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\theta}^{\top})^{\top}}f^{*}}{f^{*}}\quad\text{
and }\quad\underset{((q+2)(q+1)/2)\times
1}{\boldsymbol{s}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}(\boldsymbol{W})}:=\frac{\widetilde{\nabla}_{\boldsymbol{\theta}\boldsymbol{\theta}^{\top}}f^{*}}{f^{*}}.$
(9)

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 B.2.

Collect the relevant normalized reparameterized parameters and define
$\boldsymbol{t}(\boldsymbol{\psi},\alpha)$ as

$\boldsymbol{t}(\boldsymbol{\psi},\alpha)=\begin{pmatrix}\boldsymbol{t}_{\boldsymbol{\eta}}\\\
\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)\end{pmatrix}=\begin{pmatrix}\boldsymbol{\eta}-\boldsymbol{\eta}^{*}\\\
\alpha(1-\alpha)\boldsymbol{v}(\boldsymbol{\lambda})\end{pmatrix},$ (10)

where $v(\boldsymbol{\lambda})$ is a vector of unique elements of
$\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top}$ given by

$v(\boldsymbol{\lambda})=(\lambda_{1}\lambda_{1},...,\lambda_{q+2}\lambda_{q+2},\lambda_{1}\lambda_{2},...,\lambda_{q+1}\lambda_{q+2})^{\top},$
(11)

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

$\displaystyle 2\\{L_{n}(\boldsymbol{\psi},\alpha)-$ $\displaystyle
L_{n}(\boldsymbol{\psi}^{*},\alpha)\\}=2(\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha))^{\top}\boldsymbol{S}_{n}-(\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha))^{\top}\boldsymbol{\mathcal{I}}_{n}(\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha))+R_{n}(\boldsymbol{\psi},\alpha)$
(12)
$\displaystyle\quad=\boldsymbol{G}_{n}^{\top}\boldsymbol{\mathcal{I}}_{n}\boldsymbol{G}_{n}-\left[\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha)-\boldsymbol{G}_{n}\right]^{\top}\boldsymbol{\mathcal{I}}_{n}\left[\sqrt{n}\boldsymbol{t}(\boldsymbol{\psi},\alpha)-\boldsymbol{G}_{n}\right]+R_{n}(\boldsymbol{\psi},\alpha),$
(13)

where $\boldsymbol{\mathcal{I}}_{n}$ is the negative of the sample Hessian
defined in the proof of Proposition 3 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}]$.

######  Assumption 2.

(a) $\boldsymbol{X}$ and $\boldsymbol{Z}$ have finite $8$-th moments. (b)
$\mathbb{E}[\boldsymbol{U}\boldsymbol{U^{\top}}]$ is non-singular, where
$\boldsymbol{U}=[1,\boldsymbol{X^{\top}},\boldsymbol{Z^{\top}}]^{\top}$.

######  Proposition 3.

Suppose that Assumptions 1 and 2 hold. Then, under $H_{0}:M=1$, for
$\alpha\in(0,1)$, (a) for any $\delta>0$,
$\lim\sup_{n\rightarrow\infty}\Pr(\sup_{\boldsymbol{\psi}\in\Theta_{\boldsymbol{\psi}}:||\boldsymbol{\psi}-\boldsymbol{\psi}^{*}||\leq\kappa}|R_{n}(\boldsymbol{\psi},\alpha)|>\delta(1+||n\boldsymbol{t}(\boldsymbol{\psi},\alpha)||^{2}))\rightarrow
0$ as $\kappa\rightarrow 0$, (b)
$\boldsymbol{S}_{n}\overset{d}{\to}\boldsymbol{S}\sim
N(0,\boldsymbol{\mathcal{I}})$, and (c)
$\boldsymbol{\mathcal{I}}_{n}\overset{p}{\to}\boldsymbol{\mathcal{I}}$, where
$\boldsymbol{\mathcal{I}}$ is finite and non-singular.

The non-singularity of $\boldsymbol{\mathcal{I}}$ in Proposition 3(c)
highlights the difference between the cross-sectional normal mixture and the
panel data normal mixture models. In particular, as shown in equation (75) in
Appendix B.2, 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 Andrews, (1999) to approximate
the distribution of the supremum of the right-hand side of (13) 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})^{\prime}\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

$\hat{\boldsymbol{\psi}}=\arg\max_{\boldsymbol{\psi}\in\Theta_{\boldsymbol{\psi}}}L_{n}(\boldsymbol{\psi},\alpha)+\sum_{j=1}^{2}p_{n}(\sigma_{j}^{2}(\boldsymbol{\psi},\alpha))$
(14)

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,

$\displaystyle LR_{n}$
$\displaystyle:=\max_{\alpha\in[\epsilon,1-\epsilon]}2\\{L_{n}(\hat{\boldsymbol{\psi}},\alpha)-L_{0,n}(\hat{\boldsymbol{\gamma}}_{0},\hat{\boldsymbol{\theta}}_{0})\\}\
\text{and}\
PLR_{n}:=LR_{n}+\sum_{j=1}^{2}p_{n}(\sigma_{j}^{2}(\hat{\boldsymbol{\psi}},\alpha)).$
(15)

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 5, we also develop the EM test, which does not impose a
direct constraint on the value of $\alpha$.

With $\boldsymbol{s}(\boldsymbol{W})$ in (9), partition
$\boldsymbol{\mathcal{I}}=\mathbb{E}[\boldsymbol{s}(\boldsymbol{W})\boldsymbol{s}(\boldsymbol{W})^{\top}]$
and define

$\begin{split}\boldsymbol{\mathcal{I}}=\begin{pmatrix}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}}&\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}\boldsymbol{\lambda}}\\\
\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}&\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}\end{pmatrix},\quad\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}}=\mathbb{E}[\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{W})\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{W})^{\top}],\quad\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}=\mathbb{E}[\boldsymbol{s}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}(\boldsymbol{W})\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{W})^{\top}],\quad\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}\boldsymbol{\lambda}}=\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}^{\top},\\\
\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}=\mathbb{E}[\boldsymbol{s}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}(\boldsymbol{W})\boldsymbol{s}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}(\boldsymbol{W})^{\top}],\quad\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}=\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\lambda}}-\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}}^{-1}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}\boldsymbol{\lambda}},\quad\text{and}\quad\boldsymbol{G}_{\boldsymbol{\lambda},\boldsymbol{\eta}}:=(\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}})^{-1}\boldsymbol{S}_{\boldsymbol{\lambda},\boldsymbol{\eta}},\end{split}$

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

$\Lambda_{\boldsymbol{\lambda}}=\Big{\\{}\sqrt{n}\alpha(1-\alpha)v(\boldsymbol{\lambda}):\boldsymbol{\lambda}\in\Theta_{\boldsymbol{\lambda}}\Big{\\}}.$

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

$r_{\lambda}(\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}})=\inf_{\boldsymbol{t}_{\boldsymbol{\lambda}}\in\Lambda_{\boldsymbol{\lambda}}}r_{\lambda}(\boldsymbol{t}_{\boldsymbol{\lambda}}),\quad
r_{\lambda}(\boldsymbol{t}_{\boldsymbol{\lambda}}):=(\boldsymbol{t}_{\boldsymbol{\lambda}}-\boldsymbol{G}_{\boldsymbol{\lambda},\boldsymbol{\eta}})^{\top}\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}(\boldsymbol{t}_{\boldsymbol{\lambda}}-\boldsymbol{G}_{\boldsymbol{\lambda},\boldsymbol{\eta}}),$
(16)

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 4.

Suppose that Assumptions 1 and 2 hold. Under the null hypothesis
$H_{0}:M_{0}=1$, (a)
$\boldsymbol{t}(\hat{\boldsymbol{\psi}},\alpha)=O_{p}(n^{-1/2})$ for any
$\alpha\in(0,1)$, (b)
$LR_{n}\overset{d}{\to}(\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}})^{\top}\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}}$
and
$PLR_{n}\overset{d}{\to}(\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}})^{\top}\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}}+\text{plim}_{n\rightarrow\infty}\sum_{j=1}^{2}p_{n}(\sigma_{j}^{2}(\hat{\boldsymbol{\psi}},\alpha))$.

Proposition 4(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 (4), 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 (Chen,,
1995). 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}$.

##  4 Likelihood ratio test for $H_{0}:M=M_{0}$ against $H_{A}:M=M_{0}+1$

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}\geq 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
(17):

$f_{M_{0}}(\boldsymbol{w};\boldsymbol{\vartheta}_{M_{0}}^{*})=\sum_{j=1}^{M_{0}}\alpha^{*}_{j}f(\boldsymbol{w};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{j}),$
(17)

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

$f_{M_{0}+1}(\boldsymbol{w};\boldsymbol{\vartheta}_{M_{0}+1})=\sum_{j=1}^{M_{0}+1}\alpha_{j}f(\boldsymbol{w};\boldsymbol{\gamma},\boldsymbol{\theta}_{j}),$
(18)

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 (17). We assume
$\mu_{1}^{*}<\mu_{2}^{*},\ldots,<\mu_{M_{0}}^{*}$ in the true parameters for
identification.

The $(M_{0}+1)$-component model (18) gives rise to the true density (17) 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

$\begin{split}\Upsilon_{2h}^{*}=\\{\boldsymbol{\vartheta}_{M_{0}+1}\in\Theta_{\boldsymbol{\vartheta}_{M_{0}+1}}:\alpha_{h}=0;(\alpha_{j},\mu_{j},\sigma_{j})=(\alpha_{j}^{*},\mu_{j}^{*},\sigma_{j}^{*})\text{
for }j<h;\\\
(\alpha_{j},\mu_{j},\sigma_{j})=(\alpha_{j-1}^{*},\mu_{j-1}^{*},\sigma_{j-1}^{*})\text{
for }j>h\\}.\end{split}$

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}$.

######  Proposition 5.

$\sup_{\vartheta_{M_{0}+1}\in\Upsilon_{2h}^{*}}\mathbb{E}[\\{\nabla_{\alpha_{h}}\log
f_{M_{0}+1}(\boldsymbol{W}_{i},\vartheta_{M_{0}+1})\\}^{2}]<\infty$ if and
only if
$\max\\{\sigma^{2}:\sigma\in\Theta_{\sigma}\\}<2\max\\{\sigma_{1}^{2*},\ldots,\sigma_{M_{0}}^{2*}\\}.$

Because the restriction on the values of $\sigma_{j}^{2}$ in Proposition 5 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 3.

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

$\begin{split}\Upsilon^{*}_{1h}:=\Big{\\{}\boldsymbol{\vartheta}_{M_{0}+1}\in\Theta_{\vartheta_{M_{0}+1}}:\alpha_{h}+\alpha_{h+1}=\alpha_{h}^{*}\text{
and
}\boldsymbol{\theta}_{h}=\boldsymbol{\theta}_{h+1}=\boldsymbol{\theta}_{h}^{*};\boldsymbol{\gamma}=\boldsymbol{\gamma}^{*};\alpha_{j}=\alpha_{j}^{*}\qquad\\\
\text{ and }\boldsymbol{\theta}_{j}=\boldsymbol{\theta}_{j}^{*}\text{ for
}1\leq j<h;\alpha_{j}=\alpha_{j-1}^{*}\text{ and
}\boldsymbol{\theta}_{j}=\boldsymbol{\theta}_{j-1}^{*}\text{ for }h+1\leq
j\leq M_{0}+1\Big{\\}}\end{split}$ (19)

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
(17). 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

$\boldsymbol{\Psi}_{h}^{*}=\left\\{\begin{split}\alpha_{1},\ldots,\alpha_{M_{0}+1}\in[\epsilon,1-\epsilon];\sum_{j=1}^{M_{0}+1}\alpha_{j}=1;\gamma\in\Theta_{\gamma};\theta\in\Theta_{\theta}:\boldsymbol{\theta}_{j}\in
D_{j}^{*}\text{ for }j=1,\ldots,h-1;\\\
\boldsymbol{\theta}_{h},\boldsymbol{\theta}_{h+1}\in
D^{*}_{h};\boldsymbol{\theta}_{j}\in D^{*}_{j-1}\text{ for
}j=h+2,\ldots,M_{0}+1.\end{split}\right\\}$ (20)

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 (20)
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

$p_{n}(\sigma_{j}^{2};\hat{\sigma}_{0,j}^{2}):=-a_{n}\\{{\hat{\sigma}_{0,j}^{2}}/{\sigma_{j}^{2}}+\log({\sigma_{j}^{2}}/{\hat{\sigma}_{0,j}^{2}})-1\\},$
(21)

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 4(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 (19) 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.

######  Proposition 6.

Suppose that Assumption 1 holds. Then, under the null hypothesis
$H_{0}:M=M_{0}$,
$\inf_{\boldsymbol{\vartheta}_{M_{0}+1}\in\boldsymbol{\Psi}_{h}^{*}}|\hat{\boldsymbol{\vartheta}}_{M_{0}+1}^{h}-\boldsymbol{\vartheta}_{M_{0}+1}|\overset{p}{\to}0$
for $h=1,2,...,M_{0}$.

Consider the local PLRTS for testing
$H_{0,1h}:\boldsymbol{\vartheta}_{h}=\boldsymbol{\vartheta}_{h+1}$ defined by

$PLR^{M_{0},h}_{n}:=2\\{L_{M_{0}+1,n}(\hat{\boldsymbol{\vartheta}}^{h}_{M_{0}+1})+\tilde{p}_{n}(\hat{\boldsymbol{\vartheta}}^{h}_{M_{0}+1})-L_{M_{0},n}(\hat{\boldsymbol{\vartheta}}_{M_{0}})\\}\quad\text{for
$h=1,2,...,M_{0}$}.$

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

$PLR_{n}(M_{0}):=\max\\{PLR^{M_{0},1}_{n},PLR^{M_{0},2}_{n},...,PLR^{M_{0},M_{0}}_{n}\\}.$
(22)

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

$\tilde{\boldsymbol{s}}(\boldsymbol{W})=\begin{pmatrix}\tilde{\boldsymbol{s}}_{\boldsymbol{\eta}}(\boldsymbol{W})\\\
\tilde{\boldsymbol{s}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})\end{pmatrix},\
\text{ where }\underset{(M_{0}+p+q+1)\times
1}{\tilde{\boldsymbol{s}}_{\boldsymbol{\eta}}(\boldsymbol{W})}=\begin{pmatrix}\boldsymbol{s}_{\boldsymbol{\alpha}}(\boldsymbol{W})\\\
\boldsymbol{s}_{(\boldsymbol{\gamma},\boldsymbol{\nu})}(\boldsymbol{W})\end{pmatrix}\quad\text{
and
}\tilde{\boldsymbol{s}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})=\begin{pmatrix}\boldsymbol{s}^{1}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})\\\
\vdots\\\
\boldsymbol{s}^{M_{0}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})\end{pmatrix},$
(23)

where

$\begin{split}\boldsymbol{s}_{\boldsymbol{\alpha}}(\boldsymbol{W})&=\begin{pmatrix}f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{1})-f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{M_{0}})\\\
\vdots\\\
f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{M_{0}-1})-f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{M_{0}})\end{pmatrix}\Bigg{/}f_{M_{0}}(\boldsymbol{W};\boldsymbol{\vartheta}_{M_{0}}^{*}),\\\
\boldsymbol{s}_{(\boldsymbol{\gamma},\boldsymbol{\nu})}(\boldsymbol{W})&=\sum_{j=1}^{M_{0}}\alpha_{j}^{*}\nabla_{(\boldsymbol{\gamma},\boldsymbol{\nu})}f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{j})/f_{M_{0}}(\boldsymbol{W};\boldsymbol{\vartheta}_{M_{0}}^{*}),\\\
\boldsymbol{s}^{h}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})&=\widetilde{\nabla}_{\boldsymbol{\theta}_{h}\boldsymbol{\theta}_{h}^{\top}}f(\boldsymbol{W};\boldsymbol{\gamma}^{*},\boldsymbol{\theta}^{*}_{h})/f_{M_{0}}(\boldsymbol{W};\boldsymbol{\vartheta}_{M_{0}}^{*})\quad\text{for
$h=1,2,...,M_{0}$},\end{split}$ (24)

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

$\begin{split}\tilde{\boldsymbol{\mathcal{I}}}:=\mathbb{E}[\boldsymbol{\tilde{s}}(\boldsymbol{W})\boldsymbol{\tilde{s}}(\boldsymbol{W})^{\top}],\quad\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}}:=\mathbb{E}[\boldsymbol{\tilde{s}}_{\boldsymbol{\eta}}(\boldsymbol{W})\boldsymbol{\tilde{s}}_{\boldsymbol{\eta}}(\boldsymbol{W})^{\top}],\quad\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}:=\mathbb{E}[\boldsymbol{\tilde{s}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})\boldsymbol{\tilde{s}}_{\boldsymbol{\eta}}(\boldsymbol{W})^{\top}],\\\
\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}\boldsymbol{\lambda}}:=\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}^{\top},\quad\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda\lambda}}:=\mathbb{E}[\boldsymbol{\tilde{s}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})\boldsymbol{\tilde{s}}_{\boldsymbol{\lambda\lambda}}(\boldsymbol{W})^{\top}],\quad\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\boldsymbol{\lambda},\boldsymbol{\eta}}}:=\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda\lambda}}-\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda\eta}}\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}}^{-1}\tilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta\lambda}}.\end{split}$
(25)

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 (23), $\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

$\begin{split}r^{h}_{\boldsymbol{\lambda}}(\hat{\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}})=\inf_{{\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}}\in\Lambda_{\boldsymbol{\lambda}}}r^{h}({\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}});\quad
r^{h}_{\boldsymbol{\lambda}}({\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}}):=({\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}}-\boldsymbol{{G}}^{h}_{\boldsymbol{\lambda},\boldsymbol{\eta}})^{\top}{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}^{h}({\boldsymbol{t}}^{h}_{\boldsymbol{\lambda}}-\boldsymbol{{G}}^{h}_{\boldsymbol{\lambda},\boldsymbol{\eta}})\quad\text{for
$h=1,2,...,M_{0}$}.\end{split}$ (26)

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 3 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 (22), 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}$.

######  Assumption 3.

(a) $\alpha_{j}^{*}\in(\epsilon,1-\epsilon)$ for $j=1,\ldots,M_{0}$. (b)
$\boldsymbol{\tilde{\mathcal{I}}}$ is non-singular. (c) $a_{n}$ in (21)
satisfies $a_{n}=O(1)$.

######  Proposition 7.

Suppose that Assumptions 1–3 are satisfied. Then, under the null hypothesis
$H_{0}:M=M_{0}$,
$PLR_{n}({M_{0}})\overset{d}{\to}\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_{0}}_{\boldsymbol{\lambda}})^{\top}\boldsymbol{\mathcal{I}}^{M_{0}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}\hat{\boldsymbol{t}}^{M_{0}}_{\boldsymbol{\lambda}}\\}$.

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
(26) 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 B.2 and B.3 present an expression for the score functions using
Hermit polynomials.

##  5 EM test for $H_{0}:M=M_{0}$ against $H_{A}:M=M_{0}+1$

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

$\begin{split}w_{ij}^{(k)}&=\left\\{\begin{split}&\alpha_{j}^{(k)}f(\boldsymbol{W}_{i};\boldsymbol{\gamma}^{(k)},\boldsymbol{\theta}_{j}^{(k)})/f_{M_{0}+1}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M_{0}+1}^{h(k)}(\tau_{0})),j=1,\ldots,h-1,\\\
&\alpha_{j-1}^{(k)}f(\boldsymbol{W}_{i};\boldsymbol{\gamma}^{(k)},\boldsymbol{\theta}_{j}^{(k)})/f_{M_{0}+1}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M_{0}+1}^{h(k)}(\tau_{0})),j=h+2,\ldots,M_{0}+1,\end{split}\right.\\\
w_{ih}^{(k)}&=\tau^{h(k)}\alpha_{h}^{(k)}f(\boldsymbol{W}_{i};\boldsymbol{\gamma}^{(k)},\boldsymbol{\theta}_{h}^{(k)})/f_{M_{0}+1}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M_{0}+1}^{h(k)}(\tau_{0})),\\\
w_{i,h+1}^{(k)}&=(1-\tau^{h(k)})\alpha_{h}^{(k)}f(\boldsymbol{W}_{i};\boldsymbol{\gamma}^{(k)},\boldsymbol{\theta}_{h+1}^{(k)})/f_{M_{0}+1}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M_{0}+1}^{h(k)}(\tau_{0})),\end{split}$
(27)

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

$\displaystyle\alpha_{j}^{(k+1)}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{w_{ij}^{(k)}}\quad\text{for
}j=1,\ldots,M_{0}+1\quad\text{and}$ $\displaystyle\tau^{h(k+1)}$
$\displaystyle=\arg\min_{\tau}\left\\{\sum_{i=1}^{n}{w_{ih}^{(k)}}\log(\tau)+\sum_{i=1}^{n}{w_{i,h+1}^{(k)}}\log(1-\tau)+p(\tau)\right\\}.$

We also update $\boldsymbol{\theta}_{j}$ and $\boldsymbol{\gamma}$ as

$\displaystyle(\sigma_{j}^{(k+1)})^{2}$
$\displaystyle=\arg\min_{\sigma_{j}^{2}}\left\\{\sum_{i=1}^{n}w_{i}^{j(k)}\sum_{t=1}^{T}(y_{it}-\mu_{j}^{(k+1)}-\boldsymbol{z}_{it}^{\top}\boldsymbol{\gamma}^{(k+1)}-\boldsymbol{x}_{it}^{\top}\boldsymbol{\beta}_{j}^{(k+1)})^{2}+p_{n}(\sigma_{j}^{2})\right\\},$
$\displaystyle\boldsymbol{\gamma}^{(k+1)}$
$\displaystyle=\left(\sum_{i=1}^{n}\sum_{t=1}^{T}\boldsymbol{z}_{it}\boldsymbol{z}_{it}^{\top}\right)^{-1}\left(\sum_{i=1}^{n}\sum_{t=1}^{T}\boldsymbol{z}_{it}\left(y_{it}-\sum_{j=1}^{M_{0}+1}w_{ij}^{(k)}\tilde{\boldsymbol{x}}_{it}^{\top}\begin{pmatrix}\mu_{j}^{(k)}\\\
\boldsymbol{\beta}_{j}^{(k)}\end{pmatrix}\right)\right),\quad\text{and}$
$\displaystyle\begin{pmatrix}\mu_{j}^{(k+1)}\\\
\boldsymbol{\beta}_{j}^{(k+1)}\end{pmatrix}$
$\displaystyle=\left(\sum_{i=1}^{n}w_{ij}^{(k)}\sum_{t=1}^{T}\tilde{\boldsymbol{x}}_{it}\tilde{\boldsymbol{x}}_{it}^{\top}\right)^{-1}\left(\sum_{i=1}^{n}w_{ij}^{(k)}\sum_{t=1}^{T}\tilde{\boldsymbol{x}}_{it}(y_{it}-\boldsymbol{z}_{it}^{\top}\boldsymbol{\gamma}^{(k+1)})\right),$

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

$M_{n}^{h(k)}(\tau_{0}):=2\left\\{PL_{n}(\boldsymbol{\vartheta}_{M_{0}+1}^{h(k)}(\tau_{0}),\tau^{h(k)}(\tau_{0}))-L_{M_{0},n}(\hat{\boldsymbol{\vartheta}}_{M_{0}})\right\\}.$
(28)

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

$EM_{n}^{h}:=\max\\{M^{h(K)}_{n}(\tau_{0}):\tau_{0}\in\mathcal{T}\\}.$ (29)

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:

$EM_{n}(M_{0}):=\max\\{EM_{n}^{1(K)},\ldots,EM_{n}^{M_{0}(K)}\\}.$ (30)

######  Proposition 8.

Suppose that Assumptions 1–3 hold and $\\{0.5\\}\in\mathcal{T}$. Then, under
the null hypothesis $H_{0}:M=M_{0}$, for any finite $K$,
$EM_{n}(M_{0})\overset{d}{\to}\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_{0}}_{\boldsymbol{\lambda}})^{\top}\boldsymbol{\mathcal{I}}^{M_{0}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}\hat{\boldsymbol{t}}^{M_{0}}_{\boldsymbol{\lambda}}\\}$.

Therefore, the asymptotic null distribution of the EM test statistic
$EM_{n}(M_{0})$ is the same as that of the PLRTS.

##  6 Asymptotic distribution under local alternatives

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
(10),

$\boldsymbol{h_{\nu}}=\sqrt{n}(\boldsymbol{\nu}_{n}-\boldsymbol{\nu}^{*}),\quad\boldsymbol{h_{\lambda}}=\sqrt{n}\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda}_{n},\alpha_{n}),\quad\text{and}\quad\alpha_{n}=\alpha^{*}+o(1).$
(31)

Equivalently, the non-reparameterized contiguous local alternatives are given
by

$\boldsymbol{\theta}_{1,n}=\boldsymbol{\nu}_{n}+(1-\alpha_{n})\boldsymbol{\lambda}_{n}\quad\text{and}\quad\boldsymbol{\theta}_{2,n}=\boldsymbol{\nu}_{n}-\alpha_{n}\boldsymbol{\lambda}_{n}$
(32)

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

$\displaystyle\lambda_{j,n}=n^{-1/4}(\alpha_{n}(1-\alpha_{n}))^{-1/2}h_{\lambda,j}\quad\text{for
$j=1,...,q+2$},$

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 4.

The following proposition provides the asymptotic distribution of the PLRT and
EM test statistics under contiguous local alternatives.

######  Proposition 9.

Suppose that the assumptions in Proposition 8 hold for $M_{0}=1$. Consider a
sequence of contiguous local alternatives
$\boldsymbol{\vartheta}_{2,n}=(\alpha_{n},\boldsymbol{\theta}_{1,n}^{\top},\boldsymbol{\theta}_{2,n}^{\top})^{\top}$
given in (32), where $\alpha_{n}$ and $\boldsymbol{\lambda}_{n}$ satisfy (31).
Then, under $H_{1,n}:\boldsymbol{\vartheta}=\boldsymbol{\vartheta}_{2,n}$, we
have $PLR_{n}(1)$,
$EM_{n}(1)\overset{d}{\rightarrow}(\tilde{\boldsymbol{t}}_{\boldsymbol{\lambda}})^{\top}\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}}\tilde{\boldsymbol{t}}_{\boldsymbol{\lambda}}$,
where $\tilde{\boldsymbol{t}}_{\boldsymbol{\lambda}}$ has the same
distribution as $\hat{\boldsymbol{t}}_{\boldsymbol{\lambda}}$ in Proposition 4
but $\boldsymbol{G}_{\boldsymbol{\lambda,\boldsymbol{\eta}}}$ is replaced with
$(\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda},\boldsymbol{\eta}})^{-1}\boldsymbol{S}_{\boldsymbol{\lambda},\boldsymbol{\eta}}+\boldsymbol{h}_{\boldsymbol{\lambda}}$.

Importantly, a set of contiguous local alternatives considered in (32)
excludes a sequence such that $\alpha_{n}\rightarrow 0$ or $1$.

##  7 Sequential hypothesis testing

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 7 and 8. 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

$\displaystyle\hat{M}_{\text{PLR}}$
$\displaystyle=\min_{M\in\\{0,\ldots,\bar{M}\\}}\\{M:PLR_{n}(r)\geq\hat{c}^{r}_{1-q_{n}},r=0,\ldots,M-1,PLR_{n}({M})\
<\hat{c}^{M}_{1-q_{n}}\\},$ $\displaystyle\hat{M}_{\text{EM}}$
$\displaystyle=\min_{M\in\\{0,\ldots,\bar{M}\\}}\\{M:EM_{n}(r)\geq\hat{c}^{r}_{1-q_{n}},r=0,\ldots,M-1,EM_{n}({M})\
<\hat{c}^{M}_{1-q_{n}}\\}.$ (33)

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 (1) for
$M=1,...,\bar{M}$.

######  Assumption 4.

For $M=1,...,M_{0}-1$, (a) $Q^{M}(\boldsymbol{\vartheta}_{M})$ has a unique
maximum at $\boldsymbol{\vartheta}_{M}^{*}$ in
$\Theta_{\boldsymbol{\vartheta}_{M}}$; (b)
$\Theta_{\boldsymbol{\vartheta}_{M}}$ is compact; (c)
$\boldsymbol{\vartheta}_{M}^{*}$ is interior to
$\Theta_{\boldsymbol{\vartheta}_{M}}$; (d)
$B^{M}(\boldsymbol{\vartheta}_{M}^{*}):=\mathbb{E}\left\\{\nabla_{\boldsymbol{\vartheta}_{M}}\ln
f_{M}(\boldsymbol{w}_{i};\boldsymbol{\vartheta}_{M})\nabla_{\boldsymbol{\vartheta}_{M}^{\top}}\ln
f_{M}(\boldsymbol{w}_{i};\boldsymbol{\vartheta}_{M})\right\\}$ is non-
singular; and (e)
$A^{M}(\boldsymbol{\vartheta}_{M}^{*}):=\mathbb{E}\left\\{\nabla_{\boldsymbol{\vartheta}_{M}\boldsymbol{\vartheta}_{M}^{\top}}\ln
f_{M}(\boldsymbol{w}_{i};\boldsymbol{\vartheta}_{M})\right\\}$ has a constant
rank in some open neighborhood of $\boldsymbol{\vartheta}_{M}^{*}$; (f)
$Q^{M+1}(\boldsymbol{\vartheta}_{M+1}^{*})-Q^{M}(\boldsymbol{\vartheta}_{M}^{*})>0$.

######  Proposition 10.

Suppose that $M_{0}<\bar{M}$ and Assumptions 1–4 hold. If we choose $q_{n}$
such that $-n^{-1}\ln q_{n}=o(1)$ and $q_{n}=o(1)$, then
$\hat{M}_{\text{PLR}}-M_{0}=o_{p}(1)$ and
$\hat{M}_{\text{EM}}-M_{0}=o_{p}(1)$.

Assumptions 4(a)–(e) ensure the consistency and asymptotic normality of
$\hat{\boldsymbol{\vartheta}}_{M}$, where (c)–(e) correspond to Assumption A6
of White, (1982). Per Assumption 4(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}$.

##  8 Simulation

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$.

###  8.1 Choice of penalty function

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 11 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$:

$a_{n}=\begin{cases}\left({1+\exp\left\\{\frac{\hat{\rho}_{1}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}+\frac{\hat{\rho}_{2}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}\frac{1}{T}+\frac{\hat{\rho}_{3}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}\frac{1}{n}\right\\}}\right)^{-1},&M_{0}=1\\\
\left({1+\exp\left\\{\frac{\hat{\rho}_{1}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}+\frac{\hat{\rho}_{2}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}\frac{1}{T}+\frac{\hat{\rho}_{3}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}\frac{1}{n}+\frac{\hat{\rho}_{5}^{M_{0}}}{\hat{\rho}_{4}^{M_{0}}}\log\left(\frac{\omega(\boldsymbol{\vartheta}_{M_{0}};M_{0})}{1-\omega(\boldsymbol{\vartheta}_{M_{0}};M_{0})}\right)\right\\}}\right)^{-1},&M_{0}=2,3,4,\end{cases}$
(34)

where $\omega(\boldsymbol{\vartheta}_{M_{0}};M_{0})$ is the misclassification
probability as defined in Melnykov and Maitra, (2010) 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:

$\begin{split}&\log\left(\frac{\hat{s}}{1-\hat{s}}\right)-\log\left(\frac{{0.05}}{1-{0.05}}\right)\\\
&=\begin{cases}\rho_{1}^{M_{0}}+\rho_{2}^{M_{0}}\frac{1}{T}+\rho_{3}^{M_{0}}\frac{1}{n}+\rho_{4}^{M_{0}}\log\left(\frac{{a_{n}}}{1-{a_{n}}}\right),&M_{0}=1\\\
\rho_{1}^{M_{0}}+\rho_{2}^{M_{0}}\frac{1}{T}+\rho_{3}^{M_{0}}\frac{1}{n}+\rho_{4}^{M_{0}}\log\left(\frac{{a_{n}}}{1-{a_{n}}}\right)+\rho_{5}^{M_{0}}\log\left(\frac{\omega(\boldsymbol{\vartheta}_{M_{0}};M_{0})}{1-\omega(\boldsymbol{\vartheta}_{M_{0}};M_{0})}\right),&M_{0}=2,3,4,\end{cases}\end{split}$

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 (34) denote the corresponding estimates. Table 12
in the Appendix reports the estimates. Note that the data-dependent formula
(34) 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 11 for $M_{0}=2$ and taking
the average across the predicted values of $\hat{a}_{n}$. For $M_{0}\geq 5$,
we use the parametric bootstrap method to obtain the critical values for our
empirical application, where we set $a_{n}=0.5$.

###  8.2 Simulation results

Table 1 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\geq 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 2 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 1, 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 3 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 4 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 (33).

Table 5 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.

Table 1: Sizes (in %) of the EM Test and the PLRT of $H_{0}:M_{0}=2$ Against
$H_{A}:M_{0}=3$ at the $5\%$ Level

| EM Test | EM Test | PLRT
---|---|---|---
| Asymptotic | Parametric Bootstrap | Asymptotic
$T$ | 3 | 5 | 8 | 3 | 5 | 8 | 3 | 5 | 8
$N$ | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400 | 200 | 400
$(A,C)$ | 5.3 | 4.8 | 4 | 3.8 | 4 | 2.95 | 4.6 | 6.2 | 6.2 | 4.6 | 4.8 | 5.2 | 7.7 | 6.6 | 6.75 | 6.7 | 6.5 | 5.5
$(A,D)$ | 5.9 | 4.9 | 5 | 5 | 4.45 | 4 | 5.4 | 4.8 | 5.2 | 5.6 | 5 | 5.8 | 5.4 | 5.15 | 5.1 | 6.1 | 5.55 | 6.2
$(B,C)$ | 3.8 | 2.5 | 3.45 | 3.05 | 3.6 | 3.25 | 3.6 | 5.6 | 4.2 | 5.2 | 4 | 5.4 | 6.25 | 5.45 | 6.45 | 6.05 | 5.05 | 6
$(B,D)$ | 4.8 | 4.6 | 3.5 | 3.15 | 3.55 | 3.95 | 3.6 | 3.6 | 5.8 | 4 | 6.2 | 4.6 | 2.35 | 4.4 | 3.9 | 4.85 | 4.95 | 5.2

  * 1

$A$ and $B$ refer to, respectively, $(\alpha_{1},\alpha_{2})=(0.5,0.5)$ and
$(0.2,0.8)$, while $C$ and $D$ refer to $(\mu_{1},\mu_{2})=(-1,1)$ and
$(-0.5,0.5)$, respectively.

  * 2

The variance is set to $(\sigma_{1},\sigma_{2})=(0.8,1.2)$. The asymptotic
simulations are based on 2,000 repetitions and the bootstrap simulation is
based on 1,000 repetitions.

Table 2: Powers (in %) of the EM Test and the PLRT of $H_{0}:M_{0}=2$ Against
$H_{A}:M_{0}=3$ at the $5\%$ Level

| A | B
---|---|---
N | 100 | 500 | 100 | 500
T | 2 | 5 | 2 | 5 | 2 | 5 | 2 | 5
| EM test
$(C,G)$ | 20.9 | 81.6 | 57.6 | 100.0 | 20.5 | 82.7 | 62.6 | 100.0
$(C,H)$ | 49.2 | 99.9 | 99.9 | 100.0 | 38.4 | 98.7 | 98.8 | 100.0
$(C,I)$ | 12.1 | 20.4 | 18.0 | 62.6 | 10.6 | 20.4 | 16.8 | 65.8
$(D,G)$ | 77.9 | 100.0 | 100.0 | 100.0 | 86.5 | 100.0 | 100.0 | 100.0
$(D,H)$ | 57.4 | 100.0 | 100.0 | 100.0 | 42.8 | 100.0 | 100.0 | 100.0
$(D,I)$ | 16.0 | 59.5 | 31.8 | 99.9 | 13.8 | 70.8 | 40.3 | 100.0
$(E,G)$ | 93.0 | 100.0 | 100.0 | 100.0 | 94.0 | 100.0 | 100.0 | 100.0
$(E,H)$ | 83.8 | 100.0 | 100.0 | 100.0 | 70.7 | 100.0 | 100.0 | 100.0
$(E,I)$ | 25.7 | 97.0 | 80.2 | 100.0 | 30.7 | 96.8 | 83.1 | 100.0
$(F,G)$ | 99.9 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0
$(F,H)$ | 93.5 | 100.0 | 100.0 | 100.0 | 85.3 | 100.0 | 100.0 | 100.0
$(F,I)$ | 40.8 | 99.9 | 98.2 | 100.0 | 52.1 | 100.0 | 99.5 | 100.0
| PLRT
$(C,G)$ | 22.7 | 85.1 | 56.5 | 100.0 | 23.2 | 82.6 | 58.7 | 100.0
$(C,H)$ | 57.1 | 100.0 | 99.8 | 100.0 | 43.2 | 99.7 | 99.1 | 100.0
$(C,I)$ | 12.0 | 21.0 | 12.4 | 66.1 | 11.3 | 22.1 | 12.4 | 69.3
$(D,G)$ | 79.9 | 100.0 | 100.0 | 100.0 | 87.6 | 100.0 | 100.0 | 100.0
$(D,H)$ | 65.3 | 100.0 | 100.0 | 100.0 | 49.1 | 100.0 | 100.0 | 100.0
$(D,I)$ | 14.8 | 63.6 | 28.6 | 100.0 | 13.6 | 75.2 | 36.7 | 100.0
$(E,G)$ | 91.5 | 100.0 | 100.0 | 100.0 | 93.7 | 100.0 | 100.0 | 100.0
$(E,H)$ | 86.8 | 100.0 | 100.0 | 100.0 | 75.9 | 100.0 | 100.0 | 100.0
$(E,I)$ | 28.7 | 97.2 | 77.7 | 100.0 | 33.0 | 97.9 | 85.5 | 100.0
$(F,G)$ | 99.9 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 | 100.0
$(F,H)$ | 96.4 | 100.0 | 100.0 | 100.0 | 89.5 | 100.0 | 100.0 | 100.0
$(F,I)$ | 45.7 | 100.0 | 98.4 | 100.0 | 57.2 | 100.0 | 99.8 | 100.0

Notes: $A$ and $B$ refer to $(\alpha_{1},\alpha_{2},\alpha_{3})=(1/3,1/3,1/3)$
and $(1/4,1/2,1/4)$, respectively; $C,D,E,$ and $F$ refer to
$(\mu_{1},\mu_{2},\mu_{3})=(-0.5,0,1.5),(-1,0,1),(-1,0,2),$ and
$(-1.5,0,1.5)$, respectively; and $G,H,$ and $I$ refer to
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0.6,0.6,1.2),(0.6,1.2,0.6),$ and
$(1,1,1)$, respectively.

Table 3: Sizes (in %) of the EM Test of $H_{0}:M_{0}=3$ Against
$H_{A}:M_{0}=4$ at the $5\%$ Level

| (A,C) | (A,D) | (A,E) | (B,C) | (B,D) | (B,E)
---|---|---|---|---|---|---
100,2 | 5.95 | 5.15 | 5.05 | 5.05 | 5.85 | 4.40
500,2 | 5.60 | 5.55 | 5.25 | 5.10 | 5.65 | 4.05
100,5 | 4.30 | 6.00 | 4.20 | 5.15 | 5.10 | 5.70
500,5 | 4.20 | 4.55 | 3.95 | 4.50 | 4.15 | 4.15

Notes: $A$ and $B$ refer to $(\alpha_{1},\alpha_{2},\alpha_{3})=(1/3,1/3,1/3)$
and $(0.25,0.5,0.25)$, respectively, and $C,D,$ and $E$ refer to
$(\mu_{1},\mu_{2},\mu_{3})=(-4,0,4)$, $(-4,0,6),$ and $(-6,0,6)$,
respectively. The variance is set to
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0.75,1.5,0.75)$. The asymptotic
simulations are based on 2,000 repetitions and the bootstrap simulation is
based on 1,000 repetitions.

Table 4: Sizes of the EM Test of $H_{0}:M_{0}=2$ Against $H_{A}:M_{0}=3$ with
Conditioning Variables

| $(A,C,E)$ | $(A,C,F)$ | $(A,D,E)$ | $(A,D,F)$ | $(B,C,E)$ | $(B,C,F)$ | $(B,D,E)$ | $(B,D,F)$
---|---|---|---|---|---|---|---|---
$(N,T)$ |  |  |  |  |  |  |  |
$(200,2)$ | 8.4 | 8.2 | 7.4 | 8.8 | 8.6 | 8.2 | 7.4 | 3.6
$(500,2)$ | 4.6 | 3.2 | 3.2 | 2.2 | 4.8 | 4.8 | 3.6 | 3.6
$(200,5)$ | 4.0 | 1.8 | 3.0 | 2.6 | 2.2 | 2.0 | 2.2 | 3.2
$(500,5)$ | 2.2 | 1.2 | 1.6 | 1.4 | 3.0 | 2.0 | 1.8 | 2.0

Notes: $A$ and $B$ refer to $(\mu_{1},\mu_{2})=(-1,1)$ and $(-0.5,0.5)$,
respectively, and $C$ and $D$ refer to $(\beta_{1},\beta_{2})=(1,1)$ and
$(-1,1)$, respectively. $E$ and $F$ refer to
$(\sigma_{1},\sigma_{2})=(0.3,0.1)$ and $(0.1,0.1)$, respectively. The mixing
proportion is set to $(\alpha_{1},\alpha_{2})=(0.2,0.8)$. The asymptotic
simulations are based on 500 repetitions.

Table 5: Frequency of the Number of Components with the Simulated Data

$M$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
SHT with EM test | 0 | 0 | 0 | 0.26 | 0.72 | 0.02 | 0
AIC | 0 | 0 | 0 | 0.01 | 0.13 | 0.31 | 0.55
BIC | 0 | 0 | 0 | 0.41 | 0.58 | 0.01 | 0

  * 1

The data are generated using the estimated parameters based on the Chilean
textile industry with five components and panel length $T=3$, where
$(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5})=(0.16076522,0.32454077,0.09025875,0.35478905,0.06964622)$,
$(\mu_{1},\mu_{2},\mu_{3},\mu_{4},\mu_{5})=(-1.241241,-0.33803875,0.4480291,0.52379553,1.4139465)$,
$(\beta_{1},\beta_{2},\beta_{3},\beta_{4},\beta_{5})=(0.451833,-0.05988709,-0.2453261,-0.03106076,0.2053708)$,
and
$(\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4},\sigma_{5})=(0.9933480,0.4585760,0.9954302,0.4116855,0.1863346)$.
We use the panel length and sample size that are equal to those in the
dataset, i.e., $n=196$ and $T=3$.

  * 2

The results are based on 100 repetitions.

  * 3

Each cell indicates the proportion of times that the model selection indicates
an $M$-component model.

##  9 Empirical application

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 Kasahara et al., (2022). 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.

###  9.1 Production function and first-order condition

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

$Y_{it}=\exp{(\epsilon_{it})}F_{t}^{j}(V_{it},L_{it},K_{it},\omega_{it}),$
(35)

with

$\displaystyle
F_{t}^{j}(V_{it},L_{it},K_{it},\omega_{it}):=\exp(\gamma_{t}^{j}+\omega_{it})V_{it}^{\delta_{v,j}}L_{it}^{\delta_{\ell,j}}K_{it}^{\delta_{k,j}},$

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. Kasahara et al., (2022)).

######  Assumption 5.

(a) Each firm belongs to one of $M$ types, and the probability of being type
$j$ is given by $\alpha_{j}=P(D_{i}=j)$ with $\sum_{j=1}^{M}\alpha_{j}=1$. (b)
For the $j^{th}$ type of production technology at time $t$, the output is
expressed in terms of input as in (35), where $\epsilon_{it}\sim
N(0,\sigma_{j}^{2})$ are i.i.d across values of $i$ and $t$. $\omega_{it}$
follows an exogenous first-order stationary Markov process given by
$\omega_{it}=h^{j}(\omega_{it-1})+\eta_{it}$ where, conditional on
$\mathcal{I}_{it-1}$, $\eta_{it}$ is a mean-zero i.i.d. random variable. (c)
$(\gamma_{t}^{j},\omega_{it})\in\mathcal{I}_{it}$ and
$\epsilon_{it}\not\in\mathcal{I}_{it}$.

######  Assumption 6.

(a) Firms are price-takers in both output and input markets, where $P_{Y,t}$
and $P_{V,t}$ are, respectively, the prices of output and intermediate input
in year $t$. (b) $(P_{Y,t},P_{V,t})$ are observed by firms at the beginning of
the period before $V_{it}$ is chosen.

######  Assumption 7.

Values of $V_{it}$ are chosen at time $t$ by maximizing the expected profit
conditional on information $\mathcal{I}_{it}$ at time $t$ and conditional on
the value of $(K_{it},L_{it})$. The profit maximization problem for firms with
type $j$ technology is given by

$V_{it}=\arg\max_{V}P_{Y,t}\mathbb{E}[\exp(\epsilon_{it})|D_{i}=j]F^{j}_{t}(V,K_{it},L_{it},\omega_{it})-P_{V,t}V.$
(36)

In Assumption 5(a), each firm’s production function belongs to one of the $M$
types. Assumption 5(b) assumes that the idiosyncratic productivity shock
follows a normal distribution. Assumption 5(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 6 states that firms
observe input and output prices when deciding on $V_{it}$. Assumption 7
assumes that $V_{it}$ is chosen to maximize the current expected period profit
conditional on the value of $(K_{it},L_{it})$.444We 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 (36)
and the corresponding first-order condition characterize the intermediate
input choice once we interpret $L_{it}$ in (36) as the optimal value chosen by
firm $i$, as discussed in Ackerberg et al., (2015).

Given the above Assumptions 5, 6, and 7, we derive an empirical specification
based on the first-order condition of the profit maximization problem (36),
following the idea developed by Gandhi et al., (2020) and extending it to a
finite mixture production function modeled by Kasahara et al., (2022). 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
(36) together with the production function (35) implies that

$s_{it}=\log\delta_{v,j}+\frac{1}{2}\sigma_{j}^{2}-\epsilon_{it}\quad\text{for
$D_{i}=j$},$ (37)

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 (37), 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 (1):

$f_{M}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M})=\sum_{j=1}^{M}\alpha_{j}\prod_{t=1}^{T}\frac{1}{\sigma_{j}}\phi\left(\frac{s_{it}-\mu_{j}}{\sigma_{j}}\right).$
(38)

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

$f_{M}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M})=\sum_{j=1}^{M}\alpha_{j}\prod_{t=1}^{T}\frac{1}{\sigma_{j}}\phi\left(\frac{s_{it}-\mu_{j}-\beta_{k,j}\log
K_{it}}{\sigma_{j}}\right).$ (39)

In addition, we consider a specification in which we include not only $\log
K_{it}$ but also $\log L_{it}$ as a regressor:

$f_{M}(\boldsymbol{W}_{i};\boldsymbol{\vartheta}_{M})=\sum_{j=1}^{M}\alpha_{j}\prod_{t=1}^{T}\frac{1}{\sigma_{j}}\phi\left(\frac{s_{it}-\mu_{j}-\beta_{k,j}\log
K_{it}-\beta_{\ell,j}\log L_{it}}{\sigma_{j}}\right).$ (40)

###  9.2 Empirical results

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.555Please refer to Kasahara et al., (2021) and Kasahara and Rodrigue,
(2008) 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 6 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.

Table 6: Descriptive Statistics for the Revenue Share of Intermediate
Materials and the Log of Gross Output for the Japanese Firms and Chilean
Plants

Panel A: Japanese publicly traded firms
---
|  |  | $\frac{P_{V,t}V_{it}}{P_{Y,t}Y_{it}}$ | $\log(Y_{it})$
Industry | NObs | n | $mean$ | $sd$ | $mean$ | $sd$
Chemical | 805 | 161 | 0.34 | 0.15 | 17.52 | 1.24
Machine | 790 | 158 | 0.50 | 0.16 | 17.31 | 1.35
Electronics | 775 | 155 | 0.45 | 0.18 | 17.54 | 1.27
Panel B: Chilean plants
|  |  | $\frac{P_{V,t}V_{it}}{P_{Y,t}Y_{it}}$ | $\log(Y_{it})$
Industry | NObs | n | $mean$ | $sd$ | $mean$ | $sd$
Food products | 4645 | 929 | 0.65 | 0.15 | 10.62 | 1.66
Fabricated metal products | 1260 | 252 | 0.53 | 0.18 | 11.00 | 1.37
Textiles | 1130 | 226 | 0.58 | 0.19 | 11.01 | 1.32

  * 1

The summary statistics are based on Japanese firm-level data from 2003 to 2007
and Chilean plant-level data from 1992 to 1996. All observations with
$\log(V_{it}/Y_{it})\leq-3$ and $\log(V_{it}/Y_{it})>\log(2)$ are removed. The
dataset is a balanced panel; i.e., we keep firms/plants that are continuously
observed for these five years.

  * 2

The variable $\frac{P_{V,t}V_{it}}{P_{Y,t}Y_{it}}$ is defined as the revenue
share of the intermediate input, where $P_{V,t}$ is the average price of the
intermediate input at time $t$, $P_{Y,t}$ is the average price of the output,
$V_{it}$ is the quantity of the intermediate input, and $Y_{it}$ is the
quantity of the output.

Table 7: The EM Test for Japanese Producers Without Conditioning Variables

|  | M=1 | M=2 | M=3 | M=4 | M=5
---|---|---|---|---|---|---
|  | $T=3$
Chemical | EM | $436.37^{***}$ | $239.83^{***}$ | $130.1^{***}$ | $126.4^{***}$ | $63.24^{***}$
| BIC | 805.55 | 383.43 | 157.5 | 41.62 | -70.46
Electronics | EM | $563.94^{***}$ | $186.67^{***}$ | $115.82^{***}$ | $81.06^{***}$ | $47.76^{***}$
| BIC | 814.01 | 264.27 | 91.67 | -10.39 | -77.2
Machine | EM | $434.91^{***}$ | $194.48^{***}$ | $72.83^{***}$ | $56.94^{***}$ | $54.77^{***}$
| BIC | 458.72 | 37.85 | -142.28 | -200.74 | -242.71
|  | $T=4$
Chemical | EM | $629.22^{***}$ | $308.6^{***}$ | $181.39^{***}$ | $177.38^{***}$ | $96.35^{***}$
| BIC | 1071.45 | 456.54 | 162.15 | -4.99 | -168.01
Electronics | EM | $803.15^{***}$ | $282.32^{***}$ | $167.83^{***}$ | $106.43^{***}$ | $89.93^{***}$
| BIC | 1081.48 | 292.68 | 24.54 | -484.46 |
Machine | EM | $620.95^{***}$ | $292.52^{***}$ | $118.37^{***}$ | $102.57^{***}$ | $75.32^{***}$
| BIC | 609.1 | 2.14 | -276.04 | -380.16 | -467.96
|  | $T=5$
Chemical | EM | $818.38^{***}$ | $386.08^{***}$ | $219.13^{***}$ | $209.42^{***}$ | $118.25^{***}$
| BIC | 1331.53 | 527.48 | 155.86 | -48.53 | -243.73
Electronics | EM | $1024.86^{***}$ | $375.29^{***}$ | $226.01^{***}$ | $134.53^{***}$ | $126.36^{***}$
| BIC | 1343.12 | 332.61 | -28.32 | -239.31 | -359.17
Machine | EM | $819.98^{***}$ | $389.69^{***}$ | $156.44^{***}$ | $149.98^{***}$ | $96.32^{***}$
| BIC | 775.75 | -30.17 | -406.59 | -548.81 | -683.96

  * 1

The estimation is based on the revenue share of intermediate materials.

  * 2

$~{}^{*}$, $~{}^{**}$, and $~{}^{***}$ indicate significance at the $10\%$,
$5\%$, and $1\%$ levels, respectively.

Table 8: The EM Test for Chilean Producers Without Conditioning Variables

|  | M=1 | M=2 | M=3 | M=4 | M=5
---|---|---|---|---|---|---
|  | $T=3$
Food products | EM | $805.51^{***}$ | $637.77^{***}$ | $204.92^{***}$ | $80.54^{***}$ | $72.41^{***}$
| BIC | 422.55 | -371.13 | -991.96 | -1176.61 | -1236.82
Fabricated metal products | EM | $238.84^{***}$ | $68.91^{***}$ | $26.24^{***}$ | $24.42^{***}$ | $21.82^{***}$
| BIC | 719.74 | 496.49 | 444.02 | 433.01 | 425
Textiles | EM | $229.87^{***}$ | $146.17^{***}$ | $64.76^{***}$ | $27.06^{***}$ | $29.98^{**}$
| BIC | 635.37 | 418.28 | 288.34 | 236.9 | 223.34
|  | $T=4$
Food products | EM | $1165.08^{***}$ | $874.27^{***}$ | $257.49^{***}$ | $130.61^{***}$ | $139.59^{***}$
| BIC | 419.47 | -730.83 | -1586.11 | -1825.87 | -1938.03
Fabricated metal products | EM | $362.1^{***}$ | $120.7^{***}$ | $41.6^{***}$ | $43.68^{***}$ | $20.95^{***}$
| BIC | 905.9 | 559.3 | 453.41 | 427.34 | 399.82
Textiles | EM | $325.17^{***}$ | $222.28^{***}$ | $74.19^{***}$ | $47.58^{***}$ | $51.65^{***}$
| BIC | 821.73 | 510.98 | 303.8 | 243.51 | 210.77
|  | $T=5$
Food products | EM | $1553.9^{***}$ | $1010.31^{***}$ | $290.02^{***}$ | $172.46^{***}$ | $155.25^{***}$
| BIC | 471.66 | -1066.71 | -2057.71 | -2329.38 | -2484.82
Fabricated metal products | EM | $478.94^{***}$ | $176.5^{***}$ | $58.96^{***}$ | $59.37^{***}$ | $33.19^{***}$
| BIC | 1101.11 | 637.21 | 477.1 | 433.62 | 389.54
Textiles | EM | $428.29^{***}$ | $280.46^{***}$ | $103.41^{***}$ | $56.63^{***}$ | $53.57^{***}$
| BIC | 968.16 | 556.01 | 289.55 | 201.41 | 160

  * 1

The estimation is based on the revenue share of intermediate materials.

  * 2

$~{}^{*}$, $~{}^{**}$, and $~{}^{***}$ indicate significance at the $10\%$,
$5\%$, and $1\%$ levels, respectively.

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 (38) and the models with conditioning variables
(39)–(40).

Tables 7 and 8 report the results of the EM test for the model without
conditioning variables (38) 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 (Gandhi et al.,, 2020). 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 (39)–(40).

Table 9 presents the results of the SHT and the BIC when we estimate the
mixture regression model with $\log K_{it}$ in (39) 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 10 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 10 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 10 with those of Table 9, 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.

Table 9: The EM Test and the BIC (Dependent Variable:
$\log\frac{P_{V,t}V_{it}}{P_{Y,t}Y_{it}}$, Regressor: $\log K_{it}$)

$M_{0}$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
Japanese Chemical
EM | $459.4^{***}$ | $236.36^{***}$ | $125.42^{***}$ | $118.36^{***}$ | $87.63^{***}$ | $53.72^{***}$ | $38.69^{***}$ | $34.07^{**}$ | $36.47$ | -
BIC | 1384.76 | 943.61 | 726.53 | 620.32 | 518.86 | 449.92 | 413.49 | 394.46 | 381.48 | 366.09
Japanese Electronics
EM | $560.06^{***}$ | $213.82^{***}$ | $116.29^{***}$ | $78.81^{***}$ | $47.05^{***}$ | $40.77^{***}$ | $27.4^{**}$ | $29.02$ | - | -
BIC | 1332.14 | 788.19 | 593.44 | 495.74 | 434.15 | 406.77 | 385.45 | 372.63 | 367.31 | 351.17
Japanese Machine
EM | $433.19^{***}$ | $202.92^{***}$ | $80.42^{***}$ | $76.82^{***}$ | $53.83^{***}$ | $34.62^{**}$ | $55.65$ | - | - | -
BIC | 1355.6 | 940.49 | 757 | 696.06 | 638.48 | 617.4 | 588.94 | 568.71 | 555.15 | 544.51
Chilean Food Products
EM | $816.06^{***}$ | $489.37^{***}$ | $169.14^{***}$ | $80.88^{***}$ | $80.63^{***}$ | $52.67^{***}$ | $31.29^{***}$ | $17.16^{**}$ | $20.55^{***}$ | $-60.46$
BIC | 6759.39 | 5962.74 | 5499.3 | 5356.47 | 5301.31 | 5241.91 | 5210.27 | 5200.71 | 5210.77 | 5222.29
Chilean Fabricated Metal Products
EM | $199.35^{***}$ | $63.25^{***}$ | $49.24^{***}$ | $30.27^{***}$ | $15.73^{**}$ | $18.25^{**}$ | $10.88$ | - | - | -
BIC | 1923.64 | 1744.72 | 1699.97 | 1670.93 | 1661.03 | 1659.54 | 1665.08 | 1669.02 | 1680.96 | 1695.54
Chilean Textile
EM | $201.86^{***}$ | $95.17^{***}$ | $61.43^{***}$ | $31.17^{***}$ | $14.12^{*}$ | $17.45^{**}$ | $7.94$ | - | - | -
BIC | 1681.91 | 1499.99 | 1424.93 | 1380.93 | 1368.65 | 1364.94 | 1365.72 | 1370.72 | 1382.83 | 1392.24

  * 1

The estimation is based on the revenue share of intermediate materials using
panel data of length $T=3$.

  * 2

$~{}^{*}$, $~{}^{**}$, and $~{}^{***}$ indicate significance at the $10\%$,
$5\%$, and $1\%$ levels, respectively.

Table 10: The EM Test and the BIC (Dependent Variable:
$\log\frac{P_{V,t}V_{it}}{P_{Y,t}Y_{it}}$, Regressors: $\log K_{it}$ and $\log
L_{it}$)

$M_{0}$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
Japanese Chemical
EM | $412.35^{***}$ | $224.09^{***}$ | $141.59^{***}$ | $132.24^{**}$ | $121.56$ | - | - | - | - | -
BIC | 1294.05 | 905.44 | 705.72 | 587.3 | 490.07 | 479.74 | 389.05 | 390.29 | 382.28 | 372.69
Japanese Electronics
EM | $573.11^{***}$ | $218.38^{***}$ | $116.07^{***}$ | $94.76^{**}$ | $47.73$ | - | - | - | - | -
BIC | 1336.69 | 784.95 | 590.73 | 498.55 | 426.23 | 389.91 | 372.05 | 371.64 | 359.15 | 368.25
Japanese Machine
EM | $468.06^{***}$ | $204.01^{***}$ | $93.35^{***}$ | $81.62^{***}$ | $62.00^{***}$ | $37.04^{***}$ | $14.21$ | - | - | -
BIC | 1360.56 | 915.69 | 736.2 | 676.26 | 625.66 | 596.45 | 564.34 | 548.7 | 536.78 | 539.64
Chilean Food Products
EM | $805.09^{***}$ | $478.64^{***}$ | $177.08^{***}$ | $84.13^{***}$ | $80.96^{***}$ | $51.97^{***}$ | $32.3^{**}$ | $19.50$ | - | -
BIC | 6732.11 | 5952.7 | 5506.55 | 5362.27 | 5309.37 | 5257.78 | 5233.37 | 5229.9 | 5242.9 | 5258.41
Chilean Fabricated Metal Products
EM | $204.45^{***}$ | $63.57^{***}$ | $49.42^{***}$ | $28.61^{***}$ | $18.32$ | - | - | - | - | -
BIC | 1926.06 | 1747.29 | 1709.44 | 1685.39 | 1678.71 | 1680.54 | 1685.02 | 1696.19 | 1703.21 | 1723.56
Chilean Textile
EM | $203.69^{***}$ | $90.69^{***}$ | $58.4^{***}$ | $32.55^{***}$ | $16.19$ | - | - | - | - | -
BIC | 1673.99 | 1495.55 | 1431.18 | 1394.54 | 1382.59 | 1373.03 | 1368.8 | 1382.42 | 1394.3 | 1394.09

  * 1

The estimation is based on the revenue share of intermediate materials using
panel data of length $T=3$.

  * 2

$~{}^{*}$, $~{}^{**}$, and $~{}^{***}$ indicate significance at the $10\%$,
$5\%$, and $1\%$ levels, respectively.

##  10 Conclusion

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.

## References

  * Ackerberg et al., (2015) Ackerberg, D. A., Caves, K., and Frazer, G. (2015).  Identification properties of recent production function estimators.  Econometrica, 83(6):2411–2451.
  * Alexandrovich, (2014) Alexandrovich, G. (2014).  A note on the article ‘inference for multivariate normal mixtures’ by j. chen and x. tan.  Journal of Multivariate Analysis, 129:245–248.
  * Ando and Bai, (2016) Ando, T. and Bai, J. (2016).  Panel data models with grouped factor structure under unknown group membership.  Journal of Applied Econometrics, 31(1):163–191.
  * Andrews, (1999) Andrews, D. (1999).  Estimation When a Parameter is on a Boundary.  Econometrica, 67(6):1341–1383.
  * Andrews and Currim, (2003) Andrews, R. L. and Currim, I. S. (2003).  Retention of latent segments in regression-based marketing models.  International Journal of Research in Marketing, 20(4):315–321.
  * Azaïs et al., (2009) Azaïs, J.-M., Gassiat, E., and Mercadier, C. (2009).  The likelihood ratio test for general mixture models with or without structural parameter.  ESAIM: Probability and Statistics, 13:301—327.
  * Balat et al., (2019) Balat, J., Brambilla, I., and Sasaki, Y. (2019).  Heterogeneous firms: Skilled- labor productivity and the destination of exports.
  * Bonhomme and Manresa, (2015) Bonhomme, S. and Manresa, E. (2015).  Grouped patterns of heterogeneity in panel data.  Econometrica, 83(3):1147–1184.
  * Cameron and Heckman, (1998) Cameron, S. V. and Heckman, J. J. (1998).  Life cycle schooling and dynamic selection bias: models and evidence for five cohorts of American males.  Journal of Political Economy, 106(2):262–333.
  * Chen and Chen, (2001) Chen, H. and Chen, J. (2001).  The likelihood ratio test for homogeneity in finite mixture models.  Canadian Journal of Statistics, 29:201–215.
  * Chen and Chen, (2003) Chen, H. and Chen, J. (2003).  Tests for homogeneity in normal mixtures in the presence of a structural parameter.  Statistica Sinica, 13:351–365.
  * Chen et al., (2004) Chen, H., Chen, J., and Kalbfleisch, J. D. (2004).  Testing for a finite mixture model with two components.  Journal of the Royal Statistical Society, Series B, 66:95–115.
  * Chen, (1995) Chen, J. (1995).  Optimal rate of convergence for finite mixture models.  Annals of Statistics, 23(1):221–233.
  * Chen and Li, (2009) Chen, J. and Li, P. (2009).  Hypothesis test for normal mixture models: The EM approach.  Annals of Statistics, 37:2523–2542.
  * Chen and Tan, (2009) Chen, J. and Tan, X. (2009).  Inference for multivariate normal mixtures.  Journal of Multivariate Analysis, 100(7):1367–1383.
  * Chen et al., (2014) Chen, X., Ponomareva, M., and Tamer, E. (2014).  Likelihood inference in some finite mixture models.  Journal of Econometrics, 182(1):87–99.
  * Chernoff and Lander, (1995) Chernoff, H. and Lander, E. (1995).  Asymptotic distribution of the likelihood ratio test that a mixture of two binomials is a single binomial.  Journal of Statistical Planning and Inference, 43:19–40.
  * Dacunha-Castelle and Gassiat, (1999) Dacunha-Castelle, D. and Gassiat, E. (1999).  Testing the order of a model using locally conic parametrization: Population mixtures and stationary ARMA processes.  Annals of Statistics, 27:1178–1209.
  * Deb and Trivedi, (1997) Deb, P. and Trivedi, P. K. (1997).  Demand for medical care by the elderly: a finite mixture approach.  Journal of Applied Econometrics, 12(3):313–336.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–22.
  * Doraszelski and Jaumandreu, (2018) Doraszelski, U. and Jaumandreu, J. (2018).  Measuring the bias of technological change.  Journal of Political Economy, 126(3):1027–1084.
  * Foutz and Srivastava, (1977) Foutz, R. V. and Srivastava, R. C. (1977).  The performance of the likelihood ratio test when the model is incorrect.  The Annals of Statistics, 5(6):1183–1194.
  * Gandhi et al., (2020) Gandhi, A., Navarro, S., and Rivers, D. A. (2020).  On the Identification of Gross Output Production Functions.  Journal of Political Economy, 128(8):2973–3016.
  * Garel, (2001) Garel, B. (2001).  Likelihood ratio test for univariate Gaussian mixture.  Journal of Statistical Planning and Inference, 96:325–350.
  * Garel, (2005) Garel, B. (2005).  Asymptotic theory of the likelihood ratio test for the identification of a mixture.  Journal of Statistical Planning and Inference, 131:271–296.
  * Ghosh and Sen, (1985) Ghosh, J. K. and Sen, P. K. (1985).  On the asymptotic performance of the log-likelihood ratio statistic for the mixture model and related results.  In Le Cam, L. and Olshen, R., editors, Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer, volume 2, pages 789–806. Belmont, CA: Wadsworth.
  * Hao, (2017) Hao, J. (2017).  NormalRegPanelMixture: Finite Mixture Model with Normal Panel Data.  R package version 1.0.
  * Hartigan, (1985) Hartigan, J. (1985).  Failure of log-likelihood ratio test.  In Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer, volume 2, pages 807–810. University of California Press 2. Berkeley.
  * Heckman and Singer, (1984) Heckman, J. and Singer, B. (1984).  A method for minimizing the impact of distributional assumptions in econometric models for duration data.  Econometrica, 52(2):271–320.
  * Kamakura and Russell, (1989) Kamakura, W. and Russell, G. (1989).  A probabilistic choice model for market segmentation and elasticity structure.  Journal of Marketing Research, 26(4):379–390.
  * Kasahara and Rodrigue, (2008) Kasahara, H. and Rodrigue, J. (2008).  Does the use of imported intermediates increase productivity? plant-level evidence.  Journal of Development Economics, 87(1):106–118.
  * Kasahara et al., (2022) Kasahara, H., Schrimpf, P., and Suzuki, M. (2022).  Identification and estimation of production function with unobserved heterogeneity.  Technical report, ESRI Discussion Paper Series No.368.
  * Kasahara and Shimotsu, (2009) Kasahara, H. and Shimotsu, K. (2009).  Nonparametric Identification of Finite Mixture Models of Dynamic Discrete Choices.  Econometrica, 77(1):135–175.
  * Kasahara and Shimotsu, (2012) Kasahara, H. and Shimotsu, K. (2012).  Testing the number of components in finite mixture models.
  * Kasahara and Shimotsu, (2014) Kasahara, H. and Shimotsu, K. (2014).  Non-parametric identification and estimation of the number of components in multivariate mixtures.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 76(1):97–111.
  * Kasahara and Shimotsu, (2015) Kasahara, H. and Shimotsu, K. (2015).  Testing the number of components in normal mixture regression models.  Journal of the American Statistical Association, 110(512):1632–1645.
  * Kasahara and Shimotsu, (2019) Kasahara, H. and Shimotsu, K. (2019).  Testing the Order of Multivariate Normal Mixture Models.