The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
42,451 characters
Poisson Regression under Multivariate Sample Selection
\maketitle
\begin{abstract}
This paper develops a Poisson regression model with multivariate sample selection, in which the outcome is observed only when several potentially correlated selection conditions are satisfied. To the best of our knowledge, this is the first Poisson sample selection model that allows for an arbitrary number of selection equations. We derive the conditional mean of the observed outcome under joint normality of the outcome and selection errors and obtain a multivariate selection-correction term.
We prove identification of the model parameters and show that, under suitable support and rank conditions, the outcome parameters can be identified without an exclusion restriction. We also propose two-step estimation procedures based on nonlinear least squares and Poisson pseudo-maximum likelihood. To the best of our knowledge, this paper is the first to apply PPML to a Poisson regression model with sample selection. The consistency of both estimators is established, and a robust two-step sandwich covariance matrix is proposed to account for the estimation error from the first-step selection model.
In addition, factorial moments are used to recover the variance of the latent outcome error and the correlations between the outcome and selection errors. Monte Carlo simulations show that ignoring sample selection leads to persistent bias when the outcome and selection errors are correlated, while the proposed PPML estimator substantially reduces this bias and is more stable than nonlinear least squares, especially under moderate and strong selection dependence.
\end{abstract}
\section{Introduction}
Count outcomes frequently arise in economics, particularly in fields such as health economics, labour economics, economics of crime, economics of trade among many others \citep{winkelmann2008econometric}. Since count variables are non-negative integers, conventional linear regression models are often inappropriate for their analysis. As a result, econometric models specifically designed for count data have received considerable attention in the literature \citep{winkelmancountd,hausman,haslett2022modelling}.
The Poisson regression model remains the classic specification for the analysis of count data \citep{cameron2013regression}. However, this model assumes that the conditional variance is equal to the conditional mean, which is often inconsistent with empirical applications. Consequently, more flexible specifications such as the negative binomial model have become widely used in applied research \citep{hausman,hilbe}.
In many empirical applications, the outcome variable is observed only for a non-randomly selected subset of the population. Such sample selection mechanisms arise when the probability of observing an outcome depends on observed or unobserved characteristics that also affect the outcome of interest. For example, the number of trips taken by household members may only be observed for households satisfying specific survey participation conditions, such as vehicle ownership requirements.
Ignoring this non-random selection process may therefore lead to biased and inconsistent parameter estimates \citep{terza1998}.
In the patent-production framework \citep{hausman}, observed patent counts are modeled as the outcome of a single stochastic count-generating process. However, the actual innovation mechanism is substantially more complex. Before a patent count becomes observable, a firm must pass through several latent stages, including engagement in R\&D activity, generation of patentable innovations, submission of patent applications, and successful patent approval. These stages are driven by different economic factors and may involve distinct dependence structures. Standard count-data sample selection models typically collapse these mechanisms into a single latent participation equation, implicitly assuming that selection occurs through one binary process only \citep{greene1994,winkelmann1998}. Such a specification may be overly restrictive in empirical environments where observation of count outcomes is generated by multiple sequential or interdependent selection mechanisms. This motivates the development of generalized count-data selection frameworks with multiple latent selection equations.
The classical sample selection model for continuous outcomes was introduced by
\citet{heckman1979}. Heckman proposed a procedure that accounts for non-random
sample selection and allows consistent estimation of model parameters. However,
the original framework considered only a single outcome equation and a binary
selection mechanism. Since then, the sample selection framework has been extended
to more general settings with multiple outcomes, multiple selection equations,
and more complex selection mechanisms
\citep{Poirier, DMF,
Bourguignon,Novi,Tauchmann,DeLuca,Li,Ogundimu,Kim,Kossova2018,Kupriianova2020,Rezaee,Kossova2022}.
Nevertheless, these approaches were primarily developed for continuous outcomes
and cannot be directly applied to count-data settings.
Several approaches have been proposed to incorporate sample selection mechanisms into count-data models. \citet{greene1994} extended Poisson and negative binomial specifications to account for non-random sample selection, excess zeros, and overdispersion in count outcomes. \citet{terza1998} developed count-data models with endogenous switching, sample selection, and endogenous treatment effects, and proposed both fully parametric and partially parametric estimation methods. \citet{miranda2006} further extended maximum likelihood estimation to endogenous switching and sample selection models for count outcomes, which made these models more accessible for applied empirical work. More recently, \citet{wysz2018} proposed a flexible semiparametric copula-based framework for count-data sample selection models, allowing for non-normal dependence between the selection and outcome equations, nonlinear covariate effects, and alternative count distributions. These contributions substantially advanced the modeling of count data under non-random selection.
However, they remain based on the standard two-equation structure, in which a single binary selection equation determines whether the count outcome is observed. Thus, the existing literature does not provide a general framework that allows for an arbitrary number of selection equations. This limitation is important when sample inclusion is generated by several sequential or simultaneous selection mechanisms. The present paper addresses this gap by developing a generalized count-data sample selection model with an arbitrary number of selection equations.
\section{Count data model with multivariate sample selection}
\subsection{Model}
Throughout the rest of the paper, we assume that the errors are independent across observations. The index j denotes the selection equation.
Consider binary selection variables $z_{i,j}$ such that:
\begin{equation*}
z_{i,j} = \begin{cases}
1 ,~w_{i,j}\gamma_j+u_{i,j}\geq0\\
0 ,~\text{otherwise}
\end{cases}=
\begin{cases}
1,~-w_{i,j}\gamma_j\leq u_{i,j}\\
0 ,~\text{otherwise},
\end{cases}
\end{equation*}
where $w_{i,j}$ is a row vector of regressors, $\gamma_j$ is a column vector of coefficients and $u_j$ is a random error.
Denote by $y_i^*$ and $y_i$ a latent dependent (target) variable and its censored value, respectively. Specifically, following Heckman (1979), Kossova \& Potanin (2018). Suppose that $y_i$ is observable only for particular values of this binary variable
\begin{equation*}
\begin{aligned}
y_i
&=
\begin{cases}
y_i^*, & \text{if } z_{i,1}=1,\ldots,z_{i,m}=1,\\
\text{unobserved}, & \text{otherwise},
\end{cases}
\\[0.5em]
&=
\begin{cases}
y_i^*, & \text{if } -\eta_{i,1}\leq u_{i,1},\ldots,-\eta_{i,m}\leq u_{i,m},\\
\text{unobserved}, & \text{otherwise},
\end{cases}
\end{aligned}
\end{equation*}
We will write $z_{i,1}=1,\ldots,z_{i,m}=1$ as $Z_i = 1$ for brevity.
in contrast to previous studies on multivariate sample selection we assume that $y^*_i|~x_i, \varepsilon_i $ is taken from a discrete distribution with $\mathbb{E}[y^*|x_i,\varepsilon]=\text{exp}\{x_i\beta+\varepsilon_i\}$. Throughot the article $x_i$ is a row-vector of regressors, $\beta$ - is a column-vector m$\times$1.
For example, suppose that $y_i$ is the number of cigarettes smoked per day. We observe this variable only if the respondent smokes and answers the question about the number of cigarettes smoked. Thus, $z_{i,1}$ indicates whether the respondent smokes, while $z_{i,2}$ indicates whether the respondent answered the corresponding question.
Let
\(
u_i=(u_{i,1},\ldots,u_{i,m})~\text{and }
\eta_i= (w_{i, 1}\gamma_1, \ldots, w_{i, m}\gamma_m)
\). Then the multiple selection event can be written as
\(
u_i \geq -\eta_i,
\)
where the inequality is understood element by element.
Following \citep{Das2003,Kossova2018} suppose that ($\varepsilon_i$, \ldots, $u_{i,m})'$ are independent of ($x_i$ $w_i$). Also $\mathbb{E}[y_i^*|x_i,w_i,\varepsilon_i,u] = \mathbb{E}[y^*_i|x_i,\varepsilon_i] = \text{exp}\{x_i\beta+\varepsilon_i\}$. Under this assumption, the variables and random errors entering the selection equations do not directly affect the outcome once \(\varepsilon_i\) is fixed. However, \(\varepsilon_i\) and $u_i$ may be correlated.
Then by the law of iterated expectations:
\begin{equation}
\begin{aligned}
&\mathbb{E}\!\left[
y^*_i
\mid
x_i,w_i,Z_i=1
\right]
\\&=
\mathbb{E}\!\left[
\mathbb{E}\!\left[
y^*_i
\mid
x_i,w_i,Z_i=1,
\varepsilon_i,u_i
\right]
\mid
x_i,w_i,Z_i=1
\right]
\\
&=
\mathbb{E}\!\left[
\mathbb{E}\!\left[
y^*_i
\mid
x_i,w_i,
\varepsilon_i,u_i
\right]
\mid
x_i,w_i,Z_i=1
\right]
\\
&=
\mathbb{E}\!\left[
e^{x_i\beta+\varepsilon_i}
\mid
x_i,w_i,Z_i=1
\right]
\\
&=
\mathbb{E}\!\left[
e^{x_i\beta+\varepsilon_i}
\mid
x_i,w_i,u_i\geq -\eta_i
\right] =
e^{x_i\beta}
\mathbb{E}\!\left[
e^{\varepsilon_i}
\mid x_i, w_i,
u_i\geq -\eta_i
\right].
\end{aligned}
\end{equation}
$Z_i=1$ is fully determined by $w_i$ and $u_i$, since $z_{i,j}=1$ if and only if $u_{i,j}\geq -w_{ij}\gamma_j$. Therefore, conditional on $w_i$ and $u_i$, this event gives no additional information.
Notice that $\mathbb{E}\!\left[
e^{\varepsilon_i}
\mid
u_i\geq -\eta_i
\right]$ depends on $-\eta_i$ that depends on $w_i$, so the expectation is not a constant but the function of $\eta_i$. We consider two estimators which address this issue.
Suppose that the joint distribution of \(\varepsilon_i\) and \(u_i\) is multivariate normal:
\[
\begin{pmatrix}
\varepsilon_i \\
u_{i,1} \\
u_{i,2} \\
\vdots \\
u_{i,m}
\end{pmatrix}
\sim
\mathcal{N}
\left[
\begin{pmatrix}
0 \\
0 \\
0 \\
\vdots \\
0
\end{pmatrix},
\begin{pmatrix}
\sigma^2 & \delta_1 & \delta_2 & \cdots & \delta_m \\
\delta_1 & 1 & \rho_{1,2} & \cdots & \rho_{1,m} \\
\delta_2 & \rho_{2,1} & 1 & \cdots & \rho_{2,m} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\delta_m & \rho_{m,1} & \rho_{m,2} & \cdots & 1
\end{pmatrix}
\right].
\]
Here \(\delta=(\delta_1,\ldots,\delta_m)'\) is the covariance vector between the outcome error and the selection errors, that is \(\delta_j=\operatorname{Cov}(\varepsilon_i,u_{i,j})\). Since \(\operatorname{Var}(u_{i,j})=1\), the corresponding correlation is \(\rho_j=\dfrac{\delta_j}{\sigma}\). In the conditional-mean representation below we use \(\delta\), while the separate of \(\sigma\) and \(\rho\) is discussed later. In an ordered probit model, the scale of the latent variable is not identified from the observed ordered outcomes. Therefore, the variance of the latent error is normalized to one, $\operatorname{Var}(u_{i,j})=1$, which fixes the scale of the model; see \citet{Hansen2022}.
For calculation $\mathbb{E}\!\left[
e^{\varepsilon_i}
\mid
u_i\geq -\eta_i
\right]$ we use moment-generating function of multivariate normal truncated distribution \citep{wilhelm2021} which is provided in Appendix~\ref{app:mgf}.
By definition m.g.f of random column-vector $\xi$ is $\mathbb{E}[\text{exp}\{t'\xi\}]$ where t is same space column-vector so for truncated distribution it takes the form $\mathbb{E}[\text{exp}\{t'\xi\}|a\leq\xi\leq b]$, where $a$ and $b$ are vectors of lower and upper truncation points, respectively. In our case there is only bottom truncation, so $b$ goes to $+\infty$, $\varepsilon$ is not truncated, and all other variables are truncated by $-\eta_i$. So when $t=(1 ~0~\cdots ~0)$ and $a = (-\infty ~-\eta_{i,1} ~\ldots ~ -\eta_{i,m} )$ then we may get the conditional expected value, but for the later decomposition, it is useful to write this expression not only for \(r=1\), but for any positive integer \(r\). Taking \(t_r=(r,0,\ldots,0)'\), the moment-generating function of the truncated multivariate normal distribution gives in formula \ref{eq:truncated_mgf}
\begin{equation} \begin{aligned} &\mathbb{E}\left[ \exp\{r\varepsilon_i\} \mid u_{i1}\geq -\eta_{i,1},\ldots,u_{i,m}\geq -\eta_{i,m} \right] \\ &= \frac{\exp\left\{\frac{r^2\sigma^2}{2}\right\}} {\alpha(2\pi)^{\frac{m+1}{2}}|\Sigma_{(m+1)\times(m+1)}|^{\frac{1}{2}}}\times \\ & \int_{-\infty}^{+\infty} \ldots \int_{-\eta_{i,m}-r\delta_m}^{+\infty} \exp\left\{ -\frac{1}{2}x'\Sigma_{(m+1)\times(m+1)}^{-1}x \right\} dx_{u_m}\cdots dx_{\varepsilon} \\ &= \frac{ \exp\left\{\frac{r^2\sigma^2}{2}\right\} \cdot \mathbb{P}\{-(\eta_{i,m}+r\delta_m)<u_{i,m},\ldots, -(\eta_{i,1}+r\delta_1)<u_{i,1}\} } { \mathbb{P}\{-\eta_{i,m}<u_{i,m},\ldots,-\eta_{i,1}<u_{i,1}\} } \\ &= \exp\left\{\frac{r^2\sigma^2}{2}\right\} \cdot \frac{ \Phi_m(\eta_{i,m}+r\delta_m,\ldots,\eta_{i,1}+r\delta_1) }{ \Phi_m(\eta_{i,m},\ldots,\eta_{i,1}) }. \end{aligned} \label{eq:gencor} \end{equation}
where $\alpha$ = $\mathbb{P}\{-\eta_{i,m}<u_{i,m}, \ldots, -\eta_{i,1}<u_{i,1}\}$, $\Phi_m(\cdot)$ is a cumulative distribution function of the $m$-dimensional normal distribution with mean vector $\mathbf{0}$ and covariance matrix $\Sigma_{m\times m}$.
This expression also includes the conditional mean correction as a special case. Indeed, setting \(r=1\) gives the \(\mathbb{E}[\text{exp}\{\varepsilon_i\}\mid x_i,w_i,z_{i,1}=1,\ldots,z_{i,m}=1]\). For higher values of \(r\), the same formula will be used to construct factorial moment identities that will be used for decomposition between correlation of selection and outcome equations and variance of outcome equation.
Finally, using equation~\eqref{eq:gencor}, we obtain the following result:
\begin{equation}
\mathbb{E}[y^*_i \mid x_i,w_i, z_{i,1}=1,\ldots,z_{i,m}=1]
=
e^{x_i\beta+\frac{\sigma^2}{2}}
\cdot
\frac{
\Phi_m(w_{i,m}\gamma_m+\delta_m,\ldots,w_{i,1}\gamma_1+\delta_1)
}{
\Phi_m(w_{i,m}\gamma_m,\ldots,w_{i,1}\gamma_1)
}.
\end{equation}
For brevity, denote the selection-correction term as
\[
\lambda(w_i\gamma,\delta;R)
=
\frac{
\Phi_m(w_{i,m}\gamma_m+\delta_m,\ldots,w_{i,1}\gamma_1+\delta_1;R)
}{
\Phi_m(w_{i,m}\gamma_m,\ldots,w_{i,1}\gamma_1;R)
},
\]
where \(w_i\gamma=(w_{i,1}\gamma_1,\ldots,w_{i,m}\gamma_m)\), \(\delta=(\delta_1,\ldots,\delta_m)\) are regressors from selection equations and the covariance vector between the outcome error, and \(R\) is the correlation matrix of the selection errors.
The obtained results generalize the model of \citet{terza1998}. When there is only one selection equation, the conditional expectation is the same as in the original model.
\subsection{Identification of the generalized count data sample selection model}
~\label{sec:identification}
Under the following sufficient conditions, the outcome parameters are identified without an exclusion restriction:
\begin{itemize}
\item[\textbf{A1.}] The multivariate probit selection subsystem is identified.
\item[\textbf{A2.}] The disturbances are jointly normally distributed as specified above, with $R$ positive definite and $\sigma^2>0$.
\item[\textbf{A3.}] The support of the continuous regressor vector $Z$ contains a nonempty open subset of $\mathbb{R}^k$.
\item[\textbf{A4.}] The matrix of coefficients determining the selection indices has full row rank, $\operatorname{rank}(A)=m$.
\end{itemize}
\noindent
Under A1--A4, if two parameter pairs $(\beta,\delta)$ and
$(\widetilde{\beta},\widetilde{\delta})$ generate the same conditional mean for all admissible values of the regressors, then
$$
\widetilde{\beta}=\beta,
\qquad
\widetilde{\delta}=\delta.
$$
Hence, the outcome parameters are identified without an exclusion restriction.
The identification result follows from the nonlinear form of the multivariate normal selection correction and the support and rank conditions stated above. A formal proof is provided in Appendix ~\ref{app:identification}.
\subsection{Decomposition of structural parameters}
The conditional-mean representation identifies \((\beta,\gamma,R,\delta)\), but it does not separately recover the variance \(\sigma^2\) of the outcome error and the correlations \(\rho_j\). Getting $\rho$ is crucial for economic interpretation as power of connection between unobserved factors that affects both the count equation and selection equatinos. To recover these parameters, we use the normalized second factorial moment of $y_i$ - \(\widehat{A}_2\). Its full derivation is presented in Appendix \ref{app:decomp}, while here we report only the expressions used for calculation.
The estimation of The normalized second factorial moment is defined as
\[
\widehat{A}_2
=
\dfrac{1}{n_s}
\sum_{i=1}^{n}
\mathbb{I}\{Z_i=1\}
\dfrac{(y_i)_2}{\widehat{\mu}_i^2}
\dfrac{\widehat{\lambda}_{1,i}^2}
{\widehat{\lambda}_{2,i}},
\]
where \(n\) is the total number of observations, \(n_s=\sum_{i=1}^{n}\mathbb{I}\{Z_i=1\}\) is the number of selected observations, where $\mathbb{I}$ is indicator function. The second falling factorial is \((y_i)_2=y_i(y_i-1)\), while \(\widehat{\mu}_i\) is the fitted conditional mean. The selection-correction term of order \(r\) is defined as
\[
\widehat{\lambda}_{r,i}
=
\dfrac{
\Phi_m\left(
w_i\widehat{\gamma}
+
r\widehat{\delta};
\widehat{R}
\right)
}{
\Phi_m\left(
w_i\widehat{\gamma};
\widehat{R}
\right)
}.
\]
For the Poisson model, the normalized second factorial moment satisfies \(\widehat{A}_2=\exp\{\widehat{\sigma}^2\}\). Therefore, the variance of the latent outcome error is recovered as \(\widehat{\sigma}^2=\ln\widehat{A}_2\), while its standard deviation is \(\widehat{\sigma}=\sqrt{\ln\widehat{A}_2}\).
Since \(\delta_j=\sigma\rho_j\), the correlations between the outcome error and the selection errors are estimated as \(\widehat{\rho}_j=\dfrac{\widehat{\delta}_j}{\widehat{\sigma}}\).
If the outcome equation contains a constant, the structural constant is recovered as \(\widehat{\beta}^{\,str}_0=\widehat{\beta}_0-\dfrac{\widehat{\sigma}^2}{2}\).
\section{Estimation of count data multiple sample selection generalization}
\subsection{Nonlinear least squares}
The estimation can be implemented in two steps the same way is shown at \citet{terza1998}. At the first step, we estimate
the selection subsystem. Since the selection equations have the form
$z_{i,j}=\mathbb{I}\{w_{i,j}\gamma_j+u_{i,j}\geq 0\}$, $j=1,\ldots,m$ and we suggest that errors are normally distributed, this part of the
model is a standard multivariate probit model. Therefore, the parameters
$\gamma$ and the correlation matrix $R$ can be estimated by maximum
likelihood using all observations.
At the second step, we estimate the parameters of the outcome equation.
Then $\beta$ and $\delta$
can be estimated by minimizing the following nonlinear least squares
criterion:
\[
S(\hat{\beta},\hat{\rho})
=
\arg\min_{\beta,\rho}
\frac{1}{n_s}
\sum_{i=0}^n
\mathbb{I}\{Z_i-1\}\left[
y_i-\exp\{x_i\beta\}\lambda(w_i\hat{\gamma},\delta)
\right]^2,
\]
Since the second step minimizes the sample average of the loss function and the first-step parameters are consistent, the proposed procedure can be treated as a two-step M-estimator in the sense of Section 12 of \citet{wooldridge}. Below, we show that the sufficient conditions for consistency are satisfied in the present model.
Let $\theta=(\beta',\delta')'$ be the vector of second-step parameters. Denote the sample criterion by $S_n(\theta)$ and define the theoretical criterion as
$$
S(\theta)
=
\mathbb{E}\left[
\left(y_i-\exp\{x_i\beta\}\lambda(w_i{\gamma},\delta)\right)^2
\mid
z_{i,1}=\cdots=z_{i,m}=1
\right].
$$
Under the regularity conditions stated by
Theorem 22.2 of \citet{Hansen2022} implies that the sample criterion uniformly converges in probability to the theoretical criterion on $\Theta$, that is
$$
\sup_{\theta\in\Theta}
\left|S_n(\theta)-S(\theta)\right|
\overset{p}{\longrightarrow}
0.
$$
Therefore, asymptotically, minimization of the sample criterion is equivalent to minimization of the theoretical criterion.
By Theorem 2.7 of \citet{Hansen2022}, the theoretical mean squared error is minimized by the conditional expectation function. In the present model this function is $m_0(x_i,w_i)=\mathbb{E}[y_i\mid x_i,w_i,z_{i,1}=\cdots=z_{i,m}=1]=g(x_i,w_i;\theta_0)$. Hence, the theoretical criterion can attain its minimum only if $g(x_i,w_i;\theta)=m_0(x_i,w_i)$ almost surely on the selected sample. Since the model is identified, this equality implies $\theta=\theta_0$. Therefore, $S(\theta)$ has a unique minimum at the true parameter vector $\theta_0$.
Therefore, the proposed estimation procedure can be interpreted as a two-step M-estimator in the sense of \citet{wooldridge}. The first step provides consistent estimates of the selection-equation parameters, while the second step minimizes the nonlinear least squares criterion with generated first-step quantities. Under the standard regularity conditions for two-step M-estimators, $S_n(\theta)$ uniformly converges to $S(\theta)$. Since $S(\theta)$ is uniquely minimized at $\theta_0$, the second-step nonlinear least squares estimator is consistent:
$$
(\hat{\beta},\hat{\delta})
\overset{p}{\longrightarrow}
(\beta_0,\delta_0).
$$
\subsection{Poisson pseudo-maximum likelihood}
As an alternative to nonlinear least squares, the parameters of the outcome equation can be estimated by Poisson pseudo-maximum likelihood (PPML). \citet{santos} show that PPML can be consistently applied when the conditional mean is correctly specified, without requiring the dependent variable to follow a Poisson distribution. We adapt this approach to the present sample selection model. Thus, under standard regularity conditions, consistency of the second-step estimator requires correct specification of the conditional mean rather than the full conditional distribution.
Let \(\theta=(\beta',\delta')'\) denote the vector of second-step parameters, and let \(\eta=(\gamma,R)\) are parameters estimated at the first step. For selected observations, define the conditional mean as \(m_i(\theta,\eta)=\exp\{x_i\beta\}\lambda(w_i\gamma,\delta;R)\).
The population PPML objective function for the selected observations is
\[
Q(\theta,\eta_0)
=
\mathbb{E}\left[
m_i(\theta,\eta_0)
-
y_i\ln m_i(\theta,\eta_0)
\mid Z_i=1
\right].
\]
At the true parameter vector, the corresponding population moment condition
is equal to zero. By the law of iterated expectations,
\[
\begin{aligned}
&\mathbb{E}\left[
\left.
\dfrac{\partial \ln m_i(\theta_0,\eta_0)}{\partial\theta}
\left[
y_i-m_i(\theta_0,\eta_0)
\right]
\right|
Z_i=1
\right]
\\
&=
\mathbb{E}\left[
\left.
\dfrac{\partial \ln m_i(\theta_0,\eta_0)}{\partial\theta}
\underbrace{
\mathbb{E}\left[
y_i-m_i(\theta_0,\eta_0)
\mid x_i,w_i,Z_i=1
\right]
}_{=\,0}
\right|
Z_i=1
\right]
=0,
\end{aligned}
\]
because
\[
m_i(\theta_0,\eta_0)
=
\mathbb{E}[y_i\mid x_i,w_i,Z_i=1].
\]
The sample analogue of the population objective function is
\[
Q_n(\theta,\widehat{\eta})
=
\dfrac{1}{n_s}
\sum_{i=1}^{n}
\mathbb{I}\{Z_i=1\}
\left[
m_i(\theta,\widehat{\eta})
-
y_i\ln m_i(\theta,\widehat{\eta})
\right].
\]
The PPML estimator minimizes this sample objective function. The corresponding
first-order condition is
\[
\dfrac{1}{n_s}
\sum_{i=1}^{n}
\mathbb{I}\{Z_i=1\}
\dfrac{\partial \ln m_i(\theta,\widehat{\eta})}{\partial\theta}
\left[
y_i-m_i(\theta,\widehat{\eta})
\right]
\rightarrow
0.
\]
We have already shown at Paragraph 2.2 that the conditional mean is identified. In particular, under the identification of the multivariate probit selection subsystem, the exclusion restriction, and the full-rank condition for the outcome regressors, \(m_i(\theta,\eta_0)=m_i(\theta_0,\eta_0)\) almost surely implies \(\theta=\theta_0\).
It remains to show that the population PPML criterion is uniquely minimized at the true parameter vector. Define
\[
Q(\theta)
=
\mathbb{E}
\left[
m_i(\theta,\eta_0)
-
y_i\ln m_i(\theta,\eta_0)
\mid
Z_i=1
\right].
\]
Let \(m_{0i}=m_i(\theta_0,\eta_0)=\mathbb{E}[y_i\mid x_i,w_i,Z_i=1]\). Using the law of iterated expectations,
\[
Q(\theta)-Q(\theta_0)
=
\mathbb{E}
\left[
m_i(\theta,\eta_0)-m_{0i}
-
m_{0i}
\ln
\left(
\dfrac{m_i(\theta,\eta_0)}{m_{0i}}
\right)
\Bigm|
Z_i=1
\right].
\]
Define \(a_i(\theta)=\dfrac{m_i(\theta,\eta_0)}{m_{0i}}\). Since both conditional means are strictly positive, \(a_i(\theta)>0\), and therefore
\[
Q(\theta)-Q(\theta_0)
=
\mathbb{E}
\left[
m_{0i}
\left\{
a_i(\theta)-1-\ln a_i(\theta)
\right\}
\mid
Z_i=1
\right].
\]
For every \(a>0\), we have \(a-1-\ln a\geq 0\), with equality if and only if \(a=1\). Hence, \(Q(\theta)\geq Q(\theta_0)\), and equality is possible only if \(m_i(\theta,\eta_0)=m_{0i}\) almost surely. Since the conditional mean is identified, this equality implies \(\theta=\theta_0\). Therefore, the population PPML criterion has a unique global minimum at the true parameter vector.
The PPML estimator is an M-estimator because it minimizes a sample average of the pseudo-likelihood loss function. Since the first-step estimator is consistent and, under the standard regularity conditions for two-step M-estimators
\[
\sup_{\theta\in\Theta}
\left|
Q_n(\theta,\widehat{\eta})
-
Q(\theta,\eta_0)
\right|
\overset{p}{\longrightarrow}
0,
\]
the consistency results for two-step M-estimators in Section~12 of \citet{wooldridge} apply. Since \(Q(\theta)\) is uniquely minimized at \(\theta_0\), the second-step PPML estimator is consistent:
\[
(\widehat{\beta},\widehat{\delta})
\overset{p}{\longrightarrow}
(\beta_0,\delta_0).
\]
\subsection{Estimator of the asymptotic covariance matrix }
Since $(\lambda(w_i\widehat{\gamma},\rho;\widehat{R}))$ is constructed using first-step estimates, the usual nonlinear least squares covariance matrix is not sufficient. It treats $(\widehat{\gamma})$ and $(\widehat{R})$ as fixed, although their estimation error also affects $(\widehat{\beta}) $and $(\widehat{\rho})$. Therefore, to address for the generated regressors we proposed the robust two-step sandwich covariance matrix.
Then the covariance matrix is estimated as
\[
\widehat{\mathrm{Var}}(\widehat{\theta})
=
\frac{1}{n}
\widehat{H}^{-1}
\widehat{\Omega}
(\widehat{H}^{-1})',
\]
where
\[
\widehat{H}
=
\frac{1}{n}
\sum_{i=1}^{n}
\frac{\partial \psi_i(\widehat{\theta})}{\partial \theta'},
\qquad
\widehat{\Omega}
=
\frac{1}{n}
\sum_{i=1}^{n}
\psi_i(\widehat{\theta})\psi_i(\widehat{\theta})'.
\]
The vector $\psi_i(\theta)$ contains the first-order conditions from the first-step likelihood and the second-step nonlinear least squares criterion. The detailed derivation of this covariance matrix is given in Appendix D. For inference on the outcome equation, we use the block of $\widehat{\mathrm{Var}}(\widehat{\theta})$ corresponding to $(\widehat{\beta}',\widehat{\rho}')'$.
The covariance matrix can also be used to test whether the selection correction is necessary by testing the joint null hypothesis \(H_0:\rho_1=\cdots=\rho_m=0\), for example, using a standard Wald test.
\section{Simulation example}
Unlike a conventional Monte Carlo design based on a single fixed parameterization, we allow several nuisance features of the data-generating process to vary across replications. This choice is deliberate. The purpose of the simulation is not to evaluate the estimators at one particular calibration, but to examine whether their relative performance is robust across a broader class of admissible data-generating processes. Similar simulation strategies have been used to evaluate estimator performance across multiple or randomly parameterized DGPs rather than at a single point in the parameter space \citep{li2024local,reuvers2024sparse}. Accordingly, the regression coefficients and the correlation structure of the regressors are allowed to vary across replications, while the main parameters of interest for the comparison, including the sample size and the strength of selection dependence, are controlled across simulation scenarios.
Thus, the reported Monte Carlo measures should be interpreted as average finite-sample performance over the specified distribution of DGPs rather than as performance conditional on a single fixed parameter vector.
\subsection{Monte Carlo design}
We use a Monte Carlo simulation with a Poisson outcome and two selection equations. In each replication, five regressors are generated as \(x_i=(x_{i,1},\ldots,x_{i,5})'\sim N(0,R_x)\). The correlation matrix \(R_x\) is generated separately in each replication using the method of \citet{joe2006}.
The latent outcome is generated from
\[
y_i^*\mid x_i,\varepsilon_i
\sim
\operatorname{Poisson}
\left(
\exp\left\{
\beta_0^{str}
+\beta_1x_{i,1}
+\beta_2x_{i,2}
+\beta_3x_{i,3}
+\varepsilon_i
\right\}
\right).
\]
The coefficients \(\beta_0^{str},\beta_1,\beta_2,\beta_3\) are independently drawn from the uniform distribution \(U[-1,1]\) in each replication.
The selection equations are
\[
z_{i,1}
=
\mathbb{I}
\left\{
\gamma_{1,0}
+\gamma_{1,1}x_{i,1}
+\gamma_{1,4}x_{i,4}
+u_{i,1}
\geq 0
\right\},
\]
and
\[
z_{i,2}
=
\mathbb{I}
\left\{
\gamma_{2,0}
+\gamma_{2,2}x_{i,2}
+\gamma_{2,5}x_{i,5}
+u_{i,2}
\geq 0
\right\}.
\]
The outcome is observed only when \(z_{i,1}=z_{i,2}=1\). The variables \(x_{i,4}\) and \(x_{i,5}\) enter only the selection equations and are used as exclusion restrictions. The intercepts in the selection equations are drawn from \(U[0,1]\), while the remaining nonzero coefficients are drawn from \(U[-1,1]\).
The errors are jointly normally distributed:
\[
\begin{pmatrix}
\varepsilon_i\\
u_{i,1}\\
u_{i,2}
\end{pmatrix}
\sim
\mathcal{N}\left[
\begin{pmatrix}
0\\
0\\
0
\end{pmatrix},
\begin{pmatrix}
\sigma^2 & \sigma r & \sigma \rho\\
\sigma \rho & 1 & 0.6\\
\sigma \rho & 0.6 & 1
\end{pmatrix}
\right].
\]
We set \(\sigma=0.75\) and \(\operatorname{Corr}(u_{i,1},u_{i,2})=0.6\). The correlations between the outcome error and both selection errors are equal:
\[
\operatorname{Corr}(\varepsilon_i,u_{i,1})
=
\operatorname{Corr}(\varepsilon_i,u_{i,2})
=
\rho,
\qquad
\rho\in\{0,0.4,0.8\}.
\]
Therefore, \(\delta_1=\delta_2=\sigma \rho\).
We compare six estimators. The first two are naive Poisson and negative binomial regressions, which ignore sample selection. The next two use one selection equation for the combined indicator \(s_i=\mathbb{I}\{z_{i,1}=z_{i,2}=1\}\) and estimate the second step by NLS or PPML. The last two estimators use two separate selection equations and estimate the second step by NLS or PPML.
Each scenario is repeated \(B=100\) times. Let \(\widehat{\theta}_b\) and \(\theta_b\) denote the estimated and true values of a parameter in replication \(b\). Bias is calculated as
\[
\operatorname{Bias}(\widehat{\theta})
=
\dfrac{1}{B}
\sum_{b=1}^{B}
\left(
\widehat{\theta}_b-\theta_b
\right).
\]
The root mean squared error is
\[
\operatorname{RMSE}(\widehat{\theta})
=
\sqrt{
\dfrac{1}{B}
\sum_{b=1}^{B}
\left(
\widehat{\theta}_b-\theta_b
\right)^2
}.
\]
The root median squared error is
\[
\operatorname{RMedSE}(\widehat{\theta})
=
\sqrt{
\operatorname{median}_{b=1,\ldots,B}
\left\{
\left(
\widehat{\theta}_b-\theta_b
\right)^2
\right\}
}.
\]
The win rate is calculated as the share of common successful replications
in which an estimator has the smallest absolute error for a given parameter.
Formally,
\[
\operatorname{WinRate}_m
=
\frac{1}{B}
\sum_{b=1}^{B}
\frac{
\mathbb{I}\left\{
\left|\hat{\theta}_{m,b}-\theta_b\right|
=
\min_j
\left|\hat{\theta}_{j,b}-\theta_b\right|
\right\}
}{K_b},
\]
where \(K_b\) is the number of estimators
that attain the smallest absolute error in replication \(b\).
Thus, in case of ties, the win is divided equally among the tied estimators.
\subsection{Monte Carlo results}
Figures 1--4 summarize the main Monte Carlo results. The complete numerical values of Bias, RMSE, RMedSE, and win rate are reported in Appendix \ref{app:mntc}.
Figure 1 shows that the importance of the selection correction depends directly on the correlation between the outcome and selection errors. When $\rho=0$, the naive Poisson and negative binomial estimators have almost zero bias. In this case, adding the selection terms does not correct any systematic error and only introduces additional estimation noise. Moreover, including the selection-correction term may introduce additional finite-sample bias. When $\rho$ increases to $0.4$ and $0.8$, the situation changes. the bias of the estimator of the intercept if naive models is approximately $0.15$ for $\rho=0.4$ and $0.28$ for $\rho=0.8$, and it does not decrease as $n$ grows. Thus, this bias is caused by sample selection and cannot be removed by increasing the sample size.
Both PPML selection estimators substantially reduce the bias caused by sample selection. For example, at \(\rho=0.8\), the intercept bias of the two-selection PPML estimator remains close to zero for all sample sizes, while the naive Poisson and negative binomial estimators remain strongly biased. A similar pattern is observed for the slope coefficients, although the bias is smaller than for the intercept. In contrast, the two-selection NLS estimator is less stable: its bias may increase or even change sign as \(n\) grows.
\begin{figure}[H]
\centering
\includegraphics[width=0.72\textwidth]{figures/bias.pdf}
\caption{Bias of the Monte Carlo estimators by sample size and dependence level}
\label{fig:mc_bias}
\end{figure}
Figure 2 gives a clearer comparison of the total estimation error. When $\rho=0$, the naive negative binomial estimator generally has the lowest RMSE. Therefore, estimating a selection correction when no selection dependence is present leads to an efficiency loss. For $\rho=0.4$, the two-selection PPML estimator improves rapidly as the sample size increases. For example, At $n=50{,}000$, it has a lower RMSE than the naive estimators for all four coefficients.
The advantage of the two-selection PPML estimator is strongest when $\rho=0.8$. For the intercept, its RMSE decreases from approximately $0.16$ at $n=5{,}000$ to $0.06$ at $n=50{,}000$. Over the same range, the RMSE of the naive estimators remains close to $0.28$. For the slope coefficients, the RMSE of the two-selection PPML estimator also has stable decreasing and becomes the lowest or nearly the lowest at $n=50{,}000$.
The NLS estimators are considerably less stable. For example, the RMSE of the two-selection NLS estimator for the intercept exceeds one in several scenarios and can increase even when the sample size grows. This instability is especially important for the structural intercept. It is recovered as
$\widehat{\beta}^{\,str}_0=\widehat{\beta}_0-\widehat{\sigma}^2/2$, so estimation error in $\widehat{\sigma}^2$ creates an additional source of error in the final intercept estimate.
One possible explanation for the instability of NLS is the strong heteroskedasticity of the count outcome. Since the variance of the outcome increases with its conditional mean, observations with large counts can have a strong effect on the squared-error criterion. The selection-correction term may further increase this effect by generating large fitted values for some observations. This problem is especially important for the structural intercept, because estimation error in the second factorial moment is additionally transmitted through \(\widehat{\sigma}^2\).
\begin{figure}[H]
\centering
\includegraphics[width=0.72\textwidth]{figures/rmse.pdf}
\caption{RMSE of the Monte Carlo estimators by sample size and dependence level}
\label{fig:mc_rmse}
\end{figure}
Figure 3 shows that the median performance of the NLS estimators is much better than their RMSE performance. For example, at \(\rho=0.8\) the RMedSE of the two-selection NLS intercept decreases from approximately \(0.15\) to \(0.05\) as \(n\) increases, while its RMSE remains very large in some scenarios. The large difference between RMSE and RMedSE means that the NLS results are driven by a limited number of replications with extremely large estimation errors. This pattern is consistent with the instability discussed above: because NLS minimizes squared errors, a small number of observations with large counts or large fitted values may have a strong effect on the estimates. Thus, the poor RMSE performance of NLS appears to be driven mainly by rare but very large estimation errors rather than by systematically poor performance across replications.
The PPML estimators do not show the same degree of instability. Their RMSE and RMedSE decline more regularly with the sample size. At $\rho=0.8$, the RMedSE of the two-selection PPML intercept falls from approximately $0.07$ at $n=5{,}000$ to $0.02$ at $n=50{,}000$. In contrast, the RMedSE of the naive intercept remains close to $0.27$. Thus, the improvement of the PPML estimator is present not only in the mean squared error but also in the typical Monte Carlo replication.
\begin{figure}[H]
\centering
\includegraphics[width=0.72\textwidth]{figures/rmedse.pdf}
\caption{RMedSE of the Monte Carlo estimators by sample size and dependence level}
\label{fig:mc_rmedse}
\end{figure}
Figure 4 confirms these conclusions using the share of replications in which each estimator has the smallest absolute error. When $\rho=0$, the naive negative binomial estimator has the highest win rate for most coefficients. As $\rho$ increases, its win rate falls, while the win rate of the PPML selection estimators rises.
For $\rho=0.4$, the win rate of the two-selection PPML estimator generally increases with $n$. For example, its win rate for $\beta_1$ increases from approximately $12\%$ at $n=5{,}000$ to $37\%$ at $n=50{,}000$. For $\rho=0.8$, the two-selection PPML estimator becomes the most frequent winner for most coefficients. Its win rate reaches approximately $45\%$ for the intercept, $43\%$ for $\beta_2$, and $37\%$ for $\beta_3$ at $n=50{,}000$.
Overall, the simulations show a clear trade-off. When selection dependence is absent, the naive estimators are more efficient. When the dependence is moderate or strong, their bias does not disappear with the sample size. In these cases, the two-selection PPML estimator becomes increasingly preferable as $n$ grows. The NLS version corrects part of the selection bias in typical replications, but its sensitivity to a small number of extreme estimates makes it less reliable than PPML.
\begin{figure}[H]
\centering
\includegraphics[width=0.72\textwidth]{figures/winrate.pdf}
\caption{Win rate of the Monte Carlo estimators by sample size and dependence level}
\label{fig:mc_winrate}
\end{figure}
\section{Conclusion}
This paper proposes a generalization of the count data sample selection model to the case of Poisson regression and multiple selection equations. The Monte Carlo results show that the proposed estimator outperforms both the naive count-data models and the standard single-selection Poisson model when the outcome is affected by more than one non-random selection. When the outcome and selection errors are not correlated, the naive estimators remain more efficient because no selection correction is required. However, when this dependence is moderate or strong, the bias of the naive estimators does not disappear as the sample size increases, while the proposed estimator provides substantially more accurate estimates.
The comparison of the two estimation methods used in the second step, NLS and PPML, gives a clear advantage to PPML. The accuracy of nonlinear least squares improves as the sample size grows, but NLS requires considerably more observations to obtain results comparable to PPML. It is also more sensitive to individual replications with very large estimation errors. In contrast, PPML provides stable and accurate estimates already in moderate samples. Therefore, PPML is the preferred second-step procedure in the considered simulation designs.
Moreover, the proposed approach is not limited to the Poisson distribution. For another discrete outcome distribution, the main estimation procedure remains unchanged, while only the moment condition used to recover the variance and covariance components must be adjusted.
The proposed method also provides information that is not available from naive count-data models. In addition to the regression coefficients, it allows the researcher to estimate the variances, covariances, and correlations between the errors of the selection equations and the error of the outcome equation. The estimated correlations show whether unobserved factors that increase the probability of selection also increase or decrease the expected outcome. Therefore, the proposed model allows the researcher to measure both the direction and the strength of selection on unobservables. Overall, the PPML-based estimator corrects the regression coefficients for multiple non-random selection processes and, at the same time, measures the dependence between these processes and the outcome equation. This makes it a useful tool for empirical applications in which the outcome is observed only after several selection decisions and gives researchers more tools for economic processes modeling.
A natural extension of the framework is to consider endogenous switching models. Terza (1998) studies count-data models with a binary switch generated by a single latent equation, while the multiple-selection structure developed here can be extended to settings in which switching depends on several correlated latent processes. Future research may also consider PPML estimation in such models, including PPML for regime-specific conditional means.
\bibliographystyle{apalike}
\bibliography{references}
\pagebreak{}