Estimation of Multivariate Sample Selection Models via a Parameter-Expanded Monte Carlo EM Algorithm

This paper develops a parameter-expanded Monte Carlo EM (PX-MCEM) algorithm to perform maximum likelihood estimation in a multivariate sample selection model. In contrast to the current methods of estimation, the proposed algorithm does not directly depend on the observed-data likelihood, the evaluation of which requires intractable multivariate integrations over normal densities. Moreover, the algorithm is simple to implement and involves only quantities that are easy to simulate or have closed form expressions.


Introduction
Sample selection models, pioneered in [1]- [3], are indispensable to researchers who use observational data for statistical inference.Among the many variants of these types of models, there is a growing interest in multivariate sample selection models.These are used to model a system of two or more seemingly unrelated equations, where the outcome variable for each equation may be non-randomly missing or censored according to its own stochastic selection variable.Applications range from modeling systems of demand equations [4] [5] to household vehicle usage [6]- [8].A common specification is to assume a correlated multivariate normal distribution underlying both the outcomes of interest and the latent variables in the system.
There are two dominant approaches in the current literature to estimate these models.One approach is to use maximum likelihood (ML) estimation.However, as noted in the literature, a major hurdle in evaluating the likelihood is that it requires computations of multivariate integrals over normal densities, which do not generally have closed form solutions. [9] discusses the ML estimation of these models and proposes to use the popular Geweke, Hajivassiliou, and Keane (GHK) algorithm to approximate these integrals in a simulated ML framework.While this strategy works reasonably well, the GHK algorithm can be difficult to implement.Another popular approach is to use two-step estimation (see [10] for a survey).In general, there is a tradeoff in the statistical properties and the computational simplicity for these estimators.If efficiency and consistency are of primary concern, then ML estimation should be preferred over two-step estimation.The objective of this paper is to develop a simple ML estimation algorithm for a commonly used multivariate sample selection model.In particular, this paper develops a parameter-expanded Monte Carlo expectation maximization (PX-MCEM) algorithm that differs from [9] in a few important ways.First, the PX-MCEM algorithm does not use the observed-data likelihood directly, so it avoids the aforementioned integrations.Second, the proposed iterative algorithm does not require the evaluations of gradients or Hessians, which become increasingly difficult to evaluate with more parameters and equations.Third, the algorithm is straightforward to implement.It only depends on quantities that are either easy to simulate or have closed form expressions.This last point is especially appealing when estimating the covariance matrix parameter since there are non-standard restrictions imposed onto it for identification.
This paper is organized as follows.The multivariate sample selection model (MSSM) is formulated in Section 2. Section 3 begins with a brief overview of the EM algorithm for the MSSM and continues with the development of the PX-MCEM algorithm.Methods to obtain the standard errors are discussed.Section 4 offers some concluding remarks.

Multivariate Sample Selection Model
The MSSM is , , , , , where the prime symbol in i s , i s * , and in the rest of this paper is used to denote matrix transpose.Furthermore, , i j x and , i j w are column vectors of exogenous covariates, and j β and j γ are conforming vectors of parameters.Define ( ) , , , w must contain at least one exogenous covariate that does not overlap with , i j x (refer to [11] for these exclusion res- trictions).The unobserved errors ( ) , , , , , , are jointly distributed as a J 2 -dimensional multivariate normal with a mean vector of zeros and an unknown covariance matrix of Ω .
Formally, ( ) The submatrix νν Ω is restricted to be in correlation form to identify the parameters corresponding to the latent variables [9].The other elements of Ω are restricted such that the matrix is symmetric and positive definite.
The covariates and binary selection variables are always observed.Without loss of generality, assume that the outcomes for any observation i are only missing for the first i m equations, where 0 , denote the observed data.The observed-data likelihood derived from (1) through ( 5) is denoted as ( ) obs , , f y β γ Ω .See [9] for an exact expression of this likelihood.

Overview of the EM Algorithm
The PX-MCEM algorithm is based on the EM algorithm of [12].The basic idea behind the EM algorithm is to first augment obs y with a set of "missing data" mis y such that the observed-data likelihood is preserved when the missing data are integrated out of the complete-data likelihood.Formally, the missing data must satisfy ( ) ( ) , , , f y y β γ Ω is the complete-data likelihood to be defined later.
The EM algorithm then proceeds iteratively between an expectation step (E-step) and a maximization step (M-step) as follows.In iteration ( ) where the expectation is taken with respect to the conditional predictive distribution for the missing data, y y π β γ Ω , and in the M-step, find arg max , , , , Denote the maximal values as ( ) γ + , and ( ) Ω , and continue on with the algorithm until convergence.The final maximal values are at least local maxima of the observed-data likelihood function.
For the MSSM, mis y consists of all the missing outcomes and latent variables.Specifically, , Equation ( 10) is a degenerate density since conditioning on i s * in ,com i y determines i s from (3).Note that the observed-data likelihood from [9] is obtained when mis y is integrated out of ( 9), hence the condition in (6) holds.

PX-MCEM Algorithm
The standard EM algorithm using ( 7) and ( 8) is difficult to implement for the MSSM as the E-step and M-step are intractable.The PX-MCEM algorithm addresses this issue by modifying the E-step in two ways and leads to an M-step that can be evaluated with closed form quantities. Stated succinctly, the PX-MCEM algorithm is as follows.
Each step is described in more detail in the subsequent sections.

PX-MC E-Step
Following [13], the first modification is to expand the parameter space of the complete-data likelihood function from ( ) , , , , α δ Σ .The expanded parameters play similar roles as the original parameters, however Σ is expanded into a standard covariance matrix without the correlation restrictions.The parameter-expanded complete-data likelihood function is , , , , with ( ) ( ) , , , are defined analogously to β and γ .The advantage of using ( 12) instead of ( 9) is that Σ is easier to work with in the PX-MC M-step.
Second, instead of computing analytically, it is approximated as (11) with Monte Carlo methods and Gibbs sampling.To draw from y y π β γ Ω , simply draw ,mis i y and i s * from the conditional distribution 9), we have that . For the missing outcomes, it is easy to see from ( , where , Similarly, for the latent variables, The Gibbs sampler recursively samples from the full conditional distributions in ( 14) and (15) in the usual way.After a sufficient burn-in period, the last G draws are used in (11).

PX-MC M-Step and Reduction
Step By recognizing that ( 11) is proportional to the log-likelihood function of a seemingly unrelated regression model with NG observations and 2J equations, the maximization can be performed with IGLS.IGLS utilizes the quantities Σ  removed, which amounts to estimating Θ equation by equation, and then evaluate (17) based on Θ  .Proceed by iterating ( 16) and (17) recursively until convergence.Denote the converged values as ( ) δ + , and ( ) In the reduction step, set ( ) ( ) ( ) diagonal matrix with the first J diagonals equal to 1 and the remaining J diagonals equal to the square root of the last J diagonals of ( ) 1 t + Σ .The previous transformations are referred to as the reduction functions, and they are needed because (12) is used instead of (9) in the algorithm [13].

Standard Errors
The observed information matrix is where ( ) Ξ , and Ξ is a column vector denoting the unique elements in Ω .Evaluate (18) at the ML estimates, and take the expectation and variance with respect to ( ) mis obs ˆ, , , y y π β γ Ω .These moments are estimated by taking additional draws from the Gibbs sampler and constructing their Monte Carlo analogs.The standard errors are the square roots of the diagonals of the inverse estimated quantity in (18).

Concluding Remarks
A new and simple ML estimation algorithm is developed for multivariate sample selection models.Roughly speaking, the implementation of this algorithm only involves iteratively drawing sets of missing data from wellknown distributions and using IGLS on the complete data, both of which are inexpensive to perform.By using parameter expansion and Monte Carlo methods, the algorithm only depends on quantities with closed form expressions, even when estimating the covariance matrix parameter with correlation restrictions.This algorithm is readily extendable to other types of selection models, including extensions to various types of outcome and selection variables with an underlying normal structure, and modifications to time-series or panel data.
respectively the conditional mean and variance of , i j y * given all other elements in ,com i y , i X as a block-diagonal matrix with the rows of covariates corresponding to the elements of i X θ and covariance Ω , and and the number of Gibbs sampling draws G .