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.
103,226 characters · 9 sections · 101 citation commands
Semi-Nonparametric Models of Multidimensional Matching: an Optimal Transport Approach
In two-sided markets, agents form optimal matches based on their preferences and characteristics, generating a joint surplus sharable between partners according to their relative bargaining power. Matching models are widely used to analyze these dynamics, as seen in the labor market (workers matching with jobs) and the marriage market (spouses matching with each other). However, existing models often fail to capture real-world patterns due to restrictive assumptions. They often assume that agents match exclusively based on a single attribute or a scalar index aggregating multiple characteristics.\footnote{Assortative spousal matching on income, wages, education, risk aversion, and preference for childbearing are investigated by becker1991treatise, grossbard1993theory, pencavel1998assortative, choo2006marries, chiappori2016matching, legros2007beauty and chiappori2008birth among many others. becker1973 and chiappori2012fatter investigate spousal matching that hinges on “ability indices”.} These assumptions are not innocuous. In the labor market, for instance, workers develop highly specialized skills – cognitive for mathematicians, and manual for gymnasts. The single index model fails to capture this specialization. To address these limitations, we need matching models that can directly accommodate multidimensional heterogeneity.\footnote{Studies such as willis1979education and papageorgiou2014learning favor the multidimensional setup over a single index model in the labor market context. Spousal choices are also based on a variety of attributes as shown in becker1991treatise, weiss1997match, qian1998changes, silventoinen2003assortative, hitsch2010matching, and oreffice2010anthropometry.}
A seminal paper, lindenlaub2017, proposes a parametric model for worker-job matching. Her model imposes joint normality of characteristics and does not allow three or more attributes. While the normality assumption enables tractable closed-form solutions for equilibrium assignment and wage functions, it may lead to misleading implications if the true distributions deviate from Gaussian. Moreover, applying the model requires data to conform to normality. Lindenlaub transformed data to standard normal and employed a Gaussian copula to introduce dependence. However, this transformation can distort the underlying relationships between attributes and their post-transformation joint distributions may remain non-normal. Hence, the estimated assignment mechanism may not accurately reflect the true matching process. Another challenge is that the theoretical model predicts deterministic matching patterns, whereas real-world matching processes involve randomness. Lindenlaub introduced error terms to bridge this gap, relying on another strong assumption: the errors are normally distributed and uncorrelated, which, if violated, makes parameters estimated by maximum likelihood unreliable.
We overcome these challenges by developing empirically tractable semi-nonparametric multidimensional matching models. Our theoretical contributions are four-fold. First, we generalize lindenlaub2017’s model by accommodating any arbitrary distributions of attributes. Our general framework is not limited to two-dimensional cases. Second, we introduce an optimal transport approach villani2003,villani2008, dephilippis2014 to derive unique solutions for the equilibrium assignment and wage functions. Applying optimal transport theory provides an elegant solution to our problem.\footnote{Applications of optimal transport have been proven very successful in multiple fields of economics (e.g., ekeland2010notes, chiappori2010hedonic, chiong2016duality, lindenlaub2017, galichon2022cupid, and many more). galichon2017survey provides a comprehensive survey of the literature.} Third, we relax the Gaussian error assumption and allow for arbitrary correlation between errors. Lastly, we derive conditional moments from which production technology parameters, and the equilibrium assignment and wage functions are computed via efficient sieve estimators. While we focus on worker-job matching, our method can be more generally applied to other matching problems such as couple matching in the marriage market.
In the worker-job matching context, each worker possesses distinct skills, and each job requires specific skills to produce output according to production technology. The social planner's problem is to optimally assign workers to jobs to maximize total output in the economy. This problem can be formulated as an optimal transport problem, which yields the unique equilibrium wage and matching functions, accommodating multidimensional heterogeneity with arbitrary distributions. As optimal transport-based matching models predict deterministic matching patterns, following lindenlaub2017, we introduce error terms into the equilibrium assignment and wage functions to maintain the empirical model consistent with optimal transport theory. We estimate production technology, equilibrium assignment, and wage using the semiparametric M-estimation techniques proposed by ai2003 and chen2007. To our knowledge, this paper is the first in the literature to introduce semiparametric M-estimators to multidimensional matching models. Depending on assumptions on error terms, the model is estimated by sieve maximum likelihood (SML), least squares (SLS), or generalized least squares (SGLS). These estimators are efficient, asymptotically normal, and easy to implement. Our estimators perform very well in extensive simulation experiments for a wide class of data generating processes.
We apply our models to lindenlaub2017's matched worker-job data from the U.S. and estimate production technology, equilibrium assignment, and wage functions to investigate the technological shift and its effects on wage inequality between 1990 and 2010. We first estimate the models using the Gaussian transformed data. Our results show a larger technological progress favoring cognitive skills than Lindenlaub's estimates. Furthermore, more flexibility introduced in our models provides a much greater explanation power for the evolution of wage inequality, particularly the `wage polarization' phenomenon featuring stronger wage growth in the bottom and upper tails of the wage distribution relative to the median. The Gaussian model fails to predict wage polarization because, on top of misspecification bias, it restricts the equilibrium wage function to a quadratic form. We also conduct the same analysis on the original data that sharply differ from Gaussian. A virtue of our methods is that they can be directly applied to data without any transformation. The results from the original data show an even more significant technological shift in favor of cognitive abilities than those from the transformed data.
This paper is organized as follows. The remainder of this section discusses the related literature. Section (ref) proposes the optimal transport approach for multidimensional matching. Section (ref) proposes the empirical matching models and establishes the identification. Section (ref) presents the sieve estimators. Section (ref) derives the asymptotic properties of our sieve GLS estimator. Section (ref) conducts simulation experiments. Section (ref) revisits lindenlaub2017's empirical analysis. Section (ref) concludes. Technical proofs and additional theoretical details are provided in the appendix.
choo2006marries (CS henceforth) introduces an empirical transferable utility (TU) model considering discrete characteristics and multidimensional unobserved heterogeneity.\footnote{See galichon2019 for the imperfectly transferable utility model with unobserved heterogeneity.} Their discrete choice framework assumes that unobserved heterogeneity follows the extreme value type I distribution, under which the systemic match surplus is identified by logit formulae. dupuy2014 extend this framework to continuous types. galichon2022cupid allow for non-logit parametric distributions of unobserved heterogeneity and gualdani2023partial show partial identification of the systemic surplus under nonparametric assumptions. These papers rely on the separability assumption that unobserved heterogeneity does not have interactions in generating the surplus.
Our paper takes a different approach closely related to lindenlaub2017 and bojilov2016. Unlike the CS framework, both papers focus on models where agents form matches given their multidimensional continuous attributes that are assumed to be joint normally distributed.\footnote{Alternatively, lise2020 consider a search-theoretic model in which workers are matched to firms in a dynamic setup.} We further extend their models by dispensing with distributional assumptions, thereby offering a more flexible and robust framework for multidimensional matching. We propose efficient econometric procedures that jointly estimate both finite-dimensional parameters and infinite-dimensional functions using conditional moments implied by the model equilibrium, leveraging on the huge literature on the sieve M-estimation.\footnote{shen1997 establishes asymptotic properties of smooth functionals of sieve MLE. newey2003, ai2003, ai2007, and blundell2007 propose efficient sieve IV and sieve minimum distance (SMD) estimators. chen2009 further show that the SMD estimator under proper penalization is consistent and efficient when residuals are potentially nonsmooth. chen2007 provides an extensive overview of sieve estimation of semi-nonparametric models.} In particular, we employ the sieve GLS estimator proposed in chen2007 for our most flexible model specification. This estimator is efficient and computationally simpler than the SMD estimator.
We establish the convergence, efficiency, and asymptotic normality of our sieve estimators relying on the smoothness of unknown nonparametric components (optimal transport maps). This smoothness condition can be verified by applying the results in the mathematical literature on optimal transport maps. caffarelli1992CPAM,caffarelli1992JAMS,caffarelli1996 show the smoothness of transport maps when the distributions of characteristics on both sides are compactly supported. cordero-erausquin2019 further extend the earlier result to the cases where the distributions may have unbounded supports. The degree of the smoothness of a transport map depends on how smooth the densities are.
We consider an environment where every worker with a bundle of skills sorts into a job demanding specific combinations of those skills. Let $\mathcal{X}\subset\mathbb{R}^{d}$ and $\mathcal{Y}\subset\mathbb{R}^{d}$ be spaces of worker and job characteristics endowed with probability measures $P$ and $Q$ respectively. Workers and jobs are described by the corresponding vectors of characteristics $x\in\mathcal{X}$ and $y\in\mathcal{Y}$. Every matched pair produces a single homogeneous good measured in money value according to production technology, $s\left(x,y\right).$ They share the produced quantity through a negotiation process that allows both parties to exploit mutually beneficial outcomes, resulting in wages and profits becoming endogenous at equilibrium. The model shares similarities with frameworks used in family economics,\footnote{For recent reference books on theoretical matching with transferable utility, see browning2014 and chiappori2017.} where matches between individuals create a surplus that involves unobservable utility transfers. Unlike the marriage market, transfers between workers and jobs are observed in the labor market through wages and profits.
At the individual level, workers and jobs have many potential partners. Their choice to form a match depends on the entire set of opportunities to maximize their wages and profits, as stated below:
$w(x)$ can be interpreted as the price that the firm offering job $y$ must pay to match with worker $x$. This price is specific to the worker but not to the job. Similarly, $v(y)$ can be interpreted as the price of job $y$. The second equation in (ref) represents the firm's profit maximization problem given the wage schedule $w(x).$ We analyze such choices and the sharing of surplus from matching, within a market framework that relies on the equilibrium concept of stability. A matching is stable if no individual worker or job, nor any pair of them, prefers to deviate. Formally, if $w\left(x\right)+v\left(y\right)<s\left(x,y\right)$ for some $\left(x,y\right)\in\mathcal{X}\times\mathcal{Y}$, the matching is not stable.
As a part of the stable matching, $w^{*}:\mathcal{X}\to\mathbb{R}$ and $v^{*}:\mathcal{Y}\to\mathbb{R}$ are the solution to the following cost-minimization problem subject to stability constraints:
Here, $\mathcal{W}$ and $\mathcal{V}$ are function spaces that contain all integrable functions with respect to $P$ and $Q$, respectively. A solution to (ref), $\left(w^{*},v^{*}\right)$, also satisfies (ref), meaning that we can interpret $w^{*}(x)$ as the equilibrium wage function for worker $x$ and $v^{*}(y)$ as the equilibrium profit function for firm $y$ in terms of Walrasian equilibrium galichon2017survey. Using the expression of $v(Y)$, we reformulate the optimization problem as:
If $w^{*}\left(x\right)$ is a solution to (ref), then $w^{*}\left(x\right)+c$ is also a solution for any constant $c$. Thus, appropriate normalization is required to obtain a unique solution. We may impose a location constraint, $w\left(x_{0}\right)=0$, for some $x_{0}\in\mathcal{X}$ or a zero integration constraint, $\int_{\mathcal{X}}w\left(X\right)dX=0$.
The total masses of workers and jobs are normalized to one with distributions $P$ and $Q$, respectively. We define a matching as a probability measure $\pi$ on $\mathcal{X}\times\mathcal{Y}$. If worker $x$ is matched with job $y$, $\pi\left(x,y\right)>0$ and the stability constraint holds with equality. Furthermore, if we sum up $\pi\left(x,y\right)$ for all $y$, it should be the total mass of worker $x$. Technically, $\pi$ should satisfy the following feasibility constraints:
The above optimization problem is a linear programming problem in $w$ and $v$ since the minimand and constraints are linear in these two functions. We apply the duality to this problem.\footnote{For the duality, it is assumed that (i) $s:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}\cup\left\{-\infty\right\}$ is an upper-semicontinuous function, and (ii) there are two lower semicontinuous functions $a\in\mathcal{W}$ and $b\in\mathcal{V}$ such that $s\left(x,y\right)\leq a\left(x\right)+b\left(y\right)$ for all $\left(x,y\right)\in\mathcal{X}\times\mathcal{Y}$. See Theorem 5.10 in villani2008 or Theorem 1 in chiappori2010hedonic for details.} Indeed, (ref) is the dual problem of the well-known Monge-Kantorovich optimal transport problem:
where $\mathcal{M}\left(P,Q\right)$ is the set of all probability measures on $\mathcal{X}\times\mathcal{Y}$ satisfying feasibility constraints (ref). We can see $w$ and $v$ as dual variables for feasibility constraints (ref). Thus, the equilibrium wage, $w^{*}\left(x\right)$, should clear the market for worker $x$, and every worker $x$ is employed with $w^{*}\left(x\right)$.
The primal optimal transport problem (ref) can be understood as a social planner's problem whose solution, $\pi^{*}$, associates each $x$ to $y$ with a measurable function $T^{*}$ such that $y=T^{*}(x)$. If there exists such a function $T^{*}$, then we can reformulate (ref) using a deterministic matching function, $T:\mathcal{X}\to\mathcal{Y}$ as follows:
This is called the Monge problem proposed by monge1781histoire. The solution to this problem, $T^{*}$, that maximizes the average overall surplus is called an optimal transport map. The primal problem (ref) always has a solution that is not necessarily deterministic, while the Monge problem monge1781histoire can be ill-posed, meaning that there exists no solution $T^{*}$ satisfying the constraint $T^{*}(P)=Q$ in some cases.\footnote{For example, let $\mathcal{X}=\{0\}\times[-1, 1]$ and $\mathcal{Y}=\{-1, 1\}\times[-1,1]$ with uniform marginal distributions $\mathcal{U}(\mathcal{X})$ and $\mathcal{U}(\mathcal{Y}).$ Given the surplus function $s(x,y)=x'y$, there exists a unique optimal coupling that assigns one half of the mass at $(0,c)$ matches with $(-1,c)$ and the other half with $(1,c)$ for all $c\in[-1, 1].$}
Even when a deterministic solution exists, identifying the solution $\pi^*$ in (ref) is computationally very challenging except in a few cases where analytically tractable solutions exist peyre2019computational. In contrast, the dual problem (ref) involves optimizing over measurable and integrable functions. Under the following conditions, we can establish the unique existence and differentiability of the solution to the dual problem, say $w^{*}:\mathcal{X}\to\mathbb{R}$, which is closely related to $T^{*}$.
Assumption (ref)(ii) is the twist condition, equivalent to the injectivity of $\nabla_{x}s\left(x,y\right)$ for each fixed $x$. The classical Spence-Mirrlees condition can be viewed as a twist condition in a one-dimensional context (see carlier2003duality for more discussions). Assumption (ref) is a smoothness condition on the conjugate function. There are several ways to ensure this assumption, including some restrictions on the supports of $x$ and $y$, the function $s$, or the absolute continuity of $P$ (see Theorem 10.28 and the following remarks in villani2008). Under the above assumptions, the following proposition holds.
This proposition is a direct application of Theorem 10.28 in villani2008. Given the surplus function that satisfies Assumption (ref), one can pin down the unique equilibrium wage function and hence the matching function. For instance, for the surplus function $s(x,y) = x'y,$ the matching function $T^{*}$ is obtained by the gradient of $w^{*}$ i.e. $T^{*}\left(x\right) =\nabla w^{*}\left(x\right)$.
From now on, we specify the surplus function in the form of bi-linear technology following lindenlaub2017:
where $A$ is a $d\times d$ matrix and $b$ is a $d\times1$ vector. Here the elements of $A$ represent worker-job complementabilities and substitutabilities. The diagonal elements capture within-task complementarities and the off-diagonal elements indicate between-task complementarities. $x'b$ represents non-interaction skill terms. Matching assortativity depends crucially on the properties of the surplus function. To fix ideas, as the simplest example of Proposition (ref), consider becker1973's spousal matching model where men and women are endowed with “ability indices”, $x\in\mathcal{X}\subset\mathbb{R}$ and $y\in\mathcal{Y}\subset\mathbb{R}$, respectively. If $\partial^{2}s\left(x,y\right)/\partial x\partial y\geq0$, then the stable matching function $T^{*}:\mathcal{X}\to\mathcal{Y}$ is defined by $T^{*}\left(x\right)=F_{y}^{-1}\left(F_{x}\left(x\right)\right)$ where $F_{x}$ and $F_{y}$ are the cumulative distribution functions of $x$ and $y$. Furthermore, $w^{*}\left(x\right) =\int_{\min\left(\mathcal{X}\right)}^{x}{\nabla_{x}s\left(x,y\right)}|_{y=T^{*}\left(x\right)}dx+c$ is unique up to a constant $c$. $T^{*}$ having this property is defined as positive assortative matching (PAM) in the sense that high-type males match high-type females i.e., the matching function is strictly increasing in $x$. Negative assortative matching (NAM) is the opposite.
In our specification, the properties of $A$ are pivotal to the assortativity of the equilibrium assignment. Proposition 2 of lindenlaub2017 implies that the equilibrium matching function $T^{*}$ satisfies PAM if $A$ is a diagonal matrix with all positive principal minors. The assignment is unaffected by non-interaction terms because $\mathbb{E}_{\pi}\left[X'AY+X'b\right]=\mathbb{E}_{\pi}\left[X'AY\right]+\mathbb{E}_{P}\left[X'b\right]$ and the latter does not depend on the choice of $\pi$. The original problem (ref) and its dual problem (ref) with $s\left(x,y\right)$ can be rewritten in terms of $s^{o}\left(x,y\right)=x'Ay$ as follows
Note that the solution $w^{*}$ to the problem with original $s$ is obtained by $w^*\left(x\right)=w^{o*}\left(x\right)+x'b+c$ with any constant $c$ where $w^{o*}$ is the solution to the problem with $s^{o}$.
Now we further impose conditions on two probability measures, $P$ and $Q$ as well as $A$.
Assumptions (ref)--(ref) serve as primitive conditions to satisfy Assumptions (ref)--(ref), given our production technology (ref). The following statement derives the equilibrium assignment and wage in terms of the solution to the dual Monge-Kantorovich problem (ref).
We can interpret this problem as assigning from $\mathcal{X}$ to $A\mathcal{Y}:=\left\{AY:Y\in\mathcal{Y}\right\}$. Assumption (ref) guarantees the existence of the convex solution to the dual problem (ref).\footnote{Proposition (ref) presents the unique transformation of $x$ to $y$ without requiring this assumption. This proposition implicitly requires conditions in footnote 7 for duality, which are all satisfied under Assumption (ref)(i).} $w^{o*}\left(x\right)$ (and so $w^{*}\left(x\right)$) implicitly depends on $A$.
We now consider the case where attributes of firms and workers are two-dimensional ($d=2$), tailoring the model to our data. Every worker is endowed with a bundle of cognitive and manual skills, $x=\left(x_{C},x_{M}\right)\in\mathcal{X}\subset\mathbb{R}^{2}$. In turn, each firm is endowed with both cognitive and manual skill demands, $y=\left(y_{C},y_{M}\right)\in\mathcal{Y}\subset\mathbb{R}^{2}$. $y_{C}$ ($y_{M}$) corresponds to the productivity or skill requirement of cognitive task $C$ (manual task $M$). With $A=\left(\left(\alpha_{CC},\alpha_{MC}\right)',\left(\alpha_{CM},\alpha_{MM}\right)'\right)$ and $b=\left(\beta_{C},\beta_{M}\right)'$, define $\delta:=\frac{\alpha_{MM}}{\alpha_{CC}}$ that represents the relative level of complementarities across cognitive and manual tasks. When both complementarity parameters are positive, the value of $\delta$ smaller than 1 indicates whether worker-firm complementary in the cognitive task is stronger than in the manual task.
When $x$ and $y$ follow standard joint Gaussian distributions, the nonparametric component of the wage function, $w^{o*}(x)$, becomes a quadratic function of $x$. lindenlaub2017 derives $T^{*}$ and $w^{*}$ in closed-form assuming joint normality of $x$ and $y$ and estimates the production technology to investigate how the technology in the U.S. has evolved. In practice, however, $x$ and $y$ are non-normal as natural skills tend to have skewed distributions. To align the data with the model, lindenlaub2017 converts each element of $x$ and $y$ into a standard normal variable using inverse transform. Their dependence is then modeled using a Gaussian copula. Figure (ref) illustrates how she derives the equilibrium assignment and wage function from transformed data. If the transformed data, $\tilde{x}$ and $\tilde{y}$, follow the bivariate normal distribution, this transformation provides a way of studying dependence independent of the marginals by removing the marginal characteristics. However, the joint distribution of two Gaussian random variables is not in general normal, and hence, the model based on the closed form expression for $T^{*}$ and $w^{*}$ can be misspecified. To avoid such a potential misspecification, we allow $x$ and $y$ to have any arbitrary distributions in our model.
The actual impact of technological shifts on wage distribution may differ from the prediction based on the Gaussian model. lindenlaub2017 shows that (i) wage distributions are positively skewed for any pairs of $\alpha_{CC}$ and $\alpha_{MM}$, (ii) the variance of wage distribution increases as cognitive or manual skill complementarities increases, and (iii) wage skewness is minimized when $\alpha_{CC}=\alpha_{MM}$. However, in our simulations using non-normal distributions (detailed in Appendix (ref)), the obtained wage distribution's skewness does not reach its minimum when $\alpha_{CC}=\alpha_{MM}$.
This section describes an empirical model using the theoretical results in the previous section. Let $\left\{\left(w_{i},x_{i}',y_{i}'\right)\right\}_{i=1}^{n}$ represent an independent and identically distributed (i.i.d.) sequence of $n$ matched observations on the worker $i$'s wage $w_{i}$, her bundle of skills $x_{i}$, and the matched job's skill demands $y_{i}$. Optimal transport-based matching models are not directly applicable to empirical analysis because they produce deterministic predictions. To address this, the models are regularized by introducing unobserved heterogeneity, search frictions, or measurement errors.\footnote{Notice that we could avoid a situation in which unobserved heterogeneity affects the assignment by assuming that it involves non-interaction terms only.} Here we introduce measurement error in the equilibrium functions to keep the model in line with optimal transport theory.
lindenlaub2017 also introduces normally distributed classical measurement errors in her equilibrium solutions and estimates the model using maximum likelihood. Based on the joint normality of $x$ and $y$, the closed-form expression for $w^{*}\left(x\right)$ involves the productivity correlation. However, lindenlaub2017 uses a correlation of error contaminated $y=\left(y_{C},y_{M}\right)\in\mathbb{R}^{2}$, which is different from the actual productivity correlation. Her estimation results show that the estimated variances of measurement errors are greater than one, which is not desirable in the setting where $y_{C}$ and $y_{M}$ are assumed standard normal. If measurement errors are introduced in the assignment equation, then the productivity correlation should be reformulated. Furthermore, if measurement errors are neither normally distributed nor homoskedastic, such a reformulation of correlation does not work. Our methods do not rely on the closed-form solution under bivariate normality, thereby free from this problem.
The introduction of measurement errors can be motivated by the construction of skill measures in data. For instance, sanders2014, lindenlaub2017, and lise2020 use the U.S. Department of Labour Occupational Characteristics Database (O*NET) to determine the levels of skills required to perform each categorical task. O*NET data provides rich information (more than 270 descriptors) on skill requirements for a large number of occupations. There could be measurement errors in three possible ways. First, the researchers conventionally classify the descriptors into predetermined skill categories e.g., “cognitive”, “manual”, and “interpersonal”. However, this decision may be far from clear-cut for many descriptors. Second, the descriptors are aggregated within each category using principal component analysis. This procedure produces inevitable measurement errors even if the descriptors are correctly classified. Lastly, there may be unobserved factors not included in O*NET for skill requirements. Hence measurement errors in skill requirements can be interpreted as (additively separable) unobserved heterogeneity.
The empirical model with measurement errors is defined by:
Here, $\varepsilon_{wi}$ is a scalar measurement error in the observed wage, but it could also be understood as the match-specific idiosyncratic shock added to the equilibrium wage. $\varepsilon_{yi}$ is a $d\times1$ vector of measurement errors in the firm's skill demands. This may arise due to search friction or asymmetric information. Unlike lindenlaub2017's approach, we do not have to impose distributional assumptions on measurement errors, which are vulnerable to misspecification. Instead, we impose the following moment conditions assuming the exogeneity of $x_i$:
Let $\theta=\left(\text{vec}\left(A^{-1}\right)',b'\right)'$ denote a vector of unknown finite-dimensional parameters and $\theta\in\Theta$ where $\Theta$ is a compact subset of $\mathbb{R}^{d^{2}+d}$. The normalizing constant is not in our parameters of interest, and henceforth, we refer to $w\left(x\right):=w^{o*}\left(x\right)+c$ as the constant added infinite-dimensional parameter. We denote $z_{i}=\left(w_{i},x_{i}',y_{i}'\right)'$ and $\rho\left(z_{i};\theta,w\right)=\left(\rho_{w}\left(w_{i},x_{i};\theta,w\right),\rho_{y}\left(y_{i},x_{i};\theta,w\right)'\right)'$, where
For each observation $i$, the model (ref) satisfies the moment conditions (ref). This implies that the following conditional moments hold:
at a true parameter $\left(\theta_{0},w_{0}\right)$. Then $\left(\theta_{0},w_{0}\right)$ are identified via the model (ref) by Proposition (ref) and the exogeneity of $x_i$ as well as the following assumption on $\mathcal{Y}$.
Assumptions (ref) and (ref), combined with Proposition (ref), imply the existence of a deterministic equilibrium characterized by a unique convex function $w_{0}$. When there is no non-interaction term with $b_{0} = 0$, it follows that $\nabla w_{0}^{*} = \nabla w_{0}$. The strict convexity of $w_{0}$ further implies that $\mathbb{E}\left[\nabla w_{0}\left(x_{i}\right)\nabla w_{0}\left(x_{i}\right)^{\prime}\right]$ has full rank, thus identifying $A_{0}$. Additionally, Assumption (ref) is sufficient to identify the nonzero vector $b_{0}$, as stated in the following theorem.
We can further identify $w_{0}^{o*}$ and $c_{0}$ separately under the normalization such as $w_{0}^{o*}\left(x_{0}\right)=0$ for some $x_{0}\in\mathcal{X}$ or $\int_\mathcal{X}w_{0}^{o*}(x)dx=0$ when $\mathcal{X}$ is bounded. With the former constraint, $c_{0}$ and $w_{0}^{o*}\left(x\right)$ are identified with $w_{0}^{*}\left(x\right)=w_{0}^{o*}\left(x\right)+x'b$ since $c_{0}=w_{0}^{o*}\left(x_{0}\right)+c_{0}=w_{0}^{*}\left(x_{0}\right)-x_{0}'b$.
The model parameters are identified by the semiparametric conditional moment restrictions (ref). If the function $w$ is parametrically specified, these moment conditions lead to standard GMM estimation. As $w$ is infinite-dimensional in our specification, we approximate it using sieves. The unknown function $w\in\mathcal{W}$ is approximated by $w_{n}\in\mathcal{W}_{n}$ where $\mathcal{W}_{n}$ is an approximating multivariate function space becoming dense in $\mathcal{W}$ as $n\rightarrow\infty$. We generate $\mathcal{W}_{n}$ via tensor-product construction.
From now on, we assume that there are sets of firms and workers with $d=2$ without loss of generality. Every worker is endowed with a bundle of cognitive and manual skills, $x=\left(x_{C},x_{M}\right)$. In turn, each firm is endowed with both cognitive and manual skill demands, $y=\left(y_{C},y_{M}\right)$. $y_{C}$ ($y_{M}$) corresponds to the productivity or skill requirement of cognitive task $C$ (manual task $M$). In our case, {
} where $\left\{p_{j_{C}}\left(x_{C}\right)\right\}_{j_{C}=0}^{k_{Cn}}$ and $\left\{p_{j_{M}}\left(x_{M}\right)\right\}_{j_{M}=0}^{k_{Mn}}$ are known basis functions of $x_{C}$ and $x_{M}$. Tensor-product space is simple to extend with higher dimensions and easy to implement. For our second and third conditional moment restrictions, we approximate $\partial w_{0}\left(x\right)/\partial x_{C}$ and $\partial w_{0}\left(x\right)/\partial x_{M}$ with same parameter values $\left\{\gamma_{j_{C}j_{M}}\right\}$ used to approximate $w_{0}\left(x\right)$ in $\mathcal{W}_{n}$:
We first consider the model (ref) with normally distributed mean-zero measurement errors that are uncorrelated with each other and a diagonal matrix $A={\rm diag}\left(\alpha_{CC},\alpha_{MM}\right)$ which rules out between-task complementarities. Then this model is effectively the lindenlaub2017 model without joint normality of $x_i$ and $y_i$. We can estimate the parameters using sieve maximum likelihood (SML). Assuming $\varepsilon_{i}\sim N\left(0,\Sigma\right)$, we write the log-likelihood function of model (ref) as \[ L^{*}\left(\theta,\Sigma,w\left(\cdot\right)\right) =-\frac{n}{2}\log\det\left(\Sigma\right) +\sum_{i=1}^{n}\log\left|\det\left(\partial\rho_{i} /\partial\left(w_{i},y_{Ci},y_{Mi}\right)\right)\right| -\frac{1}{2}\sum_{i=1}^{n}\rho_{i}'\Sigma^{-1}\rho_{i}, \] where $\rho_{i}=\rho\left(z_{i};\theta,w\right)$. Solving $\partial L^{*}/\partial\Sigma=0$ for $\Sigma$, we get $\Sigma=\frac{1}{n}\sum_{i=1}^{n}\rho_{i}\rho_{i}'$, which yields the concentrated log-likelihood function of our model
The value of $\left(\theta,w\right)$ maximizing (ref) is the sieve nonlinear full information maximum likelihood estimator of $\left(\theta,w\right)$.
The normality assumption on measurement errors has no theoretical or empirical ground. Without any distributional assumptions on measurement errors, we can still estimate $\left(\theta,w\right)$ using several sieve M-estimators. As $\rho\left(z;\theta,w\right)-\rho\left(z;\theta_{0},w_{0}\right)$ does not depend on $y$ under Assumption (ref), we can apply the sieve generalized least squares (GLS) procedure chen2007 that minimizes the following objective function with respect to $(\theta, w)$: \[ \min_{\left(\theta,w\right)}\sum_{i=1}^{n}\rho\left(z_{i};\theta,w\right)^{\prime}\left[\hat{\Sigma}_{0}\left(x_{i}\right)\right]^{-1}\rho\left(z_{i};\theta,w\right), \] where $\hat{\Sigma}_{0}\left(x\right)$ is a consistent estimator of the optimal weighting matrix $\Sigma_{0}\left(x\right):={\rm Var}\left[\left.\rho\left(z_{i};\theta,w\right)\right|x_{i}=x\right].$ In addition, if $A$ is diagonal, we can rewrite the last two moment conditions as $\rho_{C}\left(y_{C},x;\kappa_{C},w\right)=y_{C}-\kappa_{C}\nabla_{C}w\left(x\right)$ and $\rho_{M}\left(y_{M},x;\kappa_{M},w\right)=y_{M}-\kappa_{M}\nabla_{M}w\left(x\right)$, where $\kappa_{C}=\alpha_{CC}^{-1}$ and $\kappa_{M}=\alpha_{MM}^{-1}$. Table (ref) outlines the three-step procedure to compute the SGLS estimator.
The SGLS estimator allows for arbitrary correlation between measurement errors and heteroskedasticity. We can impose homoskedasticity by assuming $\Sigma_{0}\left(x\right)=\Sigma_{0}$ so that the optimal weighting matrix does not vary with $x$. If we further assume that the measurement errors are uncorrelated i.e., $\Sigma_0$ is diagonal, we can use the sieve least squares (SLS) estimator from step 1 of the three-step procedure. We summarize the key differences in assumptions imposed in different estimation procedures in Table (ref).
We implement the sieve estimators using finite-dimensional Bernstein polynomials to construct the approximating space $\mathcal{W}_{n}$ of $\mathcal{W}$ on $\left[0,1\right]^{2}$. The basis functions are $p_{j_{C}}\left(x_{C}\right) =\binom{k_{C_{n}}}{j_{C}}\left(x_{C}\right)^{j_{C}}\left(1-x_{C}\right)^{k_{C_{n}}-j_{C}}$ and $p_{j_{M}}\left(x_{M}\right) =\binom{k_{M_{n}}}{j_{M}}\left(x_{M}\right)^{j_{M}}\left(1-x_{M}\right)^{k_{M_{n}}-j_{M}}$, where $j_{C}=0,1,\ldots,k_{C_{n}}$, $j_{M}=0,1,\ldots,k_{M_{n}}$, and $\binom{k}{j}$ is a binomial coefficient.\footnote{$x$ does not lie in $\left[0,1\right]^{2}$ in many applications. To satisfy the domain restriction for our simulation studies and empirical application, we use the following linear transformation when $\left(x_{C},x_{M}\right)\in\left[\underline{x}_{C},\overline{x}_{C}\right]\times\left[\underline{x}_{M},\overline{x}_{M}\right]$: $p_{j_{C}}\left(x_{C}\right) =\binom{k_{C_{n}}}{j_{C}}x_{1}^{j_{C}}\left(1-x_{1}\right)^{k_{C_{n}}-j_{C}}$ and $p_{j_{M}}\left(x_{M}\right) =\binom{k_{M_{n}}}{j_{M}}x_{2}^{j_{M}}\left(1-x_{2}\right)^{k_{M_{n}}-j_{M}}$, where $x_{1}=(x_{C}-\underline{x}_{C})/(\overline{x}_{C}-\underline{x}_{C})$ and $x_{2}=(x_{M}-\underline{x}_{M})/(\overline{x}_{M}-\underline{x}_{M})$.}
If $\gamma_{j_{C}j_{M}}=w\left(j_{C}/k_{C_{n}},j_{M}/k_{M_{n}}\right)$, the Bernstein polynomial $w_{n}\left(x;\gamma\right)$ converges uniformly to $w(x)$ by the Stone-Weierstrass approximation theorem (see, e.g., lorentz1986). This provides an approach to imposing shape restrictions on the sieve estimator with a linear constraint which can be solved easily.\footnote{compiani2022market uses linear constraints for the function $w$ to impose monotonicity restrictions and a so-called “diagonal dominance” constraint.} Without any constraint, the equilibrium wage function, $w\left(x\right)+x'b$, is unique and convex. To obtain a more stable estimator, without loss of generality, we impose linear constraints on the Bernstein polynomials, which are necessary for the function to be convex. A detailed description of implementing this convexity constraint in the estimation procedures is provided in Appendix (ref).
To understand technological changes in the production function, the parametric components of the model are of primary interest. The SGLS estimator is ideal in this case because for $\theta$ (i) it is easy to use, $\sqrt{n}$-consistent, and asymptotically normal; (ii) it is semiparametrically efficient; and (iii) the asymptotic variance estimator of $\hat{\theta}$ is consistent and easy-to-compute.\footnote{The sieve minimum distance estimator can be considered. However, when $\rho\left(z;\theta,w\right)-\rho\left(z;\theta_{0},w_{0}\right)$ does not depend on $y$, the SGLS estimator is simpler to implement and computationally faster.} We formally derive its asymptotic properties in the following section.
We establish consistency, convergence rate, asymptotic normality, and semiparametric efficiency of our SGLS estimator using results in chen1998, ai2003, and chen2007. Define $\lambda := (\theta, w(\cdot)).$ Let $\hat{\lambda}_n$ and $\lambda_0$ denote our sieve GLS estimator and the true parameter values, respectively. We first show that $\hat{\lambda}_{n}$ converges to $\lambda_{0}$ at a rate faster than $n^{-1/4}$ under a pseudo norm $\lVert\cdot\rVert$. For any $\lambda_{1}=\left(\theta_{1},w_{1}\left(\cdot\right)\right),\lambda_{2}=\left(\theta_{2},w_{2}\left(\cdot\right)\right)\in\Lambda$, $\lVert\cdot\rVert$ is defined as \[ \lVert\lambda_{1}-\lambda_{2}\rVert^{2} =\mathbb{E}\left[\left(\frac{d\rho\left(z_{i};\lambda_{0}\right)}{d\lambda} \left[\lambda_{1}-\lambda_{2}\right]\right)' \Sigma\left(x_{i}\right)^{-1} \left(\frac{d\rho\left(z_{i};\lambda_{0}\right)}{d\lambda} \left[\lambda_{1}-\lambda_{2}\right]\right)\right], \] where \[ \frac{d\rho\left(z;\lambda_{0}\right)}{d\lambda} \left[\lambda_{1}-\lambda_{2}\right] =
. \] The pseudo metric is comparatively weaker than the standard sup or $L_{2}$ metric, wherein convergence of $\hat{\lambda}_n$ to $\lambda_0$ based on the standard metric implies convergence of $\hat{\lambda}_n$ using the pseudo metric. ai2003 show that $\hat{\lambda}_n$ converging at a rate faster than $n^{-1/4}$ under the weaker metric $||\cdot||$ suffices to derive the $\sqrt{n}$-asymptotic normality of the parametric component, $\hat{\theta}_n$.
Let $\Lambda=\Theta\times\mathcal{W}$ be equipped with a norm $\lVert\lambda\rVert_{s}=\left|\theta\right|_{e}+\lVert w\rVert_{\infty} +\lVert\nabla_{C}w\rVert_{\infty}+\lVert\nabla_{M}w\rVert_{\infty}$, where $\left|\cdot\right|_{e}$ denotes the Euclidean norm and $\lVert w\rVert_{\infty}=\sup_{x\in\mathcal{X}}\left|w\left(x\right)\right|$ is the supremum norm. We introduce the H\"{o}lder class of functions. Let $\left[m\right]$ be the largest nonnegative integer such that $\left[m\right]<m$. A real-valued function $w$ on $\mathcal{X}$ is said to be in H\"{o}lder space $\Lambda^{m}\left(\mathcal{X}\right)$ if it is $\left[m\right]$ times continuously differentiable on $\mathcal{X}$ and \[ \max_{\ell_{1}+\ell_{2}\leq\left[m\right]}\sup_{x} \left|\frac{\partial^{\ell_{1}+\ell_{2}}w\left(x\right)} {\partial x_{C}^{\ell_{1}}\partial x_{M}^{\ell_{2}}}\right|+ \sup_{m_{1}+m_{2}=\left[m\right]}\sup_{x,x'} \left|\frac{\partial^{\left[m\right]}w\left(x\right)}{\partial x_{C}^{m_{1}}\partial x_{M}^{m_{2}}} -\frac{\partial^{\left[m\right]}w\left(x'\right)}{\partial x_{C}^{m_{1}}\partial x_{M}^{m_{2}}}\right| /\left|x-x'\right|_{e}^{m-\left[m\right]} \] is finite. We provide the following assumptions for convergence.
Assumptions (ref)--(ref) are typical conditions imposed in the estimation of conditional mean functions with the tensor product of finite-dimensional linear sieves. We do not explicitly require identification of $\lambda$ here as Assumptions (ref)--(ref) guarantee it by Theorem (ref). Assumption (ref) quantifies the deterministic approximation error of functions in $\Lambda^{m}\left(\mathcal{X}\right)$ by the linear sieve basis functions. Most papers in the literature require $m>d_{\mathcal{X}}/2$, where $d_{\mathcal{X}}$ is the dimension of $\mathcal{X}$. However, our objective function involves $\nabla_{C}w\left(x\right)$ and $\nabla_{M}w\left(x\right)$, so we need a higher order of $m$. Note that the smoothness of $w$ can be verified using the theory of optimal transport. The degree of the smoothness of the solution function $w$ depends on how smooth the densities of $x$ and $y^{*}$ are.\footnote{Assuming the densities are bounded away from zero and infinity, if the densities of variables $x$ and $y^{*}$ belong to the space $\Lambda^{m-2}$, the function $w_{0}$ is a member of $\Lambda^{m}\left(\mathcal{X}\right)$. For a more comprehensive understanding, refer to caffarelli1992CPAM,caffarelli1992JAMS,caffarelli1996 which covers the case of compactly supported $\mathcal{X}$ and $\mathcal{Y}^{*}.$ cordero-erausquin2019 provides an extended result for distributions with unbounded supports.} The following proposition establishes the convergence rate of $\hat{\lambda}_n.$
We now derive the asymptotic normality of the parametric components of the SGLS estimator, $\hat{\theta}_{n}$. Define $D_{v}\left(x\right):=\left(D_{v_{1}}\left(x\right),D_{v_{2}}\left(x\right),D_{v_{3}}\left(x\right),D_{v_{4}}\left(x\right)\right)$ where
Let $v^{*}=\left(v_{1}^{*},v_{2}^{*},v_{3}^{*},v_{4}^{*}\right)$, where $v_{j}^{*}$ solves
Under Assumptions (ref)--(ref), it is clear to see from Lemma B.1 in ai2003 that $|\hat{\theta}_{n}-\theta_{0}|_{e}=o_{p}\left(n^{-1/4}\right)$, $\lVert\hat{w}_{n}-w_{0}\rVert_{2}=\left(\mathbb{E}\left[\left(\hat{w}_{n}\left(x_{i}\right)-w_{0}\left(x_{i}\right)\right)^{2}\right]\right)^{1/2}=o_{p}\left(n^{-1/3}\right)$, and $\lVert\nabla_{C}\hat{w}_{n}-\nabla_{C}w_{0}\rVert_{2},\lVert\nabla_{M}\hat{w}_{n}-\nabla_{M}w_{0}\rVert_{2}=o_{p}\left(n^{-1/4}\right)$. Now the following theorem provides the asymptotic normality of $\hat{\theta}_n.$
The asymptotic variance $V_{1}^{-1}V_{2}V_{1}^{-1}$ can be consistently estimated (see, Remark 4.2 in chen2007) and the standard errors of $\left(\hat{\alpha}_{CC},\hat{\alpha}_{MM}\right)=\left(1/\hat{\kappa}_{C},1/\hat{\kappa}_{M}\right)$ are obtained by using the delta method. Furthermore, if all conditions of Theorem (ref) are satisfied with $\Sigma\left(x\right)=\Sigma_{0}\left(x\right)$, $\hat{\theta}_{n}$ achieves semiparametric efficiency with a consistent estimator $\hat{\Sigma}_{0}\left(x\right)$ of $\Sigma_{0}\left(x\right)$. The estimation of $\Sigma_{0}\left(x\right)$ is straightforward through series least square estimation, using the initial consistent SLS estimator $(\tilde{\theta}_{n},\tilde{w}_{n})$. To ensure the efficiency of the SGLS estimator, $\hat{\Sigma}_{0}\left(x\right)$ is required to exhibit the following uniform convergence rate.
Let $v_{0}=\left(v_{01},v_{02},v_{03},v_{04}\right)$, where $v_{0j}$ solves (ref) with $\Sigma\left(x\right)$ replaced by $\Sigma_{0}\left(x\right)$. Now the following theorem establishes the semiparametric efficiency of $\hat{\theta}_n$.
This section evaluates the finite sample performances of our sieve estimators using known data-generating processes (DGPs). We first generate Monte Carlo samples from lindenlaub2017's quadratic-Gaussian model. Workers' skill bundle, $x$, and occupations' skill requirements, $y$, follow joint Gaussian distributions: $$
\sim N\left(
,
\right),\quad
\sim N\left(
,
\right),$$ from which $\left\{x_i\right\}_{i=1}^n$ are drawn with sample size $n=3000$. Given the production technology \eqref{eq: bi-linear production function} assuming $A$ is diagonal, the equilibrium assignment $y^*$ is provided by $$ y^*=
=\underbrace{
}_{:=J}
,$$ where $$ J=\frac{1}{\sqrt{1+2\delta(\rho_x\rho_y+\sqrt{1-\rho_y^2}\sqrt{1-\rho_x^2})+\delta^2}}
.$$ Recall that $\delta=\frac{\alpha_{MM}}{\alpha_{CC}}$. We generate $\left\{y_i^*\right\}_{i=1}^n$ using the equilibrium assignment function. Then we generate the equilibrium wage $\left\{w_i^*\right\}_{i=1}^n$ using $$ w_{0}^{*}(x)=\frac{\alpha_{CC}}{2}(J_{11}x_C^2 + 2J_{12}x_C x_M + \delta J_{22}x_M^2) +\beta_C x_C + \beta_M x_M + c.$$ Lastly, we draw measurement errors from mean-zero Gaussian distributions: $$ \varepsilon_w \sim N(0,\sigma_w^2),\quad \varepsilon_C \sim N(0,\sigma_C^2),\quad \varepsilon_M \sim N(0,\sigma_M^2),$$ and add them to $w_i^*$, $y_{Ci}^*$, and $y_{Mi}^*$ respectively to generate the observable data $(w_i, y_i, x_i)_{i=1}^n$ following \eqref{eq:OTmodel}. The true parameter values used in simulations are: $$ (\alpha_{CC},\alpha_{MM},\beta_C,\beta_M,c,\rho_x,\rho_y,\sigma_w,\sigma_C,\sigma_M) =(0.5,0.2,1.7,-0.4,30,-0.4,-0.5,2,1,1),$$ which are close to the ML estimates in lindenlaub2017.
We estimate the production technology parameters using lindenlaub2017's parametric ML estimator and our sieve estimators (SML, SLS, and SGLS) across 1000 Monte Carlo samples. As we mentioned earlier, the ML estimator in lindenlaub2017 suffers from bias because the solution uses measurement error contaminated productivity correlation $\tilde{\rho}_y = corr(y)$ which is different from the true productivity correlation $\rho_y = corr(y^*).$ If the measurement errors in $y$ are negligible e.g. $(\sigma_C, \sigma_M)$ are close to 0, the ML estimator works well for this DGP. However, given the current parameter specification, the measurement errors are substantial so the ML estimator can be highly inconsistent. To address this issue, we define a corrected ML estimator (referred to as `ML$^*$') that uses the corrected productivity correlation: $$ \rho_y=\frac{\tilde{\rho}_y \sqrt{var(y_1)var(y_2)}} {\sqrt{var(y_1)-\sigma_C^2}\sqrt{var(y_2)-\sigma_M^2}}.$$ This correction in turn yields much more precise estimates than the original ML estimator. The sieve estimators do not share this problem.
The box plots in Figure (ref) summarize the distributions of parameter estimates delivered by the 5 estimators we consider. For the linear coefficients $(\beta_C, \beta_M),$ all the estimators equally work well. On the other hand, it is clear that the original ML estimator is heavily biased for the complementarity parameters $(\alpha_{CC}, \alpha_{MM})$ as expected. The other 4 estimators perform very well for $(\alpha_{CC}, \alpha_{MM})$ as the distributions of their parameter estimates are centered around the true parameter values. We document the estimators' bias and root-mean-squared errors (RMSE) in Table (ref). It is surprising that the sieve estimators, while more robust than ML and ML$^*$ estimators, tend to be not less efficient than the parametric ML estimators. The SML estimator's root-mean-squared errors (RMSE) are a tad larger than the ML$^*$ estimator for most parameters. The SLS and SGLS estimators, perform very similarly as the measurement errors are uncorrelated, are slightly less efficient than the SML estimator as expected.
Now we consider DGPs for which the Gaussian model is moderately misspecified. We first generate $\{x_{1i},x_{2i}\}_{i=1}^n$ and $\{y_{1i},y_{2i}\}_{j=1}^n$ separately from the Gumbel copula with the shape parameter values $1.3$ and $1.4$ respectively. Then we transform them into standard normally distributed variables: $$ x_{Ci}=\Phi^{-1}(x_{1i}),\quad x_{Mi}=\Phi^{-1}(1-x_{2i}),\quad y_{Cj}=\Phi^{-1}(y_{1j}),\quad y_{Mj}=\Phi^{-1}(1-y_{2j}),$$ so that $(x_{Ci},x_{Mi})$ and $(y_{Cj},y_{Mj})$ are negatively correlated. Define matrices $x$ and $y$ by $$ x:=
,\quad y:=
.$$ The skill demand and supply bundles are standard normally distributed in this DGP but their joint distributions are not Gaussian. The Gumbel copula exhibits asymmetric tail dependence (the upper tail has stronger dependence than the lower tail), whereas the Gaussian copula has symmetric dependence. However, given the current parameter setup, the Gumbel copula does not drastically differ from the Gaussian copula. The production technology is specified the same as before.
There exists no closed-form solution for the equilibrium assignment in this case. We, therefore, numerically solve the equilibrium matching through linear programming for each Monte Carlo sample. To do so, we first compute the pairwise surplus of each possible match between $x$ and $y$ and construct the surplus matrix $S$ whose $ij$ entry is the surplus generated by worker $i$ and firm $j.$ Let $\mathbb{I}_n$ denote a $n \times n$ identity matrix and $\textbf{1}_n$ be a $n \times 1$ one vector. Let $f$ be a vector generated by flattening $S$ by column. Then the solution ($x^*$) to the following linear programming problem provides the equilibrium assignment:
where $A_1 := \mathbb{I}_n \bigotimes \textbf{1}_n'$, $A_2:= \textbf{1}_n' \bigotimes \mathbb{I}_n,$ and $b:= \textbf{1}_{2n}.$ Reshaping $x^*$ into a $n \times n$ matrix gives the optimal transportation matrix, $T$. Then the optimal matching for $x$ is given by $y^*:=Tx$. The wage $w^*$ is computed by the solution to the dual problem (ref).\footnote{Even with the moderate sample size $n=3000$, the constraint matrix is enormous ($6000 \times 9,000,000).$ Solving the linear program (ref) over many Monte Carlo samples is computationally demanding. We employ Gurobi Optimizer 10.0 to solve it efficiently.}
For measurement errors, we consider two different specifications. In the first case, the errors are independently drawn from a gamma distribution, $\Gamma(a,b)$, with $a = 1$, $b=2$. They are demeaned and scaled to have the same means and variances specified in the Gaussian DGP. Under this specification, the ML$^*$ and SML estimators are further misspecified as measurement errors are non-normal. In the second case, we generate the errors from a joint normal distribution in which the errors are correlated as follows: $$
\sim N\left(
,
\right).$$ The ML$^*$ and SML estimators are still misspecified as measurement errors are correlated. The SLS estimator is consistent but not as efficient as the SGLS estimator as it does not take the correlation structure of errors into account. The observable data $(y_i,x_i,w_i)_{i=1}^n$ are generated by adding the measurement errors to $w_i^*$, $y_{Ci}^*$, and $y_{Mi}^*$ respectively.
The estimation results are provided in Table (ref). In both cases, the ML$^*$ estimator is misspecified for both the distributions of $X,$ $Y,$ and measurement errors so that it performs the worst for all the parameters. It exhibits especially large biases and RMSE for complementarity parameters. Even for the linear productivity parameters, the ML$^*$ estimator shows much larger RMSEs than the sieve estimators. In contrast, all the sieve estimators equally work well for the linear coefficients. The SML estimator is misspecified for the distributions of measurement errors but it produces accurate estimates for the complementary parameters $(\alpha_{CC}, \alpha_{MM})$ in both cases. The SLS and SGLS estimators perform similarly to the SML estimator when the measurement errors are drawn from the Gamma distributions. In the case of correlated errors, the SLS estimator performs similarly to the SML estimator. The SGLS outperforms the other estimators as it takes into account the correlations between measurement errors, resulting in more efficient estimation.
Lastly, we consider a DGP in which lindenlaub2017's model is more severely misspecified. Specifically, we draw $x$ and $y$ from finite Gaussian mixture distributions. Each mixture distribution has two Gaussian components with one-half weight for each. For both $x$ and $y$, the Gaussian components, $K_1$ and $K_2$, are specified as follows. $$ K_1 \sim N\left(
,
\right),\quad K_2 \sim N\left(
,
\right).$$ We set $\rho$ equal to $0.4$ for $x$ and $0.5$ for $y$. The equilibrium assignments and wages are solved via linear programming as before. We also generate the measurement errors from a joint Gaussian mixture distribution which has two components: $$ M_1 \sim N\left(
,
\right),\quad M_2 \sim N\left(
,
\right).$$ In this case, both $(x,y)$ and measurement errors have bi-modal distributions that are far from a normal distribution.
As the margins of $x$ and $y$ are not standard normal, the ML estimator is not directly applicable. Therefore, we use the inverse transform method to convert $x$ and $y$ to standard normal variables for the ML$^*$ estimator. Our sieve estimators can be applied without this transformation so we use untransformed data for the sieve-based estimators. The estimates of technology parameters are reported in Table (ref). Not surprisingly, the ML$^*$ estimator delivers parameter estimates that are very different from the true values. On top of misspecification, the transformation procedure introduces additional bias as the actual assignments are determined on the original data. On the contrary, the sieve estimators still perform extremely well in this case. Estimators relying on fewer assumptions deliver more accurate estimates. The SML estimator produces the least precise estimates among the sieve estimators. The SGLS estimator incorporates the correlation structure among measurement errors so it possesses substantial efficiency gains compared to the SLS estimator.
Our simulation exercises provide evidence that the Gaussian model can be misleading when the model is misspecified. Even with moderate misspecification, the ML estimator does not produce reliable estimates. Furthermore, Gaussian transformation is necessary if the marginal distributions of skill supply and requirements are not standard normal. This transformation may not precisely recover the underlying assignment mechanism between workers and jobs even if the model is correctly specified after transformation. On the contrary, our sieve estimators do not suffer from these problems. Therefore, our estimators can be more generally applicable regardless of underlying distributions of $x$, $y$, and measurement errors with no need for transformation. They also show excellent finite sample performances under correct specifications.
We revisit the dataset constructed by lindenlaub2017 to learn how production technology in the US has evolved. We estimate the production technology parameters in the model using the dataset and sieve estimators. The National Longitudinal Survey of Youth (NLSY) data and U.S. Department of Labor Occupational Characteristics Database (O{*}NET) are used to construct workers' cognitive and manual skills as well as the occupational skill requirements of firms. To assess the effect of technological changes on wage inequality, we compare estimation results based on two cohorts: the first cohort commencing in 1979 (referred to as NLSY79) and the second commencing in 1997 (referred to as NLSY97). Following Lindenlaub's main specification, we focus on employed workers aged 27 to 29 during the years 1990-91 and 2009-10, sourced from the NLSY79 and NLSY97 cohorts, respectively.\footnote{The dataset excludes military samples and oversamples of special demographic/racial groups to give primary focus on the core sample of the NLSY.} The wage, $w$, is defined as the CPI-adjusted hourly rate.
Firms' skill demands, $(y_C,y_M),$ is constructed from the O{*}NET, which contains information on skill requirements for each occupation. sanders2014 classifies occupational skill requirements into cognitive and manual categories and constructs two task scores for over 400 occupations. These scores are employed to obtain skill demands for individuals' currently matched jobs.\footnote{For instance, the occupation `dancer' has a normalized cognitive score ($y_C$) of $0.34$ and a normalized manual score ($y_M$) of 1, indicating the job is highly manual. On the contrary, the highly cognitive job `physicist' has a skill demand bundle of $(y_C=1,y_M=0.11).$ We use the normalized task scores for illustration purposes following lindenlaub2017. The original supports of worker skills and firms' skill demands are provided in Table (ref).} To construct the skill supply bundle $(x_C,x_M)$, survey responses in the NLSY on education and training are used. Given college education, apprenticeships, government training degrees, and occupational training history, workers are qualified for specific occupations. lindenlaub2017 matches individuals to their qualified occupations and obtains the value of $\left(x_{C},x_{M}\right)$ using the normalized skill requirements $\left(y_{C},y_{M}\right)$ of the matched jobs.\footnote{For example, a worker who studied economics at a university is qualified for the `economist' job. Then the worker possesses a normalized skill bundle of $\left(x_{C}=0.615,x_{M}=0.034\right)$, that is required to be an economist.} It is important to note that the workers' skills are independent of their current occupation because the skill bundles are constructed using qualified occupations.
Table (ref) presents summary statistics of workers' skills and firms' skill demands. In 1990/91, workers had higher cognitive skills on average than manual skills, and firms also required more cognitive skills than manual skills. Two decades later, workers had increased cognitive skills and decreased manual skills on average compared to 1990/91, with firms also showing a similar trend. The skill correlation ($\rho_x$) shifted from $-0.40$ to $-0.52$, indicating increased worker specialization. In contrast, the productivity correlation ($\rho_y$) remained stable at $-0.49$. Initially, jobs were more specialized than workers, but skill supply caught up, resulting in slightly greater worker specialization in 2009/10.
lindenlaub2017 transforms each element of $x$ and $y$ into a standard Gaussian variable, and their dependence is modeled using the Gaussian copula. However, while their margins are standard normally distributed, the transformed variables are not guaranteed to be joint normally distributed. Let $\tilde{x}$ and $\tilde{y}$ be Gaussian-transformed $x$ and $y$ respectively. As shown in Figure (ref), the joint distributions of $\tilde{x}$ and $\tilde{y}$ are not normal. Especially, the joint density of $\tilde{y}$ is multi-modal in both periods.
We employ Mardia's test mardia1970 to formally test the joint normality of $\tilde{x}$ and $\tilde{y}$. We first create a $n \times n$ matrix for $\tilde{x}$: \[ C=\left(c_{ij}\right)=x^{*}S^{-1}\left(x^{*}\right)', \] where the $i$-th row of $x^{*}$ is $x_{i}^{*}=\tilde{x}_{i}-\sum_{i=1}^{n}\tilde{x}_{i}/n$, and define multivariate measures of skewness and kurtosis as follows: \[ b_{1}=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}c_{ij}^{3},\quad b_{2}=\frac{1}{n}\sum_{i=1}^{n}c_{ii}^{2}. \] Under bivariate normality, the limiting distribution of $\frac{nb_{1}}{6}$ is a chi-square distribution with $d\left(d+1\right)\left(d+2\right)/6$ degrees of freedom and the limiting distribution of $\frac{\sqrt{n}\left(b_{2}-d\left(d+2\right)\right)}{\sqrt{8d\left(d+2\right)}}$ is the standard normal distribution where $d$ is the dimensionality of $\tilde{x}$. We conduct the same procedure for $\tilde{y}$. Table (ref) shows the test statistics. The test strongly rejects the bivariate normality of $\tilde{x}$ and $\tilde{y}$ in both periods. The normality assumption imposed in the Lindenlaub model is not satisfied even after transforming $x$ and $y$. Therefore, our semiparametric approach is more appropriate in this case.
We estimate the production technology in each period separately. First, we modify lindenlaub2017's MLE procedure to accommodate the cases where the true skill requirements $\tilde{y}^*$ have variances not equal to $1$. As the inverse transform method converts the measurement error contaminated $\tilde{y}$, not $\tilde{y}^*$, it is essentially the case that $var(\tilde{y}^*)\neq1$. We also allow the Gaussian measurement errors to be correlated with each other. We refer to this generalized and corrected ML procedure as `ML$^*$'. Then, we relax the normality of $X$ and $Y$ and estimate the technology parameters using SML. Next, we further relax the normality of the measurement errors and estimate the parameters using SLS. Finally, we allow measurement errors to be correlated with each other and estimate the parameters through SGLS. Sieve estimation is conducted using Bernstein polynomial basis functions of degree $3$, which performs the best in terms of information criteria and model fit. We compare our semiparametric estimates of the technology parameters to Lindenlaub's original estimates and the `ML$^{*}$' estimates.
The estimation results are provided in Table (ref). All the models clearly show a huge shift in the relative importance between manual and cognitive tasks over the two decades. Our results are qualitatively consistent with Lindenlaub's but quantitatively very different. In 1990/91, the estimated complementarity in manual tasks was much larger than in cognitive tasks. The Gaussian model (ML$^*$) indicates that the complementarity in manual tasks is roughly 1.7 times as large as that of cognitive tasks. Dispensing with the normality of skill demand and supply, the ratio becomes larger than 3. When we further generalize the model by removing the Gaussian assumption on measurement errors, the complementarity in cognitive tasks shrinks close to 0 (but significantly larger than 0), whereas that of manual tasks is still close to 1. The estimates of linear productivity coefficients $\beta_C$ and $\beta_M$ are similar across all the specifications as shown in simulations.
This pattern becomes the opposite in the 20 years. All the complementarity estimates for 2009/10 indicate a substantial increase in the complementarity between cognitive worker and job attributes, whereas the complementarity in manual tasks heavily decreased. In the Gaussian model (ML$^*$), the ratio of estimated complementarities in manual tasks to cognitive tasks ($\frac{\alpha_{MM}}{\alpha_{CC}}$) is around $0.43$, which is similar to the estimated value by SLS. However, SML and SGLS deliver much smaller values close to $0.1$. The linear coefficients are similar across specifications. These patterns imply substantial changes in the relative complementarities across tasks because of technological advances. Lindenlaub describes this phenomenon as “task-biased technological change in favor of cognitive tasks”. The cognitive dimension became much more important in labor market sorting. Our semiparametric models suggest that the “task-biased technological change” favoring cognitive tasks in the last two decades may have been much larger than previously found. The increases in $\beta_C$ and $\beta_M$ indicate that both cognitive and manual skill productivity have risen. However, the estimated manual skill productivity in both periods is insignificant in most specifications. Therefore, following Lindenlaub's description, we can conclude that the U.S. economy has also experienced “skill-biased technological change” in favor of cognitive skills.
We now investigate the effect of estimated technological changes on wage inequality. Wage inequality in the U.S. labor market until the late 2000s is well characterized by wage polarization that is defined as stronger wage growth in the bottom and upper tails relative to the median. The red solid line in Figure (ref) plots how wages changed relative to the median wage between 1990/91 and 2009/10 by wage percentile, implying that the U.S. labor market experienced a spike in the upper-tail wage inequality, while the lower-tail inequality declined. This phenomenon is surprisingly well-predicted in our models as shown in Figure (ref) (a). All the semiparametric models predict substantial wage polarization once the estimated parameter values are fed in. On the contrary, the Gaussian model fails to account for wage polarization in both tails. We can also observe that the model fit improves as the model becomes more flexible.
To further explore why the matching model requires greater flexibility to account for wage polarization, we compare the curvature of estimated wage functions across different models in Figures (ref)--(ref). In 1990/91, all models produce almost linear wage functions, indicating a relatively uniform relationship between wages and cognitive skills. However, in 2009/10, our semiparametric models predict a significantly steeper curvature, particularly at high cognitive skill levels, whereas the Gaussian model still generates a more linear wage function. The Gaussian model constrains the wage function to a quadratic form of standard normal variables, limiting its shape to a low-degree polynomial. In contrast, our models do not impose such constraints, allowing for greater flexibility in curvature that better fits the data. Notably, our most flexible model predicts a sharply increasing slope in the wage function for 2009/10, which is relatively flat at low cognitive skill levels and very steep at high skill levels, generating substantial wage polarization.
To understand the driving forces behind wage polarization, we isolate the effects of technological and distributional changes in Figure (ref). We only keep task-biased technological change (shutting down changes in linear productivity coefficients) in panel (b), skill-biased technological change (shutting down changes in complementarity parameters) in (c), and distributional change (shutting down both changes in linear productivity and complementarity parameters) in (d). We find that task-biased technological change explains wage polarization remarkably well. Especially, all three semiparametric models exhibit an excellent fit in the lower tail in panel (b), while the Gaussian model shows only a slight decline in lower tail inequality. In contrast, skill-biased technological change exacerbates wage inequality in the lower tail as shown in panel (c). The distributional change has a negligible impact on wage inequality. In summary, despite their parsimony, our matching models effectively account for wage inequality's evolution over the past 20 years in the U.S., with task-biased technological change being the primary driver of the observed pattern.
Lastly, we estimate our semiparametric models on the original data. Unlike the Gaussian model, our models can be applied directly to the data without any transformation. As we described in Figure (ref), the matching occurs based on original distributions. Hence transforming marginal distributions to standard normal can lead to a different solution from $T(x).$ The marginal distributions of workers' skill supply $x$ and firms' skill demand $y$ in Figure (ref) are skewed and multi-modal, quite different from standard normal. Therefore, it is crucial to investigate the robustness of the Gaussian model using the original data. We report the estimated parameters in Table (ref), using Bernstein polynomial basis functions of degree 4 to accommodate the less well-behaved original data. The estimated parameters reveal similar patterns across our models. Notably, the complementarity in cognitive tasks increased significantly from near 0 in 1990/91 to around 6 in 2009/10, while the complementarity in manual tasks decreased from near 2 to almost 0. Additionally, linear skill productivity improved, with a larger increase in cognitive skill productivity. These findings confirm that the U.S. economy experienced substantial task-biased and skill-biased technological changes favoring cognitive skills over the two decades.
The estimated models on the original data effectively capture the patterns of wage polarization, particularly in the upper tail as in Figure (ref). While the models slightly under-predict wage polarization in the lower tail, they confirm that task-biased technological change was the primary driver of wage polarization. In the absence of skill-biased technological change, the model shows a significant relative wage increase in the lower tail. However, skill-biased technological change had a negative impact on wage inequality, exacerbating lower tail inequality. Distributional change improved upper tail inequality but worsened lower tail inequality.
In summary, the estimated semiparametric models on the original data exhibit similar patterns to those on the Gaussian transformed data, with task- and skill-biased technological changes being more pronounced. Our models demonstrate a remarkable fit for wage polarization, highlighting the substantial changes in production technology in the U.S. over the past two decades. This exercise showcases the versatility and effectiveness of our semiparametric models and sieve-based estimators, which can accommodate any underlying joint distributions of skill supply and demand without requiring data transformation or distributional assumptions. Moreover, our approach builds upon standard sieve-based estimators, which have been proven to achieve semiparametric efficiency and are easy to implement in practice.
Theoretical matching models often face empirical challenges due to discrepancies between model assumptions and real-world data. lindenlaub2017 presents a tractable theoretical model suitable for comparative statics and qualitative analysis of multidimensional matching. However, its empirical application is limited by restrictive distributional assumptions on observed characteristics and measurement errors. We generalize this model by relaxing these key distributional restrictions, enabling our models to accommodate datasets with matched pairs, regardless of the underlying characteristic and error distributions. Our simulation results demonstrate the accuracy of our semi-nonparametric estimators across various data generating processes. Moreover, our flexible models generate significant wage polarization, aligning with U.S. data patterns, whereas the parametric Gaussian model falls short. Our estimated models indicate that task-biased technological progress, favoring cognitive abilities over manual skills, is the primary driver of wage polarization.
Our study opens up promising research avenues. First, incorporating additional dimensions like interpersonal and digital skills into our model holds great potential. By employing advanced techniques like artificial neural networks chen2023efficient to approximate equilibrium functions in high-dimensional spaces, one could enhance the model's accuracy in elucidating intricate matching patterns and their impact on wage inequality. Second, addressing measurement errors in assessing worker skills is crucial. Overcoming this challenge requires innovative econometric approaches that can accurately estimate models amidst multidimensional measurement errors. Finally, our framework's application extends beyond the worker-job matching problem, offering insights into matching problems in diverse contexts like the marriage market.
While our approach offers valuable insights, it is essential to acknowledge its limitations. Frictionless matching models, grounded in optimal transport, assume efficient equilibrium, eliminating concerns of unemployment or skill mismatch. Introducing randomness in assignment through factors like search frictions and unobserved heterogeneity poses a fruitful challenge, especially when extending existing one-dimensional theories e.g., in eeckhout2011, to multidimensional settings.