Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
117,526 characters · 20 sections · 50 citation commands
Identification and Estimation of Semiparametric Multilayered Sample Selection Models
heckman1974shadow, heckman1979sample pioneered the econometric analysis of selection bias by modeling the joint determination of participation and outcomes. Heckman's insight has become one of the most influential ideas in microeconometrics. Yet Heckman's framework treats selection as a binary event: an individual either participates or not. In many empirical settings, the selection structure is considerably richer. For instance, conditional on entering the labor force, workers sort into specific occupations, firms, or industries. The outcome of interest is shaped by two layers of selection: participation and sorting among alternatives. Ignoring the sorting layer conflates within-occupation wage effects with between-occupation composition effects, potentially distorting policy implications.
This paper develops semiparametric models for multilayered selection that achieve point identification by leveraging nonlinearity in the selection mechanism rather than exclusion restrictions. I consider two distinct selection architectures. Vertical sorting arises when categories can be meaningfully ordered (for example, by job quality, firm size, or amenity provision) and can be modeled through ordered threshold-crossing processes. Horizontal sorting arises when categories are unordered (such as STEM vs. non-STEM jobs) and requires a multinomial choice framework. For each architecture, I characterize the resulting selection bias, establish conditions under which the outcome equation parameters are point identified, and propose computationally tractable sieve-based estimators.
The central identification challenge in multilayered selection is that the selection bias function generally depends on multiple selection indices: the threshold functions delineating categories or the utility indices governing multinomial choice. This substantially complicates identification relative to the binary selection case, where the bias can be a function of scalar selection probability. I establish identification results for both selection architectures. First, when the ordered selection process is governed by a single index, the selection bias reduces to a function of that index, so that a single continuously distributed covariate together with nonlinearity in the selection index suffices for point identification without an exclusion restriction. When the ordered selection process is fully nonparametric, the selection bias becomes a function of two indices simultaneously. I show that at least three continuous covariates are required to identify the outcome parameters, and that the requirement is binding by exhibiting an explicit non-identification result when fewer continuous covariates are available. This sharp increase in the identification requirement is a consequence of the richer index structure and is, to my knowledge, new in the literature.
Second, for horizontal sorting modeled as multinomial choice with $K+1$ categories (including an outside option like unemployment), the bias correction function generally depends on $K$ indices. I show that additional structural restrictions on the preference heterogeneity reduce the dimensionality of the selection bias. Under a multinomial logit selection, the bias collapses to a single-index control function. A distinctive feature of the multinomial logit specification is that the nonlinearity condition for identification is automatically satisfied even with a linear utility specification. This contrasts with the ordered case, where the control function is approximately linear under Gaussian errors, making identification fragile without exclusion restrictions. Under a weaker exchangeability condition on the taste shocks, I exploit the theory of symmetric polynomials to approximate the bias by a function of a small number of elementary symmetric polynomials of the choice probabilities, providing a practical dimensional reduction that makes semiparametric estimation feasible even with moderately many choice categories.
For estimation, I propose two-step sieve plug-in estimators that can be implemented using standard software. The first step estimates the selection equation nonparametrically using sieves; the second step includes the estimated control functions as nonparametric regressors in the partially linear outcome regression, with heteroskedasticity-robust standard errors. I establish $\sqrt n$-consistency and asymptotic normality under an $o_p(n^{-1/4})$ rate condition on the first-stage sieve estimation. Monte Carlo simulations across seven data-generating processes confirm that the sieve estimators achieve near-oracle performance with correct coverage, and the corrections are robust to weak nonlinearity.
I apply the proposed framework to estimate the gender wage gap among college graduates in South Korea using the Graduates Occupational Mobility Survey (GOMS). Three selection architectures are implemented: an ordered model for firm-size sorting, and multinomial models for field (STEM vs.\ non-STEM) and sector (public vs.\ private) sorting. The uncorrected female hourly-wage penalty is around 5--6 log points in SMEs, non-STEM jobs, and the private sector, around 4 log points in large firms, and near zero (about 1 log point) in STEM and the public sector. After correction in the ordered model, the large-firm female coefficient turns positive, indicating a small conditional premium rather than a penalty, while the SME penalty is little changed. In the field and sector models the correction is modest. The corrected gap has also narrowed over the sample period (2008--2019). These empirical findings connect to a large literature on the gender wage gap, which has long recognized that selection into employment and across occupations is a first-order concern for measuring the gap.\footnote{See, for example, neal2004measured, olivetti2008unequal, blau2017gender. blau2024selection provide recent evidence that selection into employment substantially affects measured gender wage gaps in the United States, finding that correcting for selection narrows the gap by 15--20% over a four-decade period.} Existing corrections in this literature typically either impose parametric distributional assumptions mulligan2008selection or settle for partial identification blundell2007changes, lee2009bounds. The current application adds the intensive-margin sorting layer that these studies abstract from.
This paper also connects to several strands of the econometrics literature. In the binary selection setting, heckman1974shadow, heckman1979sample established the foundational control function approach under joint normality, while chamberlain1986asymptotic, ahn1993semiparametric, powell1989semiparametric, newey1990semiparametric, newey2009two, and das2003nonparametric developed semiparametric/nonparametric alternatives, all requiring exclusion restrictions. Recent work has pursued identification without excluded variables through partial identification lee2009bounds, honore2020selection, heteroscedasticity lewbel2007endogenous, klein2010heteroskedastic, and functional form variation: escanciano2016identification showed that nonlinearity in the selection mechanism can substitute for an exclusion restriction, and kim2025point applied this to the semiparametric selection model with a linear outcome equation. pan2024locally develop a related strategy using debiased machine learning. My paper generalizes the existing frameworks to multilayered selection, where fundamentally new identification arguments are needed when the bias depends on multiple indices.
For multinomial selection, lee1983generalized coupled a logit selection specification with joint normality, dubin1984econometric relaxed the outcome error distribution, and dahl2002mobility introduced a semiparametric polynomial correction. bourguignon2007selection compared and extended these approaches. All require exclusion restrictions or parametric distributional assumptions. sheng2025social exploit exchangeability and elementary symmetric polynomials for dimensionality reduction in social interaction models; I adapt this device to the multinomial sample selection setting. The Roy model tradition roy1951some, heckman1990varieties, french2011identification, heckman2018unordered provides the theoretical foundation for sorting across sectors. dhaultfoeuille2013inference obtain an identification result for a binary extended Roy model in which selection is driven by both potential earnings and a non-pecuniary cost. Exploiting that additive structure together with continuity of at least one covariate, they point-identify the non-pecuniary component without exclusion restrictions or large-support conditions. Identifying the covariate effects on sector-specific earnings, however, still requires either an exclusion restriction or an identification-at-infinity argument on potential outcomes. Most directly related is kroft2024lee, who extend lee2009bounds's nonparametric bounds to multilayered settings under conditional selection monotonicity. Their bounds impose minimal structural assumptions but can be wide in practice. My approach trades that nonparametric generality for point identification.
The remainder of the paper is organized as follows. Section (ref) introduces the multilayered selection framework and establishes identification. Section (ref) presents the estimators and their asymptotic properties. Section (ref) evaluates finite-sample performance through simulations. Section (ref) applies the method to gender wage gaps among Korean college graduates. Section (ref) concludes. Technical proofs and additional simulation and empirical results are provided in the Appendix.
Consider a population of individuals indexed by $i$, each characterized by observable covariates $X_i \in \mathcal{X} \subseteq \mathbb{R}^{d_X}$, a discrete selection indicator $D_i \in \mathcal{C} := \{0, 1, \ldots, K\}$, and potential outcomes $\{Y_{ik}^*\}_{k=1}^K$. $D_i$ encodes two layers of choice: $D_i = 0$ indicates non-participation (e.g., unemployment), while $D_i = k$ for $k \geq 1$ indicates participation in category $k$. The potential outcomes for each $k$ and the observed outcome $Y_i$ are determined by
where $\beta_k \in \mathbb{R}^{d_X}$ is a category-specific parameter vector, $\alpha_k$ is a category-specific intercept, $V_{ik}$ is unobserved heterogeneity with $E[V_{ik}|X_i] = 0$, and $\mathbf{1}[\cdot]$ is the indicator function. Conditional on selection into category $k$, the expected observed outcome is
where the selection bias captures the systematic difference in unobservable characteristics between individuals who select into category $k$ and the population average. The object of interest is $\beta_k$, the effect of covariates on potential outcomes within category $k$.
The key modeling challenge is to specify the selection process generating $D_i$ in a way that (i) allows the selection bias to be characterized by a tractable control function, (ii) permits point identification of $\beta_k$ without exclusion restrictions, and (iii) accommodates the economic structure of occupational sorting. I consider two main architectures in turn: vertical and horizontal sorting.
Suppose the categories in $\mathcal{C}$ are vertically differentiated, so that $D_i = k$ if the individual's “quality” index falls in the $k$-th interval of an ordered partition:
where $Z_i$ is a row vector of covariates affecting selection, $\gamma$ is a parameter vector, $\varepsilon_i$ is a mean-zero error term, and $-\infty = c_0 < c_1 < \cdots < c_K < c_{K+1} = \infty$ are threshold parameters. When $\varepsilon_i \sim N(0,1)$, this is the standard ordered probit model. The selection bias takes a known parametric form under joint normality of $(V_{ik}, \varepsilon_i)$ that are independent of $(X_i, Z_i)$, with $\text{Corr}(V_{ik}, \varepsilon_i) = \rho_k$ and $\text{Var}(V_{ik})^{1/2} = \sigma_k$ as follows:
where $\phi(\cdot)$ and $\Phi(\cdot)$ denote the standard normal p.d.f.\ and c.d.f.\ respectively. This generalizes the inverse Mills ratio in Heckman's binary model and provides a category-specific control function $\lambda(z\gamma; c_k, c_{k+1})$ that can be plugged into the outcome regression.
In the Heckman model, it is well known that identification without an exclusion restriction is fragile because the inverse Mills ratio is approximately linear over much of the effective support of $Z\gamma$ leung1996choice. This near-collinearity problem is equally severe, and arguably worse, in the ordered case. The control function $\lambda(z\gamma; c_k, c_{k+1})$ can exhibit a nearly linear relationship with the index $z\gamma$ (see Figure (ref) in Appendix (ref)).
To relax joint normality without an exclusion restriction while retaining a tractable structure, I consider a semiparametric specification in which the distribution of $\varepsilon_i$ is parametrically specified but the selection index function is left unrestricted:
where $h: \mathcal{X} \to \mathbb{R}$ is an unknown smooth function. Under this specification, the selection bias conditional on $D_i = k$ becomes a function of the single index $h(x)$:
where $f_\varepsilon$ and $F_\varepsilon$ denote the p.d.f.\ and c.d.f.\ of $\varepsilon_i$. $\lambda_k(\cdot)$ is category-specific because the threshold constants $c_k, c_{k+1}$ are category-specific. Consequently, the conditional mean of the observed outcome takes the partial linear form:
which is precisely the structure analyzed in kim2025point. Under standard regularity conditions therein, $\beta_k$ and $\lambda_k$ are identified for each $k = 1, \ldots, K$.\footnote{The regularity conditions require continuous variation in at least one covariate, smoothness of $\lambda_k$, no perfect multicollinearity, and a nonlinearity condition on the composite selection probability $p_k(x) = F_\varepsilon(c_{k+1} - h(x)) - F_\varepsilon(c_k - h(x))$. The intercept $\alpha_k$ is not separately identified from $\lambda_k$ so is normalized to $0$.}
The semiparametric specification (ref) restricts all threshold functions to shift in parallel through the common index $h(x)$. This is a substantive restriction as it requires, for example, that a covariate that makes an individual more likely to surpass the threshold into category 2 also makes them more likely to surpass the threshold into category 3, and by the same amount in index units. To relax this restriction, I consider a fully nonparametric ordered selection model following chesher2012iv:
where $U_i$ is normalized to $\text{Unif}(0,1)$, $X_i \perp\!\!\!\perp U_i$, $h_0(x) = 0$, $h_{K+1}(x) = 1$, and $0 < h_1(x) < h_2(x) < \cdots < h_K(x) < 1$ for all $x$ in the support. $h_k(\cdot)$ is now free to depend on $x$ in an unrestricted manner. This model nests (ref) as a special case.
The threshold functions are nonparametrically identified from the choice probabilities. Defining $\pi_j(x) := P[D_i = j | X_i = x]$, we have
Under the nonparametric specification, the selection bias conditional on $D_i = k$ becomes:
Crucially, the bias function now depends on two indices rather than the single index $h(x)$ that arose in the semiparametric case. This doubled index structure fundamentally changes the identification problem. Define the threshold mapping $H_k(x) := (h_k(x), h_{k+1}(x)) \in [0,1]^2$. The conditional mean of the outcome can then be written as:
I first establish that identification fails generically when $H_k$ is injective.
Injectivity of $H_k$ allows defining a matching $\tilde{\lambda}$ for any candidate $\beta$. Since $H_k: \mathbb{R}^{d_c} \to \mathbb{R}^2$ is generically injective when $d_c \leq 2$ and generically non-injective when $d_c \geq 3$, identification requires at least three continuous covariates, a sharp contrast to binary selection, where one suffices.\footnote{A smooth map from $\mathbb{R}^{d_c}$ to $\mathbb{R}^2$ with full-rank Jacobian is locally injective when $d_c \leq 2$ but has $(d_c - 2)$-dimensional fibers when $d_c \geq 3$, making injectivity generically impossible.} The following assumption collects the conditions under which identification can be established with three continuous covariates.
Assumption (ref) (ii)--(iii) require smoothness of selection indices and selection bias functions, and (iv) requires that $h_k$ and $h_{k+1}$ respond to the three continuous covariates in linearly independent directions. This rules out the case where $h_{k+1}$ is merely a parallel shift of $h_k$ by a constant (as in the semiparametric model). Assumption (ref)(v) is the key nonlinearity condition. Each Jacobian $J_k(x)$ has a two-dimensional column space in $\mathbb{R}^3$, and hence a one-dimensional left null space. The condition requires that these one-dimensional null spaces, evaluated at three points, collectively span $\mathbb{R}^3$. This is generically satisfied when the threshold functions exhibit sufficient nonlinearity. Assumption (ref)(vii) is a mild topological regularity condition: the selected sample explores the threshold-index space as a single connected region rather than as several isolated pieces. A simple sufficient condition is that $\operatorname{supp}(X_i \mid D_i = k)$ itself is connected. Now I establish identification of $\beta_k$ and $\lambda_k$ in the following proposition.
When both threshold indices take linear form $h_k(x) = f_k(x'\gamma_k)$, the column space of $J_k(x)$ is a fixed plane in $\mathbb{R}^3$ regardless of $x_c$, violating Assumption (ref)(v). Continuous excluded variables can substitute for covariates, but the substitution rate is non-uniform. The first excluded variable replaces one continuous covariate, leaving the requirement $d_c^{\text{covariates}} + d_c^{\text{excluded}} \geq 3$ when $d_c^{\text{excluded}} \in \{0, 1\}$. The second excluded variable, however, replaces two continuous covariates: once $d_c^{\text{excluded}} \geq 2$, identification follows from the exclusion-based route in das2003nonparametric, and no continuous covariate is needed. The required number of continuous covariates is therefore three, two, and zero as $d_c^{\text{excluded}}$ moves from zero to one to two, reflecting the regime switch from nonlinearity-based identification to exclusion-based identification.
When the categories in $\mathcal{C}$ are horizontally differentiated (for example, broadly defined occupations or industries with no natural ordering), the ordered threshold-crossing framework is inappropriate. Instead, I model selection as a utility-maximizing multinomial choice:
where $u_k(X_i)$ is the deterministic utility component for category $k$, and $\varepsilon_{ik}$ captures preference heterogeneity. The individual chooses the category that maximizes utility.\footnote{The classical Roy model roy1951some is a special case of (ref) in which $u_k(X_i) + \varepsilon_{ik} = Y_{ik}^* = \alpha_k + X_i\beta_k + V_{ik}$, so that selection is driven by potential wages alone. The generalized Roy model heckman1990varieties is the case $u_k(X_i) + \varepsilon_{ik} = Y_{ik}^* - C_{ik}$, allowing a non-pecuniary cost $C_{ik}$. The present framework accommodates both without requiring utility and potential wages to coincide.} This framework is also considered in kroft2024lee as a parametric special case of their nonparametric multilayered selection model. Given the selection rule (ref), the selection bias conditional on $D_i = k$ is
where $\delta_{kj}(x) := u_k(x) - u_j(x)$ for $j \neq k$ and $\lambda_k: \mathbb{R}^K \to \mathbb{R}$ is an unknown function whose argument is the $K$-vector of pairwise differences indexed over the non-chosen alternatives. In general, the bias depends on $K$ indices (one for each pairwise comparison), creating a severe curse of dimensionality as the number of categories grows.
An equivalent representation expresses the bias as a function of the choice probabilities rather than the utility differences. This reformulation is not required for identification, but it is practically useful in estimation.\footnote{Choice probabilities are compactly supported on $[0,1]$, whereas utility differences range over $\mathbb{R}$. Sieve smoothers approximating the selection bias function in the second stage behave substantially better in finite samples with a bounded argument than with an unbounded index.} The equivalence, which holds under mild regularity conditions, provides the conceptual foundation for using estimated choice probabilities as control functions in the second stage.
The idea that choice probabilities can invert latent utility indices is well established in discrete choice and demand analysis.\footnote{For instance, hotz1993conditional derive CCP inversion in dynamic logit models, berry1994estimating establishes inversion in the multinomial logit demand system, and berry2013connected provide a more general demand-side invertibility result under connected substitutes.} Proposition (ref) establishes a general inversion result for the present static multinomial selection model under arbitrary continuous strictly positive joint densities of the preference shocks. This, in turn, justifies using estimated choice probabilities as control-function arguments in the sample-selection outcome equation and provides a formal foundation for the probability-based correction of dahl2002mobility, who derives a single-index reduction under a maintained index-sufficiency assumption.
The conditional mean of the observed outcome now takes the partially linear form with $K$ nonparametric indices:
where $p_0(x) = 1 - \sum_{j=1}^K p_j(x)$ is determined by the simplex constraint. Identification of $\beta_k$ in this $K$-index model follows from the same argument as the nonparametric ordered case (Proposition (ref)), generalized from two to $K$ indices.
This result establishes a sharp trade-off between structural restrictions and covariate requirements: without additional restrictions, identification requires at least $K + 1$ continuous covariates. When $K$ is moderate, this is feasible; when $K$ is large, the requirement becomes prohibitive. The remainder of this section develops two structural restrictions that reduce the dimensionality of $\tilde{\lambda}_k$ and correspondingly lower the covariate requirement.
Suppose $u_k(X_i)=X_i\gamma_k$ and $\varepsilon_{i0}, \ldots, \varepsilon_{iK}$ are independently and identically distributed as standard Extreme Value Type I (Gumbel). $\gamma_0$ is normalized to $0$ for identification. I further assume that $V_{ik}$ depends on $(\varepsilon_{i0}, \ldots, \varepsilon_{iK})$ only through $\varepsilon_{ik}$:
In the occupational sorting context, this assumption means that the worker's unobserved productivity in occupation $k$ only depends on their taste for that occupation. This is natural under an occupation-specific match quality interpretation: a worker with a strong affinity for a particular occupation ($\varepsilon_{ik}$ large) tends to also be productive in that occupation ($V_{ik}$ large), because taste and ability for a specific occupation are correlated. Once conditioned on $\varepsilon_{ik}$, the preference shocks for the other occupations $\varepsilon_{ij}$, $j \neq k$, carry no additional information about their productivity in $k$. While this restriction is substantive, it substantially simplifies the identification analysis.
Under these assumptions, the conditional distribution of $\varepsilon_{ik}$ given $D_i = k$ and $X_i = x$ admits a simple characterization. Let $U_{ij} = x\gamma_j + \varepsilon_{ij}$ denote the utility from category $j$, and let $M = \max_{0 \leq j \leq K} U_{ij}$. Since $\varepsilon_{ij} \sim \text{EV}(0,1)$ independently, the cumulative distribution of the maximum is $$F_M(m | x) = \prod_{j=0}^K \exp\{-e^{-(m - x\gamma_j)}\} = \exp\{-A(x) e^{-m}\},$$ where $A(x) := \sum_{j=0}^K e^{x\gamma_j} = 1 + \sum_{j=1}^K e^{x\gamma_j}$ is the logit denominator. Hence $M | x \sim \text{EV}(\ln A(x), 1)$. Conditional on $D_i = k$ and $X_i = x$, the winner's utility equals the maximum: $U_{ik} = M$. By the well-known property of Gumbel random variables that the conditional distribution of the maximum given which alternative wins depends on $x$ only through $\ln A(x)$, we have $U_{ik} | D_i = k, X_i = x \sim \text{EV}(\ln A(x), 1)$. Since $\varepsilon_{ik} = U_{ik} - x\gamma_k$, a location shift gives
The conditional density of $\varepsilon_{ik}$ given selection into $k$ thus depends on $x$ only through the scalar index $\nu_k(x) := \ln A(x) - x\gamma_k$. Applying the restriction (ref):
The selection bias reduces to a function of the single index $\nu_k(x) = \ln(1 + \sum_{j=1}^K e^{x\gamma_j}) - x\gamma_k$. The conditional mean of the outcome is therefore
which is again a partial linear model with a single-index control function.
Since (ref) has the same partially linear structure as the semiparametric ordered model, identification of $\beta_k$ and $\lambda_k$ follows from kim2025point under the same regularity conditions. A distinctive feature of the multinomial logit model is that the nonlinearity condition is automatically satisfied, even with a linear specification $u_k(x) = x\gamma_k$ for the deterministic utility. The inclusive value $\nu_k(x)$ is inherently nonlinear in $x$ because the log-sum-exp function is convex and not affine whenever at least two $\gamma_j$ differ (see Appendix (ref)). This contrasts sharply with the ordered case, where identification without exclusion restrictions is fragile. In the multinomial logit model, identification without exclusion restrictions is robust because the nonlinearity is a structural consequence of the multinomial choice mechanism. Figure (ref) in Appendix (ref) illustrates this inherent nonlinearity for several parameter configurations.
The preceding analysis relies on two restrictions: the i.i.d.\ Gumbel assumption on the preference shocks, which implies the independence of irrelevant alternatives (IIA), and the own-shock restriction (ref). IIA rules out correlation across alternatives, which is implausible when some occupations are closer substitutes than others. The own-shock restriction fails whenever unobserved ability has a general component that affects both preferences and productivity across occupations. To accommodate richer dependence among the preference shocks, and between the preference shocks and outcome errors, while maintaining a tractable selection bias structure, I introduce an exchangeability condition under the more general nonparametric utility specification (ref).
This assumption does not impose IIA, and it replaces the own-shock restriction with a weaker symmetry requirement that allows $V_{ik}$ to depend on all preference shocks $(\varepsilon_{i0}, \ldots, \varepsilon_{iK})$, provided this dependence is symmetric in the non-chosen alternatives. Under this assumption, for an individual contemplating category $k$, the alternative categories are ex ante symmetric in their unobserved preference. The exchangeability condition is more plausible with broadly defined categories (STEM vs.\ non-STEM) than with finely disaggregated ones (specific occupations within STEM), where nested structures are more natural. Under exchangeability, the selection bias function $\tilde{\lambda}_k(p_0(x), \ldots, p_K(x))$ is symmetric in its non-chosen probabilities. The following proposition exploits this symmetry to reduce the dimensionality of the probability-based bias correction.
This proposition provides a principled truncation strategy for the selection bias in terms of observable choice probabilities. For a given truncation order $L \leq K$, one can approximate $\tilde{\lambda}_k$ by a function of only $(e_1, \ldots, e_L)$, discarding higher-order elementary symmetric polynomials. When $L$ is small relative to $K$, this yields a substantial dimensionality reduction. In practice one truncates the polynomial approximation at some order $L \leq K$ in the elementary symmetric polynomials, treating $\tilde{\lambda}_k(p_0, \ldots, p_K)$ as if it were exactly a function $\breve{\lambda}_k^{(L)}(e_1, \ldots, e_L)$ of only $L$ arguments. The conditional mean of the outcome under this $L$-truncated working model is $$E[Y_i \mid X_i = x, D_i = k] \approx x\beta_k + \breve{\lambda}_k^{(L)}\bigl(e_1(p(x)), \ldots, e_L(p(x))\bigr).$$ Write $E_L(x) := (e_1(p(x)), \ldots, e_L(p(x)))$ for the vector of the first $L$ elementary symmetric polynomials of the non-chosen probabilities. The identification of $\beta_k$ under this approximation follows from the multi-index structure of the elementary symmetric polynomials. For $\ell = 1$, $e_1(p(x)) = 1 - p_k(x)$ is nonlinear in $x$ whenever $p_k(x)$ is nonlinear, which holds generically; identification then follows from this single-index nonlinearity with one or two continuous covariates. For $\ell \geq 2$, $e_\ell$ is a polynomial of degree $\ell$ in the choice probabilities, and under the maintained rank and spanning conditions of the proposition below it supplies the multi-index nonlinearity needed for identification. For $L = 2$ the resulting two-index structure parallels the nonparametric ordered case and requires three continuous covariates. The general pattern is:
Without structural restrictions, the dimensionality of the bias grows with $K$, so in practice researchers should use a small number of broadly defined categories or impose the logit or exchangeability restriction.
This section develops two-step sieve plug-in estimators for each selection architecture and establishes their asymptotic properties. All models in the previous section produce a conditional mean of the form
where $g_k: \mathcal{X} \to \mathbb{R}^L$ is a vector of $L$ indices and $\lambda_k: \mathbb{R}^L \to \mathbb{R}$ is an unknown smooth function. The number of indices $L$ and the structure of $g_k$ depend on the selection architecture. The same second stage is used across all architectures and is described first; let $\hat{g}_k(\cdot)$ denote the first-stage estimator obtained from the selection data $\{(D_i, X_i)\}_{i=1}^n$. For each category $k = 1, \ldots, K$, use the subsample $\mathcal{I}_k := \{i : D_i = k\}$ with $n_k := |\mathcal{I}_k|$. Approximate $\lambda_k(\cdot)$ by a sieve basis and estimate $\beta_k$ jointly:
where $B_{J_n}^{(L)}(\cdot)$ is a sieve basis of dimension $\kappa_n$ for the $L$-variate function $\lambda_k$, and $\delta_k$ is the vector of sieve coefficients. The estimator of $\beta_k$ is the ordinary least squares (OLS) coefficient on $X_i$ from the regression (ref).
I specify the sieve basis as follows. For $L = 1$ (single-index models), let $B_{J_n}^{(1)}(s) = (B_1(s), \ldots, B_{J_n}(s))'$ be a univariate B-spline basis of order $r$ with $J_n$ interior knots placed on the support of $\hat{g}_k(X_i)$. The sieve dimension is $\kappa_n = J_n + r$. For $L = 2$ (nonparametric ordered model), use a tensor product B-spline basis:
where $B^{[\ell]}$ denotes a univariate B-spline basis with $J_n^{[\ell]}$ interior knots in the $\ell$-th direction. The sieve dimension is $\kappa_n = (J_n^{[1]} + r)(J_n^{[2]} + r)$. For $L \geq 3$, the tensor product construction extends naturally, with dimension $\kappa_n = \prod_{\ell=1}^L (J_n^{[\ell]} + r)$. Define the augmented regressor vector for observation $i \in \mathcal{I}_k$: $$W_{ik} := \left(X_i, \, B_{J_n}^{(L)}(\hat{g}_k(X_i))'\right)' \in \mathbb{R}^{d_X + \kappa_n},$$ and write the second-stage regression compactly as $Y_i = W_{ik}' \theta_k + \eta_{ik}$ where $\theta_k = (\beta_k', \delta_k')' \in \mathbb{R}^{d_X + \kappa_n}$. The OLS estimator is
and $\hat{\beta}_k$ is the subvector of $\hat{\theta}_k$ corresponding to $X_i$.
In the following subsection, I present details of the first stage estimation procedure for each selection architecture.
Let $Q_n(x) = (q^{(1)}(x), \ldots, q^{(q_n)}(x))'$ be a sieve basis. For the semiparametric ordered model, let $\theta$ denote the parameter vector $(\alpha, c_1, \ldots, c_K)$. The sieve MLE maximizes the ordered choice log-likelihood:
yielding $\hat{h}(x) = Q_n(x)'\hat{\alpha}$. With $F_\varepsilon = \Phi$, this is a sieve ordered probit. For the nonparametric ordered model, define $\tilde{D}_{ik} = \mathbf{1}[D_i \leq k-1]$ and estimate each threshold function by sieve logistic regression:
giving $\hat{h}_k(x) = \Lambda(Q_n(x)'\hat{\alpha}_k)$ and index vector $\hat{g}_k(x) = (\hat{h}_k(x), \hat{h}_{k+1}(x))$.\footnote{The separate estimation does not automatically enforce $\hat{h}_1(x) < \cdots < \hat{h}_K(x)$; violations are rare in practice and can be corrected by rearrangement chernozhukov2010quantile.}
For the multinomial logit model, the workhorse first-stage estimator is a sieve multinomial logit: replace the linear index $X_i \gamma_k$ with a flexible sieve approximation $Q_n(X_i)' \alpha_k$ and maximize the MNL log-likelihood
yielding $\hat{\nu}_k(x) = \ln \sum_{j=0}^K e^{Q_n(x)' \hat{\alpha}_j} - Q_n(x)' \hat{\alpha}_k$ and $\hat{p}_j(x) = \text{softmax}(Q_n(x)'\hat{\alpha})_j := e^{Q_n(x)'\hat{\alpha}_j} / \sum_{\ell=0}^K e^{Q_n(x)'\hat{\alpha}_\ell}$, with the baseline normalization $\hat{\alpha}_0 = 0$. For the exchangeability model, the same sieve MNL is used; the estimated choice probabilities are then used to form the elementary symmetric polynomials of the non-chosen probabilities:
Here I derive asymptotic properties of the proposed estimators under regularity conditions. I state the assumptions in a unified manner, applying to all selection architectures.
Let $\mathcal{H}^m(\mathcal{S})$ denote the H\"older ball of order $m$ on $\mathcal{S}$: the set of functions whose partial derivatives up to order $\lfloor m \rfloor$ are bounded and whose $\lfloor m \rfloor$-th partial derivatives are H\"older continuous of order $m - \lfloor m \rfloor$.
The following lemma establishes sufficient conditions for the first stage rate condition (Assumption (ref)) in each selection architecture.
This rank condition is guaranteed by the identification conditions of Section (ref) and those of Propositions 1--3 in kim2025point.
The following theorem establishes the $\sqrt{n}$-consistency and asymptotic normality of $\hat{\beta}_k$.
The first-stage estimation error is asymptotically negligible under (ref), so the feasible estimator has the same asymptotic distribution as the oracle using the true $g_{k0}$. Under homoskedasticity, the asymptotic variance simplifies to $V_k = \sigma_k^2 \Sigma_k^{-1}$. The following corollary establishes that the sieve plug-in estimator achieves the semiparametric efficiency bound under homoskedasticity.
The FWL residualization $\tilde{X}_{ik} = X_i - \Pi_k(g_{k0}(X_i))$ is the efficient influence function for the partially linear model, and the sieve basis consistently approximates the projection $\Pi_k$. Models with fewer indices ($L = 1$) generically yield smaller asymptotic variance than those with more, reflecting the statistical cost of the structural restrictions needed to reduce the dimensionality of the bias function.\footnote{Under heteroskedasticity, $V_k = \Sigma_k^{-1} \Omega_k \Sigma_k^{-1}$ remains valid but $\hat{\beta}_k$ is no longer efficient. Efficient estimation under heteroskedasticity would require weighted least squares with estimated conditional variance $\sigma_k^2(X_i)$.}
From the asymptotic distribution, I derive the consistent variance estimator in the following theorem.
As the estimator is a plain OLS on the augmented design $W_{ik} = (X_i, B^{(L)}(\hat{g}_k(X_i))')'$, the inference object for $\beta_k$ is the heteroskedasticity-robust sandwich variance $\hat{V}_k$ from (ref), obtained directly using standard software.
This section provides Monte Carlo evidence on the finite-sample performance of the proposed estimators across selection architectures. The baseline simulations use $n = 5{,}000$ observations and $500$ replications, with cubic splines throughout. For each data-generating process (DGP) and estimator, I report the root mean squared error (RMSE), mean bias, and empirical coverage probability of the 95% confidence interval.
Two DGPs test ordered selection models with $K = 2$ occupation categories. Four estimators are compared in the main text: (i) OLS, which ignores selection; (ii) Linear, a parametric selection correction using a linear ordered probit; (iii) Oracle, an infeasible benchmark using the true correction function in the second stage; and (iv) Sieve, the proposed estimator using a nonparametric first stage and a sieve approximation to $\lambda_k(\hat{h}_k(x), \hat{h}_{k+1}(x))$ in the second stage.
\noindentDGP1: Two continuous covariates. The selection mechanism follows an ordered threshold-crossing model with a nonlinear index: $$D_i = k \quad \text{if} \quad c_k \leq 0.5 X_i - 0.5 X_i^2 + 0.2 X_i^3 + 0.5 X_i Z_i + Z_i - 0.5 Z_i^2 + U_i < c_{k+1},$$ with thresholds $c_1 = -1.5$ and $c_2 = 0.5$. The covariates $(X_i, Z_i) \sim N(\mathbf{0}, I_2)$ are independent, and the errors $(U_i, V_{i1}, V_{i2})$ are jointly normal with correlation 0.75 between the selection error and each outcome error. The outcome equations are $Y_{i1} = 0.5 + 0.5 X_i + 0.25 Z_i + V_{i1}$ and $Y_{i2} = 0.6 + 0.7 X_i + 0.5 Z_i + V_{i2}$.
\noindentDGP2: Mixed covariates. This DGP modifies DGP1 by replacing the continuous covariate $Z_i$ with a binary indicator $Z_i = \mathbf{1}[Z_i' > 0]$ where $Z_i' \sim N(0,1)$, and enriching the selection index with interaction terms $X_i^2 Z_i$ and $X_i^3 Z_i$: $$\tilde{h}(X_i, Z_i) = -0.2 X_i - 0.5 X_i^2 + 0.3 X_i^3 + 0.1 X_i Z_i + 0.5 Z_i - 0.3 X_i^2 Z_i + 0.2 X_i^3 Z_i.$$
Table (ref) reports the occupation-level RMSE, absolute bias, and coverage across DGPs. In DGP1, the OLS estimator exhibits substantial bias, confirming large selection bias. The Linear estimator fails catastrophically for Occupation 1, reflecting the severe mismatch between the linear index and the true nonlinear index. For Occupation 2, the Linear estimator is less extreme but still substantially biased. The Oracle estimator is nearly unbiased with tight dispersion and the Sieve estimator tracks the oracle very closely, with RMSE only slightly larger than the infeasible benchmark. In DGP2 the binary covariate complicates the within-category support of the first-stage probability vector, but the Sieve estimator remains near-oracle on both occupations. Coverage is at or near the nominal 95% level for Oracle and Sieve in all four cells.
Additional simulation results with $K=3$ are reported in Appendix (ref): Ordered DGP3 ($K = 3$, three continuous covariates) showcases near-oracle Sieve performance at full nonlinearity and stress-tests weak nonlinearity by scaling the higher-order terms of the selection index by $\delta$. As the selection index becomes nearly linear, even the Oracle's variance inflates; the Sieve estimator essentially matches the Oracle for moderate nonlinearity and incurs a finite-sample bias penalty only in the most adversarial near-linear regime, while remaining far better behaved than the parametric Heckman-type estimator throughout.
Four DGPs are considered for multinomial selection: a baseline $K = 2$ design under IIA, two exchangeable designs with $K = 3$, and a non-exchangeable design with $K = 3$. Table (ref) summarizes the key features of the DGPs. Five estimators are compared in the main text: (i) OLS, which ignores selection; (ii) MLogit, using the sieve-estimated inclusive value $\hat{\nu}_k(x)$ as a single-index control function under the own-shock restriction; (iii) Oracle, an infeasible benchmark using the true inclusive value as the linear correction (DGPs 1--2, where the own-shock restriction holds) or the true choice-probability vector with the Sieve second-stage (DGPs 3--4, where own-shock fails); (iv) Sieve, using the sieve MNL predicted probability vector through cubic B-spline marginals plus pairwise tensor interactions in the second stage; and (v) Exch-$L2$ (only for DGPs 2--4 with $K = 3$), using the first two elementary symmetric polynomials $(\hat{e}_1, \hat{e}_2)$ of the sieve MNL choice probabilities under exchangeability.
DGP1 defines the selection procedure with $K = 2$ following utility maximization, $D_i = \operatorname*{arg\,max}_{j \in \{0,1,2\}} \{f_j(X_i, Z_i) + \varepsilon_{ij}\},$ where $f_0 \equiv 0$ and $\varepsilon_{ij} \sim \text{Gumbel}(0,1)$. Both covariates are independently drawn from the standard normal distribution and the utility functions $f_j$ are polynomials of degree 3 in $X$ and degree 2 in $Z$ with interaction terms. The outcome errors satisfy $V_{ik} = \varepsilon_{ik} + \tilde{\varepsilon}_{ik}$ with $\tilde{\varepsilon}_{ik} \sim \text{Gumbel}(0,1)$ independent, inducing correlation between selection and outcome errors. The outcome parameters are $\beta_1 = (0.5, 0.7)'$ and $\beta_2 = (0.8, 0.5)'$. DGP1 satisfies IIA, exchangeability, and the own-shock restriction.
DGP2 features $K = 3$ and three continuous covariates $(X_i, Z_i, W_i) \sim N(\mathbf{0}, I_3)$. The inclusion of three continuous covariates is dictated by the identification theory: by Proposition (ref)(ii), the $L = 2$ truncation of the elementary symmetric polynomial correction requires $L + 1 = 3$ continuous covariates for identification. The utilities are $U_{ij} = f_j(X_i, Z_i, W_i) + \varepsilon_{ij}$ for $j \in \{0, 1, 2, 3\}$, where the utility functions $f_j$ include quadratic terms and pairwise interactions among the three covariates, and the $\varepsilon_{ij}$ are i.i.d.\ Gumbel. The outcome equation for each $k$ is: $Y_{ik} = \alpha_{k} + \beta_{k1} X_i + \beta_{k2} Z_i + \beta_{k3} W_i + V_{ik},$ where $(\alpha_1, \beta_{11}, \beta_{12}, \beta_{13}) = (0.4, 0.5, 0.7, 0.3)$, $(\alpha_2, \beta_{21}, \beta_{22}, \beta_{23}) = (0.6, 0.8, 0.5, 0.4)$, and $(\alpha_3, \beta_{31}, \beta_{32}, \beta_{33}) = (0.5, 0.3, 0.9, 0.6)$. The selection bias for each category is approximated using elementary symmetric polynomials of the choice probabilities, truncated at order $L = 2$ as in the $L$-truncated working model. This DGP tests whether the exchangeability approximation is effective when the number of occupations exceeds two and the required continuous variation for identification is available. IIA, exchangeability, and the own-shock restriction hold in this DGP.
DGP3 is designed to showcase the practical value of the exchangeability framework when the own-shock restriction and IIA fail but exchangeability holds. The preference shocks are equicorrelated normal, $\varepsilon_{ij} = \sqrt{\rho}\, c_i + \sqrt{1 - \rho}\, z_{ij},\ \forall j \in \{0,1,2,3\},$ where $c_i \sim N(0,1)$ is a common factor, $z_{ij} \sim N(0,1)$ are independent, and $\rho = 0.5$. This symmetric structure is exchangeable by construction but violates IIA. The covariates, utility functions, and outcome equations are identical to DGP2. The critical departure from DGP2 is in the outcome errors, which introduce self-reinforcing specialization: $V_{ik} = \varepsilon_{ik} + \gamma_{\text{sr}} \, \varepsilon_{ik} \left( \varepsilon_{ik} - \bar{\varepsilon}_{i,-k} \right) + \sigma_\eta \, \eta_{ik},$ where $\bar{\varepsilon}_{i,-k} = K^{-1} \sum_{j \neq k} \varepsilon_{ij}$ is the mean preference shock of competing alternatives, $\gamma_{\text{sr}} = 1.0$ controls the specialization intensity, $\sigma_\eta = 0.2$, and $\eta_{ik} \sim N(0,1)$. Workers whose preference shock $\varepsilon_{ik}$ exceeds the competition average receive amplified productivity in occupation $k$, creating positive selection. This violates the own-shock restriction because $V_{ik}$ depends on all $\varepsilon_{ij}$ through $\bar{\varepsilon}_{i,-k}$. However, the dependence on competing shocks is symmetric, preserving exchangeability.
Lastly, DGP4 breaks IIA, exchangeability, and the own-shock restriction all together to assess the boundary of the proposed methods while keeping $K = 3$. The preference shocks follow a single-factor model: $\varepsilon_{ij} = \lambda_j \, f_i + u_{ij},\ \forall j \in \{0, 1, 2, 3\},$ where $f_i \sim N(0,1)$ is a common factor, $u_{ij} \sim N(0,1)$ are independent. As the loadings $\boldsymbol{\lambda} = (0, 0.3, 0.8, -0.5)$ are heterogeneous across alternatives, the factor does not cancel in utility differences: $\text{Cov}(\varepsilon_{ij} - \varepsilon_{ik}, \varepsilon_{il} - \varepsilon_{ik}) = (\lambda_j - \lambda_k)(\lambda_l - \lambda_k)$ depends on the identities of $j$ and $l$, not just their count, violating exchangeability. Occupation 2 (loading $\lambda_2 = 0.8$) and occupation 3 ($\lambda_3 = -0.5$) have the most dissimilar loadings and thus the weakest substitutability, while occupation 1 ($\lambda_1 = 0.3$) is closer to the outside option ($\lambda_0 = 0$). The covariates, utility specifications, and outcome equations are identical to DGP2--3. The outcome error follows the same self-reinforcing specialization structure as DGP3, with $\gamma_{\text{sr}} = 2.0$ and $\sigma_\eta = 0.2$.
Table (ref) reports the RMSE, bias, and coverage averaged across all occupations and coefficients within the four DGPs. DGPs 1--2 show a consistent pattern. The OLS is severely biased, whereas the MLogit estimator (correctly specified) exhibits near-zero bias. The oracle provides a tight benchmark with the smallest RMSE that MLogit nearly matches. The Sieve estimator also removes most of the selection bias but has larger RMSE due to the higher dimensionality of its control function. The semiparametric estimators (Sieve, Exch-$L2$) all achieve coverage rates close to the nominal 95% level. In DGP2 with three covariates, the Exch-$L2$ estimator, which uses a more flexible two-dimensional control function, is almost identical to MLogit. Sieve achieves comparable but slightly inflated RMSE.
DGP3 provides the most informative test of the exchangeability framework. Although MLogit's own-shock restriction is violated in DGP3, both MLogit and the correctly specified Exch-$L2$ correct most of the bias and outperform Oracle in terms of RMSE. Sieve performs similarly to Oracle with slightly inflated RMSE. With $K = 3$, the second stage for Sieve and Oracle requires $K + 1 = 4$ continuous covariates for identification. Since only three are available, the probability correction remains under-identified for Sieve and Oracle, confirming the value of the structural restrictions in reducing the effective dimensionality of the bias correction. Coverage is comparable across MLogit, Exch-$L2$ and Sieve.
In DGP4 with non-exchangeability, MLogit achieves the lowest RMSE, with Exch-$L2$ a close second. Both substantially outperform Sieve and Oracle, which remain underidentified yet correctly specified. In terms of bias, Exch-$L2$ leads over MLogit, while Sieve and Oracle exhibit larger bias. The performance is heterogeneous across occupations. For occupation 2 ($\lambda_2 = 0.8$, the most extreme loading), both Sieve and Oracle exhibit substantial bias on $\beta_{21}$ ($-0.231$ and $-0.183$, respectively) with RMSE substantially larger than that of MLogit and Exch-$L2$. The MLogit and Exch-$L2$ approaches deliver stable performance across all occupations because they require fewer identifying covariates (one and three, respectively).
To isolate structural biases from finite-sample variance, supplementary simulations at $n = 100{,}000$ are reported in Table (ref) in Appendix (ref). With variance largely eliminated, Exch-$L2$ overtakes MLogit (RMSE 0.066 vs.\ 0.072). Sieve and Oracle remain substantially worse (RMSE 0.115 and 0.108). These results demonstrate the sharp identification requirement on continuous covariates. Appendix (ref) confirms these findings under additional conditions. Sensitivity analysis at $n = 1{,}000$ and $n = 2{,}000$ confirms that bias reduction is already substantial at moderate sample sizes. A bootstrap validation exercise confirms that the analytical standard errors are a reasonable approximation to bootstrap standard errors.
The Monte Carlo results demonstrate the near-oracle performance of the proposed sieve estimators whenever the identification conditions are met. When the own-shock restriction is violated but exchangeability holds, the Exch-$L2$ estimator, which is then correctly specified, performs almost identically to the Oracle, and MLogit remains nearly as accurate despite relying on the now-violated own-shock restriction. The structural restrictions are valuable. When IIA, exchangeability, and the own-shock restrictions do not hold, MLogit and Exch-$L2$ still perform reasonably well. Sieve remains stable in all of these settings, but it pays its own RMSE penalty when the identification conditions are not met.
I apply the proposed framework to estimate the gender wage gap among recent college graduates in South Korea, using the Graduates Occupational Mobility Survey (GOMS). It is well documented that selection into employment mulligan2008selection, blau2017gender and occupational sorting goldin2014grand are quantitatively important for estimating gender wage gaps. I extend this line of work by explicitly modeling two layers of selection (participation and sorting). South Korea offers a particularly informative setting for multilayered selection. The economy features a pronounced dualism between large conglomerates and small and medium enterprises (SMEs), which motivates the ordered selection framework. The labor market also exhibits horizontal segmentation along occupation field (STEM vs. non-STEM) and sector (public vs. private). Female labor force participation among young graduates remains lower than male participation, and oecd2024korea reports the largest gender wage gap among member countries. These institutional features can generate strong selection on both the extensive and intensive margins.
The GOMS is a nationally representative survey of college graduates. Each cohort is surveyed approximately 18 months after graduation. I pool the 2008--2019 waves and restrict the sample to graduates aged 35 or younger, yielding a full sample of 207,985. The outcome variable is the log hourly wage, constructed as the log of monthly gross earnings divided by total hours (regular plus overtime). The key covariate of interest is a female indicator. The pre-determined controls that enter both the selection and the outcome equations are age, college GPA (on a 0--100 scale), parental income (the midpoint of the reported income bracket), a four-year university indicator (versus two-year colleges), major category (seven groups), university founding type (the institution's ownership category---national, public, private, and so on; six categories), school region (17 administrative units at the city/province level), and survey year fixed effects. The outcome equation additionally controls for job tenure in months and for the realized sorting categories of the other two architectures; these are determined after selection (and are defined only for the employed), so they do not enter the selection equation. No pre-determined covariate is excluded from the outcome equation; identification rests entirely on the nonlinearity of the control function rather than on an exclusion restriction. The three continuous covariates (age, GPA, and parental income) enter the first stage through a flexible sieve specification (cubic penalized regression splines and pairwise tensor products) and enter the outcome equation linearly.
I implement the following architectures for $D_i \in \{0, 1, 2\}$: i) Ordered (firm size) where $D_i = 1$ and $2$ denote employment at a SME and at a large firm ($\geq 300$ employees), respectively; ii) Multinomial (field) where $D_i = 1$ and $2$ denote a non-STEM job and a STEM job, respectively; and iii) Multinomial (sector) where $D_i = 1$ and $2$ denote public-sector and private-sector employment, respectively. $D_i=0$ indicates non-participation. The first architecture reflects the hierarchy in the Korean labor market. In the ordered first stage, an ordered probit with a linear index (including cubic terms of continuous covariates and interactions) is estimated for the parametric control function. For nonparametric first stage, I separately estimate $P(D_i \geq 1 \mid X_i)$ and $P(D_i = 2 \mid X_i)$ which serve as the control function arguments. Unlike firm size, there is no natural ordering in multinomial architectures. In both multinomial architectures, the first stage estimates a sieve multinomial logit as in (ref), from which the inclusive values and choice probabilities are constructed.
Table (ref) reports summary statistics by gender and by the three selection classifications. Several patterns are noteworthy. Women constitute 46.4% of the sample but are underrepresented at large firms and in STEM jobs, and overrepresented in the public sector. Women have modestly higher GPAs but are less likely to hold STEM majors. The raw gender wage gap is substantial: the mean difference in log hourly wages is 12 log points. Across selection categories, large-firm workers earn more than SME workers, STEM workers earn more than non-STEM workers, and private-sector workers earn more than public-sector workers. Non-participants have similar age and parental income to the employed sample, but lower GPAs and fewer four-year university graduates.
For each architecture, I compare several estimators: i) OLS on the selected subsample; ii) a parametric correction: the ordered probit generalized inverse Mills ratio (ordered) or the inclusive value (multinomial) as a single-index control function; and iii) the proposed Sieve estimator, implemented as OLS on a tensor-product cubic B-spline basis in the estimated choice probabilities. In the ordered model, the wage equation for SME workers includes both threshold-probability indices $(\hat{p}_1, \hat{p}_2)$, while the large-firm equation includes only the upper threshold probability $\hat{p}_2$, reflecting the boundary structure of the top category. In the multinomial models, the wage regression uses the inclusive value $\hat{\nu}_k$ as a single-index control in the MLogit estimator, while the Sieve estimator uses the estimated choice probabilities $(\hat{p}_1, \hat{p}_2)$ as controls. Additionally, an Exch-$L1$ estimator uses the first elementary symmetric polynomial ($1 - \hat{p}_k$) as a control assuming exchangeability.
Table (ref) reports the estimated female coefficients across all specifications for both log hourly and log monthly wages. For hourly wage, the results reveal different selection patterns across the three architectures. For firm-size sorting, the OLS gender gap is $-$0.047 in SMEs and $-$0.038 in large firms. The Sieve estimator leaves the SME gap little changed ($-0.055$) but turns the large-firm gap positive ($+0.027$), reversing its sign. The large-firm result is the notable one: a substantial upward correction that implies negative selection bias, namely that conditional on working at a large firm, women have lower mean unobserved wage components than men. This pattern can be explained with amenity-based sorting. Large Korean firms offer substantial non-wage benefits (parental leave, regular working hours, and job security) that are particularly valued by women. If women sort into large firms partly for these amenities and are willing to accept employment even with relatively low wage draws, while men in large firms are selected primarily on the wage dimension, then the female workforce at large firms will have systematically lower unobserved productivity than the male workforce. For SMEs the selection correction is small and its sign varies across estimators: the parametric control function attenuates the gap to $-0.028$, while the Sieve estimator slightly widens it to $-0.055$. This is consistent with a weaker amenity bundle at smaller firms that provides less scope for amenity-wage tradeoffs, leaving the SME gap robustly negative at around $-0.05$.
In multinomial field sorting, the OLS female penalty in non-STEM occupations ($-0.061$) is modestly reduced by the MLogit correction ($-0.055$) but left essentially unchanged by the Sieve estimator ($-0.063$). In STEM, the OLS gap of $-0.014$ shifts to $-0.021$ under Sieve and to $-0.026$ under MLogit CF. The selection corrections on this margin are small, and the STEM gap in particular remains small in absolute terms across all estimators. Both men and women in STEM have passed through similar meritocratic screening (technical credentials, quantitative aptitude, degree requirements) that operates comparably regardless of gender, leaving less scope for gender-differential compositional effects. The single-index Exch-$L1$ estimator falls close to MLogit CF, plausibly because the single-index correction does not fully account for the selection structure.
The sectoral sorting results reflect the wage structure of the Korean public sector. The public-sector OLS female coefficient is close to $0$, consistent with the seniority-based salary schedules that leave little scope for gender-differential pay, and the Sieve estimator barely moves it (to $-0.012$). Women make up 57.9% of public-sector workers but only 41.4% of private-sector workers, confirming strong gender-differential sorting on this margin even though it translates into little within-sector wage gap. In the private sector, the OLS estimate of $-0.058$ is essentially unchanged under Sieve ($-0.057$), whereas the MLogit and Exch-$L1$ single-index corrections move it upward to $-0.047$. The remaining gender gap is the largest across all three architectures, consistent with the greater scope for discretionary wage-setting in the private sector. The divergence between the unrestricted Sieve estimate and the single-index MLogit and Exch-$L1$ estimates indicates that the single-index restriction those estimators impose is violated here.\footnote{Table (ref) in Appendix (ref) reports joint Wald tests on the second-stage control-function sieve basis terms, confirming that the single-index restrictions imposed by MLogit and Exch-$L1$ are empirically binding.}
For log monthly wage regression, the estimated gender gaps are systematically 3--4 log points larger than the hourly wage gaps. The two are linked by the log identity $\log w^{\text{hourly}} = \log w^{\text{monthly}} - \log h^{\text{total}}$, which implies that the female coefficient in the hourly wage specification differs from that in the monthly wage specification by approximately the gender gap in log total hours. Women in the wage sample work 44.3 hours per week vs.\ 47.1 for men ($\log$-ratio $\approx -0.06$), and the bulk of this gap comes from overtime: men report 5.0 overtime hours per week against women's 3.3, while regular hours are similar (42.1 vs.\ 41.0). Dividing monthly pay by total hours mechanically attributes the hours difference to women's hourly rate and shrinks the gap by roughly the log-hours ratio. Despite the level difference between the two wage measures, the selection-correction patterns are qualitatively the same: the large-firm Sieve correction is strongly upward in both specifications, turning the hourly gap positive ($-0.038 \to +0.027$) and nearly eliminating the monthly gap ($-0.074 \to -0.007$).
The overlap condition (Assumption (ref)) is verified in Table (ref) in Appendix (ref). Adequate overlap is confirmed and no trimming is required. The first-stage estimates exhibit strong nonlinearity in the continuous covariates. Table (ref) in Appendix (ref) reports that the marginal sieve terms for age and GPA are strongly nonlinear and the parental-income term is also nonlinear (effective degrees of freedom above one), with several pairwise tensor interactions entering nonlinearly as well, providing sufficient identifying variation in the absence of exclusion restrictions. Two robustness checks, reported in Appendix (ref), confirm the main findings. First, augmenting the first-stage with an additional quasi-continuous covariate (semesters completed) yields similar results (Table (ref)). Second, estimating all three architectures separately for each year produces estimates that are qualitatively consistent with the pooled results (Table (ref)).
The proposed framework permits a decomposition of the raw gender gap into a structural within-category gap, a within-category covariate-composition term, and a between-category sorting term. For gender $g$ in category $k$, let $\bar w_k^g$ denote the mean log wage, $s_k^g$ the share, and $\hat\beta_k$ the estimated female coefficient in the category-$k$ wage regression (so $-\hat\beta_k$ is the regression-adjusted within-category gap). The male-female mean log-wage difference can be written as $$\bar w^M - \bar w^F = \underbrace{\sum_k s_k^M (-\hat\beta_k)}_{\text{structural within}} + \underbrace{\sum_k s_k^M\bigl[(\bar w_k^M - \bar w_k^F) + \hat\beta_k\bigr]}_{\text{covariate composition}} + \underbrace{\sum_k \bar w_k^F (s_k^M - s_k^F)}_{\text{between-category sorting}}.$$ The structural within term is the regression-adjusted gender differential, weighted by male sorting shares; the covariate-composition term collects the part of the raw within-category gap explained by gender differences in covariates (tenure, major, and the like); and the sorting term is the part of the gap generated by women being concentrated in lower-paying categories, valued at female within-category mean wages. Table (ref) reports the decomposition for each architecture. The raw gender gap is approximately 12 log points. The structural within-category gap accounts for 36--41% of it, covariate composition for 52--63%, and pure between-category sorting for only 13% (firm size), 5% (occupation field), and $-4\%$ (sector). Women do sort disproportionately into lower-paying categories (67% at SMEs versus 55% of men, 76% in non-STEM occupations versus 68%, and 27% in the public sector versus 16%), but because the unconditional wage premia across these categories are modest, the mechanical contribution of that sorting to the aggregate gap is small. Most of the raw gap is a within-category phenomenon.
The bottom panel of Table (ref) isolates the effect of the selection correction on the structural within-category gap. Replacing the OLS coefficient with the selection-corrected Sieve coefficient changes this component appreciably only in the ordered (firm-size) architecture, where it falls from $0.043$ to $0.018$ (a selection component of $0.025$). This reflects the offsetting category-level corrections documented above: the corrected female coefficient is negative at SMEs ($-0.055$) but positive at large firms ($+0.027$), and the two nearly cancel in the male-share-weighted average. In the occupation and sector architectures the selection correction barely moves the structural within-category gap (selection components of $-0.003$ and $0.000$), consistent with the small corrections reported in Table (ref).
The between-category sorting term has a direct counterfactual reading: it is how much the gap would change if women sorted like men, holding female within-category mean wages fixed. Equalizing firm-size sorting would reduce the gap by 1.5 log points (13% of the raw gap), and equalizing occupation-field sorting by only 0.7 log points (5%). Equalizing sector sorting would slightly increase the gap (by 0.5 log points), because the public sector pays less than the private sector on average but exhibits a smaller within-sector gender differential, so reducing women's overrepresentation in public pushes them into a sector with a larger penalty. The modest size of these counterfactuals underscores that cross-category sorting, while real, accounts for a small share of the entry-level gap; the larger pieces are the structural within-category penalty and gender differences in covariates within categories.
To examine whether the selection patterns are stable over time, I re-estimate all specifications on three subperiods: 2008--2011, 2012--2015, and 2016--2019. This approach allows both the structural wage parameters and the selection correction to vary freely across periods. The GOMS is a repeated cross-section, hence the `dynamics' documented here reflect cohort-level changes in the entry-level gender wage gap. Changes across subperiods may be driven by compositional shifts in the graduating population, evolving labor market institutions, or macroeconomic conditions affecting cohort-specific labor demand. These sources cannot be disentangled in the present data. Table (ref) reports the estimated female coefficient from the Sieve estimator alongside the OLS baseline for each subperiod.
Three patterns emerge from the hourly wage results. First, the OLS gap narrows steadily over time in most categories: the SME gap from $-0.062$ to $-0.028$, the large firm gap from $-0.064$ to $-0.031$, the STEM gap from $-0.049$ to $+0.012$, and the private sector gap from $-0.078$ to $-0.045$. In non-STEM jobs and the public sector, the decline is more modest. Second, the selection-corrected gaps confirm the main findings from Table (ref). In the ordered model, the large-firm Sieve estimates are far less negative than OLS, turning positive from 2012--2015 onward, while the SME Sieve estimates track OLS closely. In the multinomial models the Sieve corrections are more modest in every period. Third, the corrected gap narrows or is stable over time. In large firms and STEM jobs, the Sieve estimate becomes positive in recent years, while the SME gap narrows toward zero ($-0.018$ by 2016--2019). The corrected gap is near zero throughout in the public sector, whereas the gap persists in non-STEM jobs and the private sector. The monthly wage results are broadly consistent with hourly wage patterns.
This paper establishes semiparametric identification of multilayered selection models without exclusion restrictions. The theoretical contribution provides a unified framework for correcting selection bias in settings where individuals first decide whether to participate and then sort into one of several categories, either vertically or horizontally. The key insight is that, when enough continuous covariates are available, nonlinearity in the selection structure generates sufficient variation to separate the structural outcome parameters from the selection bias. The empirical application shows that selection into firm size in particular materially reshapes the measured gap: the selection correction reverses the sign of the large-firm gender gap, while the corrections on the occupation and sector margins are more modest. A decomposition shows that most of the raw entry-level gap in Korea is a within-category phenomenon, a structural female penalty plus gender differences in covariates, with pure cross-category sorting contributing only modestly. In the most recent cohorts the corrected gap is closed or reversed in large firms and STEM jobs and near zero in the public sector, though it persists in non-STEM occupations and the private sector, suggesting that gender-equity policy should direct attention to the lifecycle dynamics that emerge after market entry.
Several extensions of the framework are worth pursuing. The current analysis treats the selection architecture as given; developing formal specification tests to discriminate between ordered and multinomial selection would strengthen empirical practice. The exchangeability assumption for multinomial selection, while considerably weaker than IIA, may still be restrictive in settings with strongly asymmetric substitution patterns; extending the framework to accommodate richer correlation structures while maintaining tractability is an open challenge. Finally, augmenting the bounds approach for multilayered selection (as in kroft2024lee) with structural restrictions could provide informative bounds in settings where the point identification conditions are not satisfied.