Empirical Economics (2023) 64:3167–3195
https://doi.org/10.1007/s00181-023-02406-w

DS-HECK: double-lasso estimation of Heckman selection
model
Masayuki Hirukawa1 · Di Liu2 · Irina Murtazashvili3 · Artem Prokhorov4,5,6Received: 23 July 2022 / Accepted: 2 March 2023 / Published online: 27 March 2023
© The Author(s) 2023

Abstract
We extend the Heckman (1979) sample selection model by allowing for a large numberof controls that are selected using lasso under a sparsity scenario. The standard lasso
estimation is known to under-select causing an omitted variable bias in addition to thesample selection bias. We outline the required adjustments needed to restore consistency of lasso-based estimation and inference for vector-valued parameters of interestin such models. The adjustments include double lasso for both the selection equationand main equation and a correction of the variance matrix. We also connect the estimator with results on redundancy of moment conditions. We demonstrate the effectof the adjustments using simulations and we investigate the determinants of female
labor market participation and earnings in the US using the new approach. The paper
comes with dsheckman, a dedicated Stata command for estimating double-selection
Heckman models.
Keywords Heckman · Probit · Double lasso · Post selection inference

B Artem Prokhorov

artem.prokhorov@sydney.edu.au
Masayuki Hirukawa
hirukawa@econ.ryukoku.ac.jp
Di Liuflyingliudi@gmail.com
Irina Murtazashvili
im99@drexel.edu

1Faculty of Economics, Ryukoku University, Kyoto, Japan

2

Stata Corp, College Station, Texas, USA

3

LeBow College of Business, Drexel University, Philadelphia, USA

4

Discipline of Business Analytics, University of Sydney Business School, Sydney, Australia

5

CEBA, St.Petersburg State University, St.Petersburg, Russia

6

CIREQ, University of Montreal, Montreal, Canada

123

3168

M. Hirukawa et al.

JEL Classification C13

1 Introduction
In this paper we consider an extension of the familiar Heckman (1979) sample selection
model, which can be written as follows:
α + u 1i ,y1i = x1iy2i = I(xi β + zi η + u 2i ≥ 0)

(1)
(2)

where I(·) in the selection Eq. (2) denotes the indicator function, which takes the value
one if its argument is true and zero otherwise, (u 1i , u 2i ) are the error terms and theoutcome variable y1i is observed only if the selection variable y2i = 1. The main
Eq. (1) contains a k1 × 1 vector of explanatory variables x1i and we are interested
in estimating and testing the coefficient vector α. The explanatory variables of theselection equation are separated into two parts, xi = (x1i , x2i ) and zi . In the traditionalversion of the model there is no distinction between x2i and zi – both represent the
explanatory variables that are present in the selectivity model but not in the mainequation of interest. These are the well known exclusion restrictions that facilitate
identification of α. In our setting, for reasons that will become clear shortly we wish
to differentiate between x2i and zi .
Our extension is to introduce a high-dimensional vector of the explanatory variablesin the selectivity model (2) which may or may not belong to the model. The vector x isa low-dimensional k × 1 vector of selection determinants that we wish to keep in themodel no matter what. The vector z is a high-dimensional p × 1 vector of potential
controls, where p can be as large as the (pre-selection) sample size N or larger andwhere we do not know which of these controls are important, if any. The vector β isa k × 1 vector of coefficients on x, which can be a target of inference too. The vector
η, on the contrary, is just a p × 1 nuisance parameter vector.
This extension has many empirical applications in economics where we have awell defined list of regressors for the main equation which has roots in economic
theory (e.g., consumer and labor theory) while what determines selection into thesample is less certain (see, e.g., Roy 1951; Heckman and Honore 1990). The classic
examples are the estimation of the female labor supply function and wage functions(see, e.g., Heckman 1979; Arellano and Bonhomme 2017), which may be subject
to selection bias as determinants of the sample selection are confounded with the
behavioral functions of interest. We return to women’s labor force participation andlabor supply decisions in our empirical application section.
Our objective is to consistently estimate α in the outcome Eq. (1) under a potential
sample selection bias arising from the fact that in the observed sample
α + E(u 1i |xi , zi , y2i = 1) = x1i
α,E(y1i |xi , zi , y2i = 1) = x1i

unless E(u 1i |xi , zi , y2i = 1) = 0, which is a questionable assumption in practice. Heckman (1979) assumed joint normality of (u 1i , u 2i ) and showed that123

DS-HECK: double-lasso estimation of Heckman selection model

3169

E(u 1i |xi , zi , y2i = 1) = γ λ(xi β + zi η), where λ(·) = φ(·)/(·) is known as theinverse Mills ratio. The two-step heckit procedure is (a) to run the maximum likelihood estimation (MLE) for the probit of y2i on (xi , zi ) and use the estimates (β̂, η̂)to obtain λ̂i ≡ λ(xi β̂ + zi η̂) and then (b) to regress y1i on x1i and λ̂i . Under correct
specification, the resulting estimators α̂ and γ̂ are consistent and the usual t-test onγ̂ can be used to test for selection bias. If the null of no bias is rejected, the standard
errors of the second step have to be corrected for the first step estimation error which
is done via a full MLE using normality of the errors or via an analytic correction to
the variance in the second step.
The high-dimensionality of zi poses a challenge in applying the traditional twostep procedure. First, we cannot include all the variables in zi in the first step because
there are too many variables. If p is larger than N , the probit with all xi and zi isinfeasible, and even if p is substantially smaller than N but is large then including allthese variables can cause difficulties in MLE convergence.
In order to make estimation feasible, it is common to impose a certain structure on
η, known in the literature on regularized estimation as a sparsity scenario. In particular,
we assume that only a few elements in the coefficient vector η are substantially different
from zero. Although we assume that η is sparse, we do not know which elements are
non-zero and a consistent model selection technique is required. A popular approachto regularizing linear models is the least absolute shrinkage and selection operator(lasso) developed by Tibshirani (1996). The method penalizes the objective function
with an l1 -norm of the coefficients. This shrinks the irrelevant coefficients to zero andthus serves as a model selection tool. However, even for purely linear models, thisapproach has well known challenges.
First, lasso makes mistakes. Failure to account for the fact that the covariates havebeen selected by lasso results in invalid inference. The reason is that lasso, like manyother model selection techniques, does not always find all the relevant covariates
especially when some coefficients are small. Model selection mistakes made by lasso
cause the distribution of this naive estimator to be biased and nonnormal. For example,
Leeb and Pötscher (2008a, b) and Pötscher and Leeb (2009) show that the normal
approximation for the naive lasso estimator will produce misleading inference. Belloni
et al. (2014b), Belloni et al. (2016b), and Chernozhukov et al. (2018) derive estimatorsthat are robust to the mistakes made by lasso. Such robust estimators are often referred
to as Neyman orthogonal (NO) estimators because they can be viewed as extensionsof an approach proposed by Neyman (1959).
The second challenge is choosing the lasso tuning parameter. Lasso’s ability toselect relevant covariates depends on the method used to choose the tuning parameters.
Belloni et al. (2014b) propose a plug-in method and show that NO estimators perform
well on linear models under that method. Belloni et al. (2016b) extend the linear
lasso to logit models and show good performance using a simplified version of theplug-in method. Drukker and Liu (2022) extend the plug-in method to cross-sectional
generalized linear models and provide Monte Carlo evidence that their extension works
well in finite samples.
In this paper, we develop NO estimation for the model in (1)–(2) which we calldouble-selection Heckman procedure, or DS-HECK. The DS-HECK estimator draws
upon the classical two-step heckit estimator and the double-selection lasso for the high-

123

3170

M. Hirukawa et al.

dimensional generalized linear models proposed by Belloni et al. (2016b). We detailthe steps involved in the estimation, work out the estimator properties and derive thevariance corrections. We also provide new insights into how NO estimation is linked
to results on redundancy of knowledge in moment-based estimation considered by
Breusch et al. (1999) and Prokhorov and Schmidt (2009).
The rest of the paper is organized as follows. Section 2 describes and studies the
DS-HECK estimator. In Sect. 3, we present simulation results that demonstrate an
excellent performance of DS-HECK in finite samples. In Sect. 4, we apply DS-HECK
to estimate married women’s wage using the 2013 PSID wave, in the presence of highdimensional controls and potential sample selection bias. Finally, Sect. 5 concludes.

2 The DS-HECK estimator
2.1 Settings
We maintain the standard assumption of the Heckman sample selection model.
Assumption 1 (a) (x, z, y2 ) are always observed, y1 is observed only when y2 = 1; (b)(u 1 , u 2 ) is independent of x and z with zero mean; (c) u 2 ∼ N (0, 1); (d) E(u 1 |u 2 ) =γ u2.
Assumption 1 is in essence the same as in Wooldridge (2010, p. 803). Part (a)describes the nature of sample selection. Part (b) assumes that x and z are exogenous.
Part (c) is restrictive but needed to derive the conditional expectation of y1 given that
it is observed. Part (d) requires linearity in the conditional expectation of u 1 given u 2 ,and it holds when (u 1 , u 2 ) is bivariate normal. However, it also holds under weakerassumptions when u 1 is not normally distributed.
Additionally, we impose a sparsity scenario on η.
Assumption 2 η is sparse; that is, most of the elements of η are zeros. Namely, ||η||0 ≤s, where || · ||0 denotes the number of non-zero components of a vector. We require s
22p,N ))
to be small relative to the sample size N . In particular, s log (max(−→ 0.
NThis assumption follows Belloni et al. (2016b). In the settings of generalized linear
models, it allows for the estimation of the nuisance parameter in the selection equation
at the rate o(N −1/4 ) (see their Condition IR). In our settings, this rate is needed toguarantee the consistent estimation of β.
Under Assumption 1, it is easy to show thatE(y1 |x, z, y2 = 1) = x1 α + γ λ(x β + z η),(3)

where λ(·) = φ(·)/(·) is the inverse Mills ratio. In essence, this is the classic
formulation of Heckman (1979) where the presence of x2 and z in the selection equation(but not in the main equation of interest) permits estimation of the model even when
the inverse Mills ratio is close to being linear in its argument.

123

DS-HECK: double-lasso estimation of Heckman selection model

3171

We wish to explore the behavior of this conditional expectation with respect topotential errors in the choice of z. It is easy to see from applying the mean-value
theorem to the inverse Mills ratio evaluated at x β, that we can rewrite (3) as follows:


(1)

E(y1 |x, z, y2 = 1) = x1 α + γ λ(x β) + γ λ (q)z η= x1 α + γ λ(x β) + z ω,(4)where q is a point between x β +z η and x β, λ(1) (·) is the first-order derivative of λ(·),
and ω = γ λ(1) (q)η. It is well known that λ(1) (·) is monotone and bounded between-1 and 0 (see, e.g., Sampford 1953). We note that ω depends on q, γ and η, and thatif η is sparse then ω is sparse too.
Proposition 1 The vector ω in Eq. (4) inherits the same sparsity properties as thevector η, i.e. ||ω||0 = ||η||0 .
Proof. Sketches of the proofs of all less obvious propositions are given in the Appendix.
This Proposition makes it clear that a Heckman model with sparsity in the selection
equation can be written as a heckit model with the same sparsity scenario in the mainequation of interest.
Next we derive some conditions on the linear approximation of the inverse Mills
ratio using z in the selected sample which is common for lasso-based model selection
but new in the context of the Heckman model.
Let n denote the size of the selected sample, defined as follows:n=NI(y2i = 1).

i=1Then, for the selected observations, we can write
y1i = x1i
α + g(xi , zi ) + i ,where g(xi , zi ) = γ λ(xi β + zi η), E( i |xi , zi ) = 0 and V( i |xi , zi ) = σ 2 . We follow(Belloni et al. 2014b) and write g(xi , zi ) in a linear form subject to a bound on theapproximation error:
y1i = x1i
α + γ λ(xi β) + zi δ + ri + i ,

(5)
 
n
1

where ri , i = 1, . . . , n, is the approximation error such that n i=1 Eri2 = sOn . Additionally, we assume that the selected and pre-selection sample sizes
are of the same order.
Assumption 3 Equation (5) holds in the selected sample with n : Nn → c ∈ (0, 1)
and with

 n
1 
s
.
Eri2 = O
n
n
i=1

123

3172

M. Hirukawa et al.

The assumption on the approximation error follows (Belloni et al. 2014b). Similar toAssumptions 2, 3 ensures that we can estimate the nuisance parameter in the selected
sample at the rate o(n −1/4 ) rate. In the context of the Heckman model, it implies
that ||δ||0 = ||ω||0 = ||η||0 ≤ s and that δ is also estimated at o(N −1/4 ) sincen and N are of the same order. Example-specific primitive conditions that ensure
Assumption 3 holds are discussed by Belloni et al. (2014b, Section 4.1) for parametricand nonparametric cases, with the parametric example in Section 4.1.1 being most
relevant for our setting.
Next we investigate how to consistently estimate this model accounting for thehigh-dimensional nuisance parameter in both equations.
2.2 Estimation of the selection equation
Clearly if we knew the true value of β, we could treat λ(x β) as a known variable
and we could estimate α and γ , treating δ as nuisance parameters. So we start witha consistent estimation of β in Eq. (2) using the approach of Belloni et al. (2016b)
combined with parameter tuning of Drukker and Liu (2022).
The estimation involves three steps:
Step 1 (post-lasso probit) We start by estimating a penalized probit of y2 on x andz using the lasso penalty:(β̂, η̂) = arg min E N (β,ηi (β, η)) + λ1 ||(β, η)||1 ,where E N denotes the sample mean of N observations, i (·) is the negative loglikelihood for the probit model, || · ||1 is the lasso (l1 ) norm of the parametersand λ1 is a tuning parameter chosen using the plug-in method of Drukker andLiu (2022). This produces a subset of the variables in z indexed by suppor t(η̂),
where for a p-vector v, suppor t(v) := { j ∈ {1, ..., p} : v j = 0}. Thesevariables are used in the post-lasso probit:(β̃, η̃) = arg min E N (β,ηi (β, η)) : suppor t(η) ⊆ suppor t(η̂)

As a result, we obtain the sparse probit estimates (β̃, η̃) where η̃ contains only a
few non-zero elements. Belloni et al. (2016b) propose using these estimates toconstruct weights fˆi = ŵi /σ̂i , where ŵi = φ(xi β̃ + zi η̃), and σ̂i2 = (xi β̃ +zi η̃)(1 − (xi β̃ + zi η̃)), for i = 1, . . . , N .
Step 2. We use the weights from Step 1, to run a weighted lasso regression inwhich for each variable x j in x, j = 1, . . . , k, we run the penalized regressionof fˆi xi j on fˆi zi ,θ̂ j = arg min E N ( fˆi2 (xi j − zi θ j )2 ) + λ2 ||θ j ||1 ,θj123

DS-HECK: double-lasso estimation of Heckman selection model

3173

where λ2 is chosen by the plug-in method of Drukker and Liu (2022). For
each element of x, this produces a selection from the variables in z indexed by
suppor t(θ̂ j ), j = 1, . . . , k.
Step 3 (double-selection probit). We use the variables selected from z in Steps 1
and 2 to run the probit of y2 on x and the union of the sets of variables selectedin Steps 1 and 2:(β̌, η̌) = arg min E N (β,ηi (β, η) fˆi /σ̂i ),where suppor t(η) ⊆ suppor t(η̂) ∪ suppor t(θ̂1 ) ∪ . . . ∪ suppor t(θ̂k ).
In the setting of generalized linear models, Belloni et al. (2016b) show that the
double-selection probit corrects for the omitted variable bias introduced by a naive
application of lasso to Eq. (2). The intuition is that the double selection using weights
fˆi “neymanizes” Step 3. That is, it ensures that the estimation error from the first step
does not affect the estimated parameter vector of the last step.
It follows from Belloni et al. (2016b, Theorem1) that, under Assumption 2, β̌ is aconsistent estimator of β and its variance can be obtained from Step 3 using the wellknown “sandwich” formula for probit. For example,√ in Stata it can be obtained using
the vce(robust) syntax. We obtain a regular √ N -consistent estimator of β withstandard inference, even though a penalized non- N estimator is being used to carryout model selection for the high-dimensional nuisance parameter η.
2.3 Connection to redundancy of moment conditions
Belloni et al. (2016b) use the Neyman orthogonalization to obtain their result. In thissection we show how NO argument relates to the concept of moment redundancy
pioneered by Breusch et al. (1999). This offers an alternative way of arriving at theweights derived by Belloni et al. (2016b).
The key insight of Belloni et al. (2016b) is that the weights f i ensure the validityof the main moment condition:
Egi (β0 , η0 , θ0 ) ≡ E[y2i − (xi β0 + zi η0 )]( f i xi − f i zi θ0 ) = 0,

(6)which has to hold simultaneously with the conditionE∂gi (β0 , η0 , θ0 ) = 0.
∂ηIt is easy to see that Eq. (6) holds due to E(y2i |xi , zi ) = (xi β0 + zi η0 ) and that forfi =φ(xi β0 +zi η0 ), where σi2 = V(y2i |xi , zi ) = (xi β0 + zi η0 )(1 − (xi β0 + zi η0 )),σi2

the zero expected derivative condition holds, too. See Eq. (28) of Belloni et al. (2016b).
What has not been noted is that these conditions correspond to what Prokhorov
and Schmidt (2009) call moment and parameter redundancy (M/P-redundancy), that

123

3174

M. Hirukawa et al.

is the situation when neither the knowledge of the additional moment conditions nor
the knowledge of the parameter they identify help improve efficiency of estimation.
To see this, let xi be a scalar for notational simplicity, and write the moment conditions identifying β and η as follows:
Eh 1i (β0 , η0 ) ≡ E(y2i − (xi β0 + zi η0 )) f i xi = 0
Eh 2i (β0 , η0 ) ≡ E(y2i − (xi β0 + zi η0 )) f i zi = 0

(7)
(8)

where the subscript “0” on a parameter denotes the true value. These moment conditions correspond to the first order conditions of the probit and stem from thespecification P(y2i = 1|xi , zi ) = (xi β0 + zi η0 ).
This is the system of moment conditions considered by Breusch et al. (1999) in
the Generalized Method of Moments (GMM) framework. See their Eq. (6). Theyshow that the (optimal) GMM estimation based on Eqs. (7)–(8) is equivalent to theestimation based on Eq. (8) and the error in the linear projection of Eq. (7) on Eq. (8).
Using their notation, we can write the equivalent moment system as follows:E[h 1i (β0 , η0 ) − 12 −1
22 h 2i (β0 , η0 )] = 0(9)Eh 2i (β0 , η0 ) = 0(10)where 12 and 22 are the relevant parts of the moment variance matrix11 12h 1i (β0 , η0 )=
≡V.
21 22h 2i (β0 , η0 )It is easy to see that Eq. (9) coincides with Eq. (6) subject to the additional notation
2 2

2 2
 −1
that θ0 = 12 −1
22 = Eσi f i xi zi [Eσi f i zi zi ] . It is also easy to see that the entire
estimation problem can be written in the GMM framework as follows:
Egi (β0 , η0 , θ0 ) ≡ E[y2i − (xi β0 + zi η0 )]( f i xi − f i zi θ0 ) = 0
Eh 2i (β0 , η0 ) ≡ E(y2i − (xi β0 + zi η0 )) f i zi = 0
Eh 3i (θ0 ) ≡ Eσi f i zi (σi f i xi − σi f i zi θ0 ) = 0

(11)
(12)
(13)

where the first equation is Eq. (6) above, the second equation is the moment condition
that identifes η0 and the third equation is a modified (through the inclusion of thescalar σi ) version of the OLS first-order conditions used to estimate θ0 .
We note that, due to a separability result of Ahn and Schmidt (1995), we cannot
improve on the estimation of θ0 by estimating it jointly with (β0 , η0 ) because theadditional conditions (11)–(12) only determine θ0 in terms of β0 and η0 . See also
Prokhorov and Schmidt (2009, Statement 6). We further note that by Statement 7 ofProkhorov and Schmidt (2009), joint estimation of the entire vector (β0 , η0 , θ0 ) is
equivalent to a two-step estimation where θ0 is estimated first and the second step isadjusted for the estimation error of the first step.
More importantly, because the correlation between the moment functions in
Eqs. (11) and (12) is zero and the expected derivative of Eq. (11) with respect to η is123

DS-HECK: double-lasso estimation of Heckman selection model

3175

zero, the condition of partial redundancy of Breusch et al. (1999, Theorem 7) holds(in their notation G 21 = 0 and 21 = 0). This means the moment condition (12) isredundant for the estimation of β (M-redundancy). Additionally, these conditions are
sufficient to show that the knowledge of the value of η0 is redundant (P-redundancy).
See Statement 4 of Prokhorov and Schmidt (2009). So the NO condition of Belloni
et al. (2016b) corresponds to a well established situation in GMM estimation whenneither the knowledge of the parameter η0 , nor the knowledge of the moment condition(12) helps estimate β0 .1
2.4 Choice of penalty parameter
The penalty parameters, λ1 and λ2 can be derived analytically but typically, data-driven
methods are used. Their theoretical validity and practical performance have been
well studied. For example, cross-validation or AIC typically under-penalize (overselect) by including too many variables to reduce bias, compared with BIC or plug-in
methods. Additionally, when too many variables are selected, this violates the sparsity
assumption required for the double lasso to work.
Plug-in methods have been shown to perform well in a multitude of settings (see,e.g., Drukker and Liu 2022; Belloni et al. 2014b, 2016a).
2.5 Estimation of the main equation
We can now return to the estimation of α and γ . Similar to Belloni et al. (2016b),
Belloni et al. (2014b) observe that the direct application of lasso to linear models
with a large-dimensional nuisance parameter results in a biased estimation of theparameter of interest, which in their case is a scalar treatment effect. They propose adouble selection procedure. We follow their approach subject to a few modifications
that reflect the specifics of our main equation.
First, with a consistent estimator of β, a natural estimator of the inverse Mills ratio
in Eq. (4) is as follows:





λ(x
i β) = φ(xi β̌)/(xi β̌).

It is also natural to account for the fact that this is a generated regressor when constructing the variance matrix, something we consider later.
Second, Belloni et al. (2014b) derive their results for a scalar parameter. Because thevariables of interest x form a vector, we need to extend the original double selection
lasso estimation to vectors. We provide the details of this extension using the NO
arguments in Appendix A.
1 A qualification to this observation is that z has to have fixed and moderate dimensionality so that GMM
is feasible. We thank a referee for pointing this out to us.

123

3176

M. Hirukawa et al.

We can now discuss the estimation of the main equation which combines the double
selection lasso of Belloni et al. (2014b) and parameter tuning2 by Drukker and Liu(2022). It proceeds in three steps:
Step 1 . We run the lasso regression of y1 on zθ̌ y = arg min En (y1i − zi θ y )2 + λ1 ||θ y ||1 .
θyThis produces a subset of z indexed by suppor t(θ̌ y ).
Step 2 . For each variable x1 j in x1 , j = 1, . . . , k1 , we run the lasso regression of x1 j
on z:


θ̌ j = arg min En (x1i j − zi θ j )2 + λ2 ||θ j ||1 .
θjAdditionally, we run the lasso regression of λ(xi β) on z: β) − z θ )2 + λ ||θ || .
θ̌λ = arg min En (λ(x
λ
2 λ 1
i
i
θλ

This step produces subsets of z indexed by suppor t(θ̌ j ), j = 1, . . . , k1 , andsuppor t(θ̌λ ).
Step 3 . We run the regression of y1i on x1i , λ(xi β), and the union of the sets selectedin Steps 1 and 2:



 2


(
α, 
γ ,
δ ) = arg min En (y1i − x1i
α − λ(xi β)γ − zi δ) ,α,γ ,δ

where suppor t(δ) ⊆ suppor t(θ̌ y ) ∪ suppor t(θ̌1 ) ∪ . . . ∪ suppor t(θ̌k1 ) ∪suppor t(θ̌λ )Proposition 2 Under Assumptions 1–3, the DS-HECK estimation in Steps 1–3 above
is consistent for α and γ and post-double-selection inference on α and γ is valid.
The DS-HECK estimator corrects the bias generated by applying the lasso directly
to Eq. (4). The simulation experiments we report in Sect. 3 illustrate the size of bias.
Following Belloni et al. (2014b), we can claim that inference about the vector(α  , γ ) is valid but, unlike Belloni et al. (2014b), it is valid up to the variance matrix
correction reflecting the post-lasso probit estimation of β.3
2 The notation for tuning parameters λ , λ is a generic notation, the same in this section as in Section
1 2

2.2, but the actual values of the tuning parameters will differ between each instance of their use. Also, thisnotation is different from λ(·), which denotes the inverse Mills ratio.
3 Our implementation of Steps 1–2 follows (Belloni et al. 2014b) and includes penalty loadings which we
do not include in the notation for simplicity. See Belloni et al. (2014b, p. 615–616 and Appendix A).

123

DS-HECK: double-lasso estimation of Heckman selection model

3177

2.6 Variance matrix estimation
We start with some new notation. Let λ̂i = λ(xi β) = λ(xi β̌) and define

ξi = λ̂i (λ̂i + xi β̌),
where β̌ is obtained by the double selection probit. Let e denote the vector of residuals
from the last step of the double selection lasso estimation, with typical element ei , i =1, . . . , n. That is,θ,
ei = y1i − x1i

α −
λi 
γ − zi 

where suppor t(θ ) ⊆ suppor t(θ̌ y )∪suppor t(θ̌1 )∪. . .∪suppor t(θ̌k1 )∪suppor t(θ̌λ ).
Let W denote the matrix containing x1 , the n × 1 vector of λ̂i ’s, and the variables inz that survived the double selection. Let R be a n × n diagonal matrix, with diagonalγ /σ and σ̂ 2 = (e e + 
γ 2 i ξi )/n.
elements (1 − ρ̂ 2 ξi ), where ρ̂ = 
Proposition 3 A consistent estimator of the variance matrix of the DS-HECK estimator
(
α, 
γ,
θ  ) isV = σ̂ 2 (W  W )−1 (W  RW + Q)(W  W )−1 ,where
Q = ρ̂ 2 (W  Dx)Vb (x DW )and Vb is the “sandwich” variance matrix for the double selection probit estimator β̌and D is the diagonal matrix with diagonal elements ξi .
The variance for 
α and γ is the upper (k1 + 1) × (k1 + 1) submatrix of V . The
dsheckman command implements this variance estimator.

3 Monte Carlo simulations
To evaluate the finite-sample performances of DS-HECK, we conduct a simulation
study using four estimators: (i) ordinary least squares on the selected sample (OLS), (ii)Heckman two-step estimator based on the true model (Oracle), (iii) Heckman two-stepestimator using lasso to select variables in Eq. (2) (Naive), and the proposed double
selection Heckman estimator (DS). We have implemented the DS-HECK estimatorin the Stata command dsheckman, available on the authors’ web pages along withdata sets and codes for the simulations and application, and we describe its syntax inAppendix C.
OLS is inconsistent unless there is no sample selection bias, i.e. γ = 0. Naive isinconsistent due to error made by lasso. Moreover, Naive does not provide valid inference as it is not robust to the model selection bias. In contrast, DS is expected to retain

123

3178

M. Hirukawa et al.

consistency in the presence of sample selection biases and show robustness against
the model selection bias. Oracle is expected to behave like the standard Heckman
estimator under the true model but, in practice, Oracle is infeasible since we do notknow the true model.
3.1 Setup
Our data generating process is as follows:
y1 = 1 + x1 + x2 + u 1
y2 = I(x β0 + z η0 + u 2 ≥ 0)
where u 2 ∼ N (0, 1) and where u 1 = γ u 2 + and ∼ N (0, 1) is independent of u 2 .
We vary the strength of the selection bias by setting γ to be 0.4, 0.5, 0.6, 0.7, and 0.8
and we observe y1 only when y2 = 1.
The selection equation is generated using nine non-zero variables in z of whichfour have a relative large effect and five relatively small:x β0 + z η0 = −1.5 + x1 − x2 + z 1 − z 2 + 0.046z 3 + z 5− 0.046z 10 − 0.046z 11 + 0.046z 12 − 0.046z 15 + z 20 .
The value 0.046 is chosen so that it violates the so called “beta-min” condition andcauses lasso to make model selection mistakes (see, e.g., Liu et al. 2020; Drukker andLiu 2022). The sample size is 2000. The number of replications is 1000.
We consider two scenarios for p, the dimension of z: (i) p = 1000, fewer variables
than observations; (ii) p = 2100, more variables than observations. The variables
are generated using a Toeplitz correlation structure with decreasing dependence. In
particular, let Z be the matrix of dimension N × p containing z. Then,Z = M L√
where M is N × p and has the typical element (ζi j −15)/ 30, where ζi j ∼ χ 2 (15) and
L is the Cholesky decomposition of a symmetric Toeplitz matrix V of dimension p× p
such that its elements obey the following laws: Vi, j = Vi−1, j−1 and V1, j = j −1.3 .
The variables x are also correlated and they are generated as functions of z. In
particular,
x1 = z 3 + z 10 + z 11 + z 12 + z 15 + x1
x2 = 0.5(z 3 + z 10 + z 11 + z 12 + 2z 15 ) + x2
where x1 and x2 follow a Toeplitz structure similar to Z .
As a result, for the selected sample, the true model is
y1 = 1 + x1 + x2 + γ λ(x β0 + z η0 ) + u

123

DS-HECK: double-lasso estimation of Heckman selection model

3179

where λ(·) = φ(·)/(·) and the true parameter values are β1 = β2 = 1 and γ =0.4, 0.5, 0.6, 0.7, or 0.8. The standard deviation of λ(1) (·) is about 0.3.4
3.2 Results
For each estimator, we report the following measures: (i) true value of parameter (True),(ii) mean squared error (MSE), (iii) average of estimates across simulations (Mean)(iv) standard deviation of estimates across simulations (SD), (v) average standard error
across simulations (S E) and (vi) rejection rate for the H0 that the parameter equals itstrue value against the nominal 5% level of significance (Rej. rate).
We report the simulation results for β1 , β2 , and γ in Tables 1, 2 and 3, respectively.
Several observations are clear from the tables. First, OLS is badly biased and the bias
is greater when selection is stronger. Second, Naive is also biased and it fails to provide
valid inference at any value of γ and p. This demonstrates that Naive is not robust to
the model selection errors. The rejection rate and MSE increase as γ increases, whichis expected because greater γ value indicate a greater sample selection bias. Third,
Oracle shows consistency and rejection rates close to the nominal 5% significancelevel. Fourth, DS performs similarly to Oracle for all values of γ and p. In particular,
its MSE is consistently smaller than for Naive and OLS, S E is close to SD, whichshows that the proposed variance adjustment works well. Finally, Rej. rate for DS isnear the 5% significance level, which supports that our estimator offers valid inference.

4 Application to female earnings estimation
4.1 Labor force participation and earnings
Estimation of the female earnings equation is a topic of long standing interest among
economists. Early investigations of the labor supply decisions including both participation and hours by married women date back to Gronau (1974) and Heckman(1974), who were among the first to highlight sample selection bias stemming fromthe labor supply decision. Labor market decisions by women over the later decades
have been studied and documented by Mroz (1987), Ahn and Powell (1993), Neumark
and Korenman (1994), Vella (1998), Devereux (2004), Blau and Kahn (2007), Mulligan and Rubinstein (2008), Cavalcanti and Tavares (2008), Bar et al. (2015), among
others.
Numerous applied works have extensively scrutinized the empirical aspects of thesample selection problem when estimating female labor market behavior. The determinants of the earnings equations for married women are similar to those of men
and have been mostly agreed upon. These determinants traditionally include women’s
education, experience, age, tenure and location. The hallmark of correcting for sample selection bias is finding some appropriate exclusion restriction(s) (i.e., variable(s)affecting selection but not the earnings) in order to ensure proper identification of
4 We ran simulations with higher variability of λ(1) (·). This did not affect the performance of dsheckman

for β1 and β2 . For γ , when p = 2100, we noticed slight over-rejection with no effect on key takeaways.

123

3180

M. Hirukawa et al.

Table 1 Simulation results for β1
p = 1000
True MSE

p = 2100
Mean

SD

SE

Rej. Rate MSE

Mean

SD

SE

Rej. Rate

γ = 0.4
Oracle10.0014 0.9984 0.0369 0.0374 0.0450

0.0013 0.9985 0.0362 0.0374 0.0510

DS1

0.0025 0.9977 0.0501 0.0484 0.0540

0.0022 0.9991 0.0464 0.0486 0.0380

Naive

1

0.0049 0.9424 0.0395 0.0350 0.4200

0.0052 0.9379 0.0360 0.0349 0.4310

OLS

1

0.0046 0.9422 0.0350 0.0350 0.3890

0.0045 0.9421 0.0339 0.0350 0.3610

Oracle10.0014 0.9993 0.0380 0.0387 0.0430

0.0015 0.9980 0.0389 0.0387 0.0480

DS1

0.0025 1.0000 0.0504 0.0502 0.0570

0.0027 0.9979 0.0518 0.0502 0.0600

Naive

1

0.0067 0.9289 0.0409 0.0362 0.5240

0.0075 0.9229 0.0389 0.0361 0.5840

OLS

1

0.0063 0.9291 0.0353 0.0362 0.4930

0.0066 0.9277 0.0373 0.0363 0.5190

γ = 0.5

γ = 0.6
Oracle10.0016 0.9995 0.0395 0.0401 0.0440

0.0015 0.9984 0.0393 0.0400 0.0520

DS1

0.0025 1.0001 0.0504 0.0519 0.0460

0.0027 0.9981 0.0522 0.0521 0.0500

Naive

1

0.0089 0.9183 0.0475 0.0376 0.6110

0.0101 0.9081 0.0413 0.0373 0.6880

OLS

1

0.0087 0.9148 0.0374 0.0376 0.6190

0.0087 0.9144 0.0370 0.0376 0.6150

γ = 0.7
Oracle10.0019 0.9969 0.0435 0.0416 0.0690

0.0018 1.0001 0.0419 0.0419 0.0510

DS1

0.0032 0.9956 0.0564 0.0539 0.0570

0.0028 1.0008 0.0525 0.0543 0.0410

Naive

1

0.0126 0.8999 0.0508 0.0389 0.7400

0.0135 0.8930 0.0455 0.0390 0.7640

OLS

1

0.0118 0.8993 0.0401 0.0392 0.7170

0.0115 0.9005 0.0397 0.0393 0.7110

γ = 0.8
Oracle10.0019 0.9962 0.0437 0.0434 0.0470

0.0019 0.9968 0.0438 0.0434 0.0510

DS1

0.0031 0.9965 0.0560 0.0558 0.0470

0.0034 0.9966 0.0582 0.0562 0.0530

Naive

1

0.0165 0.8835 0.0544 0.0406 0.8020

0.0173 0.8777 0.0484 0.0404 0.8360

OLS

1

0.0152 0.8838 0.0408 0.0409 0.8150

0.0148 0.8858 0.0414 0.0408 0.7870

the model parameters. Two main competing choices of such exclusion restrictions
have been exploited for estimation of labor market earnings for married women:
non-wife/husband’s income and the existence or number of (young) children. The
underlying argument is that these two sets of variables affect the labor supply decision of married women but not their earnings. Huber and Mellace (2014) provide
an overview of the related literature on sample selection bias in the female earnings
equations.
Cavalcanti and Tavares (2008) provide an alternative view and argue that the declining price and wider availability of home appliances play a crucial role in explaining therise in female labor force participation. This suggests a long list of potential exclusion
restrictions. Furthermore, the exact nature and functional form of the chosen exclusion
restriction(s) in the selection equation is uncertain. For example, should labor work
experience include years of part-time employment? Should educational attainment be123

DS-HECK: double-lasso estimation of Heckman selection model

3181

Table 2 Simulation results for β2
p = 1000
True MSE

p = 2100
Mean

SD

SE

Rej. Rate MSE

Mean

SD

SE

Rej. Rate

γ = 0.4
Oracle10.0028 1.0018 0.0527 0.0524 0.0530

0.0027 1.0014 0.0518 0.0524 0.0500

DS1

0.0033 1.0015 0.0576 0.0584 0.0390

0.0034 1.0025 0.0580 0.0587 0.0530

Naive

1

0.0145 1.1018 0.0641 0.0530 0.5210

0.0151 1.1093 0.0564 0.0532 0.5340

OLS

1

0.0057 1.0541 0.0523 0.0512 0.1820

0.0054 1.0538 0.0506 0.0512 0.1760

Oracle10.0029 0.9998 0.0536 0.0542 0.0470

0.0031 1.0029 0.0558 0.0542 0.0620

DS1

0.0036 1.0005 0.0600 0.0603 0.0380

0.0038 1.0030 0.0613 0.0606 0.0600

Naive

1

0.0200 1.1249 0.0659 0.0547 0.6420

0.0224 1.1365 0.0614 0.0548 0.7020

OLS

1

0.0070 1.0651 0.0523 0.0530 0.2460

0.0077 1.0683 0.0546 0.0530 0.2540

γ = 0.5

γ = 0.6
Oracle10.0031 1.0014 0.0555 0.0562 0.0510

0.0030 1.0020 0.0548 0.0560 0.0450

DS1

0.0038 1.0034 0.0620 0.0623 0.0450

0.0038 1.0038 0.0615 0.0628 0.0410

Naive

1

0.0271 1.1453 0.0775 0.0566 0.7150

0.0307 1.1626 0.0653 0.0567 0.8000

OLS

1

0.0095 1.0804 0.0550 0.0550 0.3100

0.0092 1.0795 0.0536 0.0549 0.2950

γ = 0.7
Oracle10.0036 1.0031 0.0596 0.0583 0.0620

0.0034 0.9997 0.0582 0.0585 0.0430

DS1

0.0042 1.0034 0.0649 0.0648 0.0590

0.0043 0.9994 0.0653 0.0657 0.0430

Naive

1

0.0377 1.1755 0.0833 0.0587 0.8170

0.0414 1.1902 0.0725 0.0591 0.8680

OLS

1

0.0121 1.0934 0.0583 0.0573 0.3740

0.0118 1.0925 0.0572 0.0575 0.3640

γ = 0.8
Oracle10.0037 1.0060 0.0606 0.0607 0.0440

0.0037 1.0016 0.0611 0.0608 0.0520

DS1

0.0049 1.0064 0.0698 0.0676 0.0630

0.0049 1.0010 0.0703 0.0683 0.0530

Naive

1

0.0505 1.2059 0.0898 0.0611 0.8690

0.0519 1.2143 0.0777 0.0613 0.8990

OLS

1

0.0158 1.1106 0.0598 0.0598 0.4510

0.0147 1.1052 0.0600 0.0598 0.4250

measured in full years of completed education or in milestones such as high school orcollege, as a replacement or complement to years of education? Similarly, should age
enter the model in a linear or quadratic form?
Our goal in this section is to illustrate the performance of the DS-HECK procedureon the following earnings equation:log (ear nings) = α0 + α1 education + α2 x1 + α3 state dummies + u 1 , (14)where log (ear nings) is the natural logarithm of the individual’s total annual labor
income, education is the person’s completed years of education, state dummies is
a vector containing a full set of state dummies, and u 1 is an idiosyncratic error. The
vector x1 varies across the exact specification we consider and can contain age and/or
work experience.

123

123
SD

SE

Rej. Rate

0.4

0.4

0.4

0.5

0.5

0.5

0.6

0.6

0.6

0.7

0.7

0.7

0.8

0.8

0.8

OracleDS

Naive

OracleDS

Naive

OracleDS

Naive

OracleDS

Naive

OracleDS

Naive0.02000.01310.01170.01720.01240.01160.01350.00980.00920.01350.01070.01020.01150.00960.00920.72600.7762

0.7912

0.63590.67650.69030.5454

0.58350.59520.45580.48430.49400.36580.38730.39480.12040.11190.10780.11450.10910.10720.1024

0.0977

0.09580.10770.10220.10110.1017

0.09700.09590.11240.10650.1064

0.10800.10290.10300.1041

0.09980.10000.10100.09690.0973

0.0982

0.09430.09480.13800.07200.05100.11700.07700.06000.07800.04200.03900.08200.06500.05400.07100.05800.04900.02020.01300.01170.01530.01110.01040.01420.01020.00980.01250.00950.00930.01080.00910.0087

MSE

Mean

True

MSE

p = 2100

p = 1000

Table 3 Simulation results for γ

0.72130.76570.78310.64530.68300.69990.54500.58100.59450.4537

0.48370.49480.36460.3877

0.3969Mean

0.11820.10880.10710.11110.10400.10200.10580.09950.09880.1017

0.0963

0.0961

0.09790.0946

0.0933SD

0.11320.10620.1066

0.10910.10250.10310.10500.09930.10000.1013

0.0964

0.0972

0.09880.0942

0.0952SE

0.12400.08200.05600.07400.05100.04100.08000.05000.04600.07600.05500.04600.05900.05100.0430Rej. Rate

3182
M. Hirukawa et al.

DS-HECK: double-lasso estimation of Heckman selection model

3183

To address the potential self-selection bias we employ DS-HECK as a data-driven
procedure for choosing the explanatory variables and functional form (among highdimensional options provided) in the following labor force participation equation:
inl f = I(x β + z η + u 2 ≥ 0),(15)

where inl f is the dummy variable that is equal to one for those women who are in thelabor force at the time of the interview and zero otherwise, and u 2 is an idiosyncraticerror. The vector x includes all the explanatory variables from Eq. (14) (both in a linear
and quadratic functional form) as well as exclusion restrictions.
In practice, x is constructed as follows. To simplify notation, denote all the
explanatory variables from Eq. (14) as x1 . First, running lasso probit of inl f on
high-dimensional controls w, where w includes x1 and some other high-dimensional
controls. Denote the selected variables as x2 . Second, x is union between all x1 andx2 . All the non-selected controls in w are used as z.
4.2 Sample construction
We obtain our sample from the 2013 wave of the Panel Study of Income Dynamics(PSID) where we focus on the sub-population of white married women. The choiceof explanatory variables reflects their availability in the PSID and the traditional setof regressors used in the existing literature on female labor force participation andearnings. Specifically, the explanatory variables we collect from the PSID include
information on the educational attainment of the individual (both as the number ofyears completed and as a set of indicators for milestone achievements in education),a set of indicators for whether the individual obtained her education in the USA,outside the USA, or both, as well as a set of indicators for the educational levels of theindividual’s parents, work experience of the individual, age and geographical locationof the individual (captured by a set of dummy variables for the current state where
the individual is located), a set of indicators reflecting the Beale-Ross rural–urban
continuum code for the individual’s current residence, and an indicator for whether
the individual owns a vehicle. Table 4 contains a description of the key explanatory
variables.
The set of (potential) exclusion restrictions includes the number of children in thehousehold under 18 years of age, an indicator for whether there are any children age 15
years old or younger in the individual’s household, annual labor income of the husband,child care expenses, and household major expenditure (i.e. expenditure on household
furnishings and equipment, including household textiles, furniture, floor coverings,
major appliances, small appliances and miscellaneous housewares). While admittedly
far from ideal, this last variable is the closest information we find in the PSID tocapture household expenditure on major household appliances, which allows us totest the argument of Cavalcanti and Tavares (2008). Finally, the dependent variables
for the earnings and selection equations are (the natural logarithm of) the individual’s
total annual labor income and the indicator for whether the individual is in the labor
force, respectively.

123

3184

M. Hirukawa et al.

Table 4 Description of key PSID variables
VariableType

Definition

Education

ContinuousYears of education

High school education

Categorical

Graduated from high school (3 categories)

Enrolled in school

Categorical

If currently enrolled in regular school (2categories)

College attendanceCategorical

If attended college (2 categories)

Other degree or certificate

Categorical

If received other degree/certificate (2categories)

US education

Categorical

If the individual obtained her education in
the USA, outside the USA, or both (3categories)

Father’s education

Categorical

Educational level of the individual’s father (8categories)

Mother’s education

Categorical

Educational level of the individual’s mother(8 categories)

If owns a vehicle

Categorical

If the individual owns a vehicle (2categories)

Current state

Categorical

Geographical location of the individual (46states)Rural–urban location

Categorical

Beale-Ross rural–urban continuum code forthe individual’s current residence (9categories)

Number of kids

Continuous

The number of children in the household
under 18 years of age

If kids less than or equal to 15 years old

Categorical

If there are any children age 15 years old or
younger in the individual’s household (2categories)

Child care expenditureContinuousChild care expenses (in thousand dollars)

Husband’s labor income

Continuous

Annual labor income of the individual’shusband (in thousand dollas)

Household major expenditureContinuous

Expenses on household furnishings andequipment, including household textiles,
furniture, floor coverings, major
appliances, small appliances andmiscellaneous housewares (in thousand
dollars)

Controls

Potential exclusion restrictions

Our sample contains 1,989 white married women, of whom 1,294 are in the labor
force and 695 are not. Table 5 reports summary statistics for key variables in the dataset.
A set of dummy variables for the current state as well as a set of indicators reflecting
the Beale-Ross rural–urban continuum code are omitted to save space. A total of 46
states are present in our sample, with Delaware, District of Columbia, Hawaii, New
Mexico, and Rhode Island omitted from our sample during data cleaning. We notethat some women report being neither employed nor (temporarily) unemployed while123

DS-HECK: double-lasso estimation of Heckman selection model

3185

also reporting non-zero labor income during that time. There are 161 such women in
the sample. We treat these individuals as not being in the labor force.
4.3 Empirical findings
We consider several specifications when estimating the earnings equation subject tosample selection bias. Table 6 reports the key results for both equations obtained using
DS-HECK. The top panel provides coefficient estimates (as well as their standard
errors) for α1 , α2 and the coefficient on the inverse Mills ratio while the bottom panel
provides estimates (and standard errors) for β. In addition to the reported estimates,
each specification contains two more sets of of estimates which we do not report in thetable to save space. First, we do not report the coefficients on the full set of state andurban-rural dummies included in both equations. Second, for each specification there
are the selected controls in both equations; the number of such controls is reported at
the bottom of the respective panels but the coefficients themselves are not reported.
We note that following the original Heckman specification, the explanatory variables present in the earnings equation (x1 ) are always kept in the labor force
participation equation. Columns (1) and (2) report the estimates when work experience enters the two equations, with and without the quadratic form, while age is notincluded. Columns (3) and (4) report the estimate when age is included but experience
is not. Columns (5) and (6) contain both but differ in whether experience squared isincluded. Child care expenditure is always selected in the post-lasso probit step andwe view it as an important exclusion restriction (part of x not in x1 ), so it is reported
separately from other controls (z). The variables If kids under 15, Number of kids,
Husband’s income, Major expenses are not selected in the post-lasso probit step but
are sometimes selected in Step 2 of the selection model estimation, in which case they
are counted under selected controls.
Column (8) reports the estimates of the traditional Heckman specification with
Experience, Experience2 , Age and Age2 in both equations and Child care expenditure,
If kids under 15, Number of kids, Husband’s income, Major expenses in the participation equation. No selection is used on additional controls. In column (7) we report a
specification that forces exactly the same variables as in the classical Heckman specification in column (8) to be included and additionally, permits the lasso to select any
additional controls.
As Table 6 suggests, the signs of all the reported coefficient estimates are as
expected, and they are highly significant for the most part. We note that from thefive potential exclusion restrictions of interest, the lasso selects child care expenditure
as the relevant covariate for the labor force participation equation. Finally, we notethat the results for the labor force participation equation reported in column (7) and(8) are similar to those reported in columns (5) and (6) except that (7) and (8) use theadditional exclusion restrictions some of which turn out to be statistically significant
while (5) and (6) drop them.

123

123
1.04

Received education in the USA (1–3 range)

0.47
1.47
61.49
1.49

If kids less than or equal to 15 years oldChild care expenditure

Husband’s labor income

Household major expenditure
45,575.65
10.37
1
1,294

Labor income

Log (Labor income)

If in the labor force

Number of observations

Outcomes

1.00

Number of kids under 18 years old

Exclusion restrictions

0.99

0.21

If received other degree/certificate

If owns a vehicle

0.82

If attended college

5.08

0.03

If currently enrolled in regular school

4.96

1.05

Graduated from high school (1–3 range)

Mother’s educational level (1–8 range)

14.78

Education (in years)

Father’s educational level (1–8 range)

43.84
9.86

Experience

In the labor force
Mean

Age

Controls

Table 5 Key sample characteristics

0

0.98

43,701.59

2.90

81.16

3.91

0.50

1.17

0.09

1.61

1.80

0.30

0.41

0.39

0.17

0.27

1.98

6.64

11.74

SD

695

0

8.52

12,439.52

1.20

62.76

0.25

0.35

0.82

0.98

4.54

4.56

1.04

0.20

0.68

0.03

1.14

13.96

12.21

53.63

0

1.58

15,746.09

2.69

195.65

1.64

0.48

1.27

0.15

1. 64

1.92

0.27

0.40

0.47

0.16

0.46

2.16

9.30

16.17

Not in the labor force
Mean
SD

1,989

0.65

10.17

41,909.04

1.39

61.93

1.04

0.43

0.93

0.99

4.81

4.90

1.04

0.20

0.77

0.03

1.09

14.49

10.68

47.26

All
Mean

0.48

1.21

42,822.50

2.83

132.85

3.35

0.49

1.21

0.12

1.64

1.86

0.29

0.40

0.42

0.17

0.35

2.08

7.75

14.24

SD

3186
M. Hirukawa et al.

19

Experience2

Experience

Education (years)− 0.001***(0.000)− 0.001***(0.000)0.012(0.013)0.012(0.013)0.066***(0.024)0.066***

19(0.024)

Labor force participation equation

Number of selected controls

34

1294

(0.484)1294− 1.463***(0.388)35

0.018***(0.025)0.074***

20

34

1294

(0.146)− 0.838***(0.025)0.074***

21

34

1294

(0.456)− 1.584***(0.015)0.042***(0.025)0.070***

21

35

1294

(0.171)(0.000)− 0.001**(0.015)0.041***(0.025)0.069***

21

34

1294

(0.455)− 1.581***(0.001)(0.055)0.001*− 0.924***(0.005)(0.001)(0.057)(0.004)− 0.0860.001*− 0.0900.019***(0.001)− 1.336***

Number of controls− 0.018(0.021)− 0.000(0.015)(0.033)0.060*(6)(0.023)0.089***(5)(0.001)(0.034)0.056*(4)0.001(0.018)(0.007)(0.023)0.089***

(3)

0.001

0.006

0.026***

0.047(0.040)

0.053

(2)

(0.034)(1)

Number of observations (n)

λ

Age2

Age

Experience2

Experience

Education (years)

Earnings equation

Table 6 Estimation results(0.000)− 0.001**(0.015)0.042***(0.025)0.068***

18

29

1294

(0.167)− 0.783***(0.000)0.000(0.027)− 0.001(0.001)0.001(0.017)− 0.012(0.026)0.078***(7)(0.000)− 0.001**(0.015)0.042***(0.026)

0.067***

1294(0.156)− 0.942***(0.000)0.000(0.022)0.002(0.000)0.000(0.014)0.001(0.016)0.075***(8)

DS-HECK: double-lasso estimation of Heckman selection model
3187

123

123

1989

35

21

Number of controls

Number of selected controls

21

34

1989

0.073***(0.025)0.073***(2)(0.025)(1)

Number of observations (N)Major expenses

Husband’s income

Number of kids

If kids under 15

Child care expenditure

Age2

Age

Table 6 continued

22

34

1989

(0.026)

22

34

1989

(0.026)0.075***(0.000)0.075***− 0.002***(0.000)− 0.002***0.168***(0.023)0.168***(4)(0.023)

(3)

23

35

1989

(0.026)0.072***(0.000)− 0.002***(0.022)0.160***

(5)

22

34

1989

(0.025)0.072***(0.000)− 0.002***(0.023)0.159***(6)(0.000)− 0.001(0.000)

20

29

1989

1989

0.014

− 0.001**(0.046)(0.012)(0.047)− 0.199***(0.020)− 0.195***(0.129)0.007− 0.288**(0.130)− 0.248*(0.016)0.085***(0.000)− 0.002***(0.020)0.164***(8)(0.025)0.067***(0.000)− 0.002***(0.023)0.152***(7)3188
M. Hirukawa et al.

DS-HECK: double-lasso estimation of Heckman selection model

3189

Next we focus on the estimates of the labor income equation. According to theresults reported in Table 6, we conclude that the educational level of the individual
plays a crucial role in explaining labor income for white married women in 2012.
When statistically significant, the estimated rate of return to education ranges from
5.6% to almost 9% depending on the specification. Furthermore, there is evidence thatthe individual’s age is more important, both statistically and practically, than work
experience for explaining the individual’s labor income in our sample. We note thatwhen the individual’s age (in any functional form) is used in the labor income equation,
the rate of return to education is statistically significant.
Most importantly, Table 6 suggests that the inverse Mills ratio is highly statisticallysignificant in all specifications implying that the correction for the sample selection
was needed. Given the economic interpretation of the estimated coefficients, theirsigns and economic as well as statistical significance, specification (5) seems most
attractive in light of the existing studies on the topic. Interestingly, the traditional
Heckman specification produces results that are close to those reported in column (5).

5 Conclusion
We have proposed an extension to the traditional Heckman sample selection modelby incorporating it into a model selection framework with many controls. A double
application of the lasso to each part of the model permits valid inference on the parts ofthe model that is of interest to empirical economists. We detail the steps involved in aconsistent estimation with valid standard errors and we provide a new Stata command
to implement it. We also investigate aspects of the Neyman orthogonality that relate
it to the concept of redundancy in the GMM estimation.
Lasso and double selection in linear models have been recently subject of scrutiny
in cases when lasso under-selects controls in finite samples even under a sparsity
scenario and the double selection estimators have severe omitted variable biases (see,e.g., Wuthrich and Zhu 2021; Lahiri 2021). This happens when the signal from thevariables lasso works on is weak and they do not get selected by either of the two
selection procedures. The solution proposed by Wuthrich and Zhu (2021) is to resort
in such cases to the regular OLS estimation using a high-dimensional variance matrix
computations which is computationally difficult and works only when p < n.
We show how large the errors committed by the naive application of lasso can be
and we provide an application to a classic problem in labor economics where using
our method leads to a few new insights. We provide a user-friendly and versatile Stata
command, which can help empirical economists use the proposed methodology. The
command as well as the simulation and application data are made available on theauthors’ web pages.

123

3190

M. Hirukawa et al.

Finally we note that the results of this paper can be extended to other consistentmethods of model selection beyond lasso, such as the Dantzig Selector proposed by
Candes and Tao (2007).
Funding Open Access funding enabled and organized by CAUL and its Member Institutions Masayuki
Hirukawa gratefully acknowledges financial support from Japan Society for the Promotion of Science(grant number 19K01595).

Declarations
Conflicts of interest Masayuki Hirukawa declares that he has no conflict of interest. Di Liu declares that hehas no conflict of interest. Irina Murtazashvili declares that she has no conflict of interest. Artem Prokhorov
declares that he has no conflict of interest.
Ethical approval This article does not contain any studies with human participants or animals performedby any of the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence,and indicate if changes were made. The images or other third party material in this article are includedin the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. Ifmaterial is not included in the article’s Creative Commons licence and your intended use is not permitted
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix
A Neyman orthogonal estimation for vectors
Consider an extension of the partial linear framework of Belloni et al. (2014a):
Y = D  θ0 + g0 (X ) + U

E(U |D, X ) = 0

D1 = m 0,1 (X ) + V1..
.

E(V1 |X ) = 0

Dk = m 0,k (X ) + Vk

E(Vk |X ) = 0

where D and θ0 are k×1-vectors, rather then scalars. The functions g0 (X ) and m 0, j (X )are unknown but we can use approximately sparse linear models to approximate thesefunctions so that the approximation errors is small enough. We know that if θ0 is
a scalar, we can apply the double selection lasso to achieve valid inference for θ0 .
Here, we will show that the same arguments can be applied when when θ0 is a lowdimensional vector.
We follow Chernozhukov et al. (2018) and show that
1. the moment condition implied by the Robinson-style “partialling-out” is Neyman
orthogonal.
2. the application of Theorem 3.1 of Chernozhukov et al. (2018) provides the asymptotic distribution (note that cross-fitting only permits to relax sparsity requirement,

123

DS-HECK: double-lasso estimation of Heckman selection model

3191

and the “partialling-out” approach still provides valid inference for θ even withoutcross-fitting).
3. the “partialling-out’ approach is equivalent to the double selection lasso.
We start with a general definition of Neyman orthogonal moments (see, e.g., Chernozhukov et al. 2018, Definition 2.1). Suppose we wish to estimate a low-dimensionalparameter θ0 using the moment conditions
E [ψ(W ; θ0 , η0 )] = 0,(16)where W contains the random variables and η0 is a nuisance parameter vector, which
can be high-dimensional. Suppose we use machine learning techniques to estimate
η0 . In essence, Neyman orthogonality means that small mistakes in the estimation ofη will not disturb the consistent estimation of θ0 .
Formally, Neyman orthogonality depends on the concept of pathwise (Gateux)
derivative. Let Dr denote the Gateux derivative of the moment condition ψ withrespect to η in direction r . Then,Dr [η − η0 ] = ∂r {E [ψ(W ; θ0 , η0 + (η − η0 )r )]}(17)

for all r ∈ [0, 1). When Dr is evaluated at r = 0, we denote it as D0 [η − η0 ].
Definition 1 The moment condition E [ψ(W ; θ0 , η0 )] = 0 is Neyman orthogonal if
D0 [η − η0 ] = 0.
Now we show Neyman orthogonality of the “partialing-out” approach. For notational simplicity, we group the equations for D j ( j = 1, . . . , k).
Y = D  θ0 + g0 (X ) + U ,E(U |D, X ) = 0,

(18)

E(V |X ) = 0,

(19)

D = m 0 (X ) + V ,where D = (D1 , D2 , . . . , Dk ) , m 0 (X ) = (m 0,1 (X ), m 0,2 (X ), . . . , m 0,k (X )) , andV = (V1 , V2 , . . . , Vk ) . Plugging Eq. (19) into (18), we obtain the reduced form forY:
Y = l0 (X ) + B,

E(B|X ) = 0,(20)

where l0 (X ) = m 0 (X )θ0 + g0 (X ) and B = U + V  θ0 . Now we can show thatthe following (Robinson 1988)-style “partialling-out” moment condition is Neyman
orthogonal.
Proposition 4 The moment condition using the function
ψ [W ; θ, η] = (Y − l(X ) − (D − m(X )) θ )(D − m(X ))(21)

where W = (Y , X , D) and η = (l, m), is Neyman orthogonal.

123

3192

M. Hirukawa et al.

Proof First we show that E [ψ (W ; θ0 , η0 )] = 0:E [ψ (W ; θ0 , η0 )] = E (Y − l0 (X ) − (D − m 0 (X )) θ0 )(D − m 0 (X ))= E (B − V  θ0 )V
= E (U + V  θ0 − V  θ0 )V= E [U V ]= 0,
where the last equality holds because E(U |D, X ) = 0
Next, we prove that D0 [η − η0 ] = 0. First note thatE [ψ(W ; θ0 , η0 + (η − η0 )r )]= E [(Y − l0 (X ) − r (l(X ) − l0 (X ))−(D − m 0 (X ) − r (m(X ) − m 0 (X ))) θ )(D − m 0 (X ) − r (m(X ) − m 0 (X ))Thus, we can compute Dr [η − η0 ] as follows:Dr [η − η0 ] = − E [(Y − l0 (X ) − r (l(X ) − l0 (X ))−(D − m 0 (X ) − r (m(X ) − m 0 (X ))) θ0 )(m(X ) − m 0 (X ))− E (l(X ) − l0 (X )) − (m(X ) − m 0 (X )) θ0 )(D − m 0 (X ) − r (m(X ) − m 0 (X ))Now, set r = 0 and evaluate D0 [η − η0 ]:
D0 [η − η0 ] = − E Y − l0 (X ) − (D − m 0 (X )) θ0 (m(X ) − m 0 (X ))− E (l(X ) − l0 (X )) − (m(X ) − m 0 (X )) θ0 ) (D − m 0 (X ))= − E [U (m(X ) − m 0 (X ))] − E (l(X ) − l0 (X )) − (m(X ) − m 0 (X )) θ0 )V= 0,

where the last equality holds because of the law of iterated expectation.
Therefore, the moment condition is Neyman orthogonal.


ψi (W ; θ, η)
Next, we provide the asymptotic distribution. Let θ be a solution to n1
= 0, given η. By Theorem 3.1 of Chernozhukov et al. (2018),
θisaconsistentestimator√of θ0 , and it is asymptotically normal with the rate of N . Its variance is = J0−1 E(ψ(W ; θ0 , η0 )ψ(W ; θ0 , η0 ) )(J0−1 ) ,(22)

where J0 = E[∂ψ/∂θ |θ=θ0 ]. We can simplify this expression to arrive at the following
formula: = (E(V V  ))−1 E(V V  U 2 )(E(V V  ))−1 .

123

(23)

DS-HECK: double-lasso estimation of Heckman selection model

3193

Equation (23) leads to the following straightforward estimator of :=

1   
Vi Vi
n
i

−1 

1     2
Vi Vi Uini

1   
Vi Vi
n

−1
,

(24)

i

i and U
i are the residuals from Eqs. (18) and (19), respectively. Apparently,
where Vthe variance estimator is the classic heteroskedasticity-consistent estimator.
The last step is to show that double selection lasso is equivalent to the “partiallingout”. We start with the equation
y = D  θ0 + g0 (X ) + U .
If we partial out X from both Y and D, the equation becomesy − l0 (X ) = (D − m 0 (X )) θ0 + U(25)

It is easy to see that the solution for θ in Eq. (25) comes from the Neyman orthogonal
moment condition using the moment function in Eq. (21).
Therefore, double selection is equivalent to the Robinson-style “partialling-out”
approach.

B Sketch of Proofs of Propositions 2 and 3
Proof of Proposition 2. The proof is similar to that of Theorem 1 of Belloni et al.
(2014b) and essentially consists in verifying that the assumptions of that theorem are
fulfilled given what we assumed in Assumptions 1–3.

Proof of Proposition 3. The construction of a variance estimator outlined in the paragraphs immediately preceding the proposition is similar to the standard correction for
the first-step estimation error in the Heckman (1979) procedure (see, e.g., Cameron
and Trivedi 2005, Section 16.10.2). An important difference is that now we include theeffect of the controls selected in the last step of the double selection lasso estimation
when computing the residuals.
C dsheckman: Stata command for DS-HECK estimator
The syntax of dsheckman isdsheckman depvar indepvars [if] [in], selection(depvar_s = indepvars_s)[selvars(varlist)]
where• depvar specifies the dependent variable in the main equation, which corresponds
to y1 in Eq. (1).

123

3194

M. Hirukawa et al.

• indepvars specifies the independent variables in the main equation, which correspond to x1 in Eq. (1).
• depvar_s specifies the dependent variable in the selection equation, which corresponds to y2 in Eq. (2).
• indepvars_s specifies the independent variables in the selection equation, whichcorrespond to x and z in Eq. (2).
• selvars() specifies x in Eq. (2). If this option is not used, x is constructed in
two steps: (a) run the lasso probit of devpar_s on indepvars_s using the plugin
penalty; denote the selected variables by x2 ; (b) construct x as the union of x1 andx2 .

References
Ahn H, Powell JL (1993) Semiparametric estimation of censored selection models with a nonparametricselection mechanism. J Econom 58:3–29
Ahn SC, Schmidt P (1995) A separability result for GMM estimation, with applications to GLS prediction
and conditional moment tests. Econom Rev 14:19–34
Arellano M, Bonhomme S (2017) Quantile selection models with an application to understanding changes
in wage inequality. Econometrica 85:1–28
Bar M, Kim S, Leukhina O (2015) Gender wage gap accounting: the role of selection bias. Demography
52:1729–1750
Belloni A, Chernozhukov V, Hansen C (2014a) High-dimensional methods and inference on structural andtreatment effects. J Econ Perspect 28:29–50
Belloni A, Chernozhukov V, Hansen C (2014b) Inference on treatment effects after selection among highdimensional controls. Rev Econ Stud 81:608–650
Belloni A, Chernozhukov V, Hansen C, Kozbur D (2016a) Inference in high-dimensional panel models with
an application to gun control. J Bus Econ Stat 34:590–605
Belloni A, Chernozhukov V, Wei Y (2016b) Post-selection inference for generalized linear models withmany controls. J Bus Econ Stat 34:606–619
Blau FD, Kahn LM (2007) Changes in the labor supply behavior of married women: 1980–2000. J Law
Econ 25:393–438
Breusch T, Qian H, Schmidt P, Wyhowski D (1999) Redundancy of moment conditions. J Econom 91:89–
111
Cameron AC, Trivedi PK (2005) Microeconometrics: methods and applications. Cambridge University
Press
Candes E, Tao T (2007) The Dantzig selector: statistical estimation when p is much larger than n. Ann Stat
35(6):2313–2351
Cavalcanti TVDV, Tavares J (2008) Assessing the “Engines of liberation”: home appliances and female
labor force participation. Rev Econ Stat 90:81–88
Chernozhukov V, Chetverikov D, Demirer M, Duflo E, Hansen C, Newey W, Robins J (2018) Double/debiased machine learning for treatment and structural parameters. Econom J 21:C1–C68
Devereux PJ (2004) Changes in relative wages and family labor supply. J Hum Resour 39:696–722
Drukker D, Liu D (2022) Finite-sample results for lasso and stepwise Neyman-orthogonal Poisson estimators. Econ Rev (Forthcoming)
Gronau R (1974) Wage comparisons—a selectivity bias. J Polit Econ 82:1119–1143
Heckman J (1974) Shadow prices, market wages, and labor supply. Econometrica 42:679–694
Heckman JJ (1979) Sample selection bias as a specification error. Econometrica 47:153–161
Heckman JJ, Honore BE (1990) The empirical content of the Roy model. Econometrica 58:1121–1149
Huber M, Mellace G (2014) Testing exclusion restrictions and additive separability in sample selection
models. Empir Econom 47:75–92
Lahiri S (2021) Necessary and sufficient conditions for variable selection consistency of the LASSO in highdimensions. Ann Stat 49(2):820–844

123

DS-HECK: double-lasso estimation of Heckman selection model

3195

Leeb H, Pötscher BM (2008) Can one estimate the unconditional distribution of post-model-selection
estimators? Econom Theor 24:338–376
Leeb H, Pötscher BM (2008) Sparse estimators and the oracle property, or the return of Hodges’ estimator.
J Econom 142:201–211
Liu H, Xu X, Li JJ (2020) A Bootstrap lasso + partial ridge method to construct confidence intervals forparameters in high-dimensional sparse linear models. Stat Sin 30:1333–1355
Mroz TA (1987) The sensitivity of an empirical model of married women’s hours of work to economic andstatistical assumptions. Econometrica 55:765–799
Mulligan CB, Rubinstein Y (2008) Selection, investment, and women’s relative wages over time. Q J Econ
123:1061–1110
Neumark D, Korenman S (1994) Sources of bias in women’s wage equations: results using sibling data. J
Hum Resour 29:379–405
Neyman J (1959) Optimal asymptotic tests of composite statistical hypotheses. In: Grenander U (ed) Probability and statistics: the Harald Cramer, pp 213–234
Pötscher BM, Leeb H (2009) On the distribution of penalized maximum likelihood estimators: the LASSO,
SCAD, and thresholding. J Multivar Anal 100:2065–2082
Prokhorov A, Schmidt P (2009) GMM redundancy results for general missing data problems. J Econom
151:47–55
Robinson PM (1988) Root-n-consistent semiparametric regression. Econometrica 56:931–954
Roy AD (1951) Some thoughts on the distribution of earnings 1. Oxf Econ Pap 3:135–146
Sampford MR (1953) Some inequalities on Mill’s ratio and related functions. Ann Math Stat 24:130–132
Tibshirani R (1996) Regression shrinkage and selection via the Lasso. J Roy Stat Soc Ser B (Methodol)
58:267–288
Vella F (1998) Estimating models with sample selection bias: a survey. J Hum Resour 33:127–169
Wooldridge JM (2010) Econometric analysis of cross section and panel data. MIT Press
Wuthrich K, Zhu Y (2021) Omitted variable bias of lasso-based inference methods: a finite sample analysis.
Rev Econ Stat 20:1–47
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps
and institutional affiliations.

123