The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
103,225 characters
Semi-Nonparametric Models of Multidimensional Matching: an Optimal Transport Approach
\begin{frontmatter}
\title{Semi-Nonparametric Models of Multidimensional Matching: an Optimal Transport Approach}
\runtitle{Semi-Nonparametric Multidimensional Matching}
\begin{aug}
\author[id=au1,addressref={add1}]{\fnms{Dongwoo}~\snm{Kim}\ead[label=e1]{[email removed]}}
\author[id=au2,addressref={add2}]{\fnms{Young Jun}~\snm{Lee}\ead[label=e2]{[email removed]}}
\address[id=add1]{
\orgdiv{Department of Economics},
\orgname{Simon Fraser University}}
\address[id=add2]{
\orgname{Korea Institute for International Economic Policy}}
\end{aug}
\support{This paper is based on the third chapter of Lee's doctoral dissertation at University College London. We thank Dennis Kristensen, Krishna Pendakur, Martin Weidner, and Daniel Wilhelm for their helpful suggestions. We also greatly benefited from the seminar participants at UCL, Bocconi University, KIEP, and the KAEA micro virtual seminar. All errors are our own. The authors gratefully acknowledge support from the Social Sciences and Humanities Research Council of Canada under the Insight Grant (435-2024-0322) and PRIN PROJECT 2017 (prot.2017TMFPSH).}
\begin{abstract}
This paper proposes empirically tractable multidimensional matching models, focusing on worker-job matching. We generalize the parametric model proposed by Lindenlaub (2017), which relies on the assumption of joint normality of observed characteristics of workers and jobs. In our paper, we allow unrestricted distributions of characteristics and show identification of the production technology, and equilibrium wage and matching functions using tools from optimal transport theory. Given identification, we propose efficient, consistent, asymptotically normal sieve estimators. We revisit Lindenlaub's empirical application and show that, between 1990 and 2010, the U.S. economy experienced much larger technological progress favoring cognitive abilities than the original findings suggest. Furthermore, our flexible model specifications provide a significantly better fit for patterns in the evolution of wage inequality.
\vspace{1em}\\
Keywords: Multidimensional matching, transferable utility, optimal transport, sieve extremum estimation, technological progress, wage polarization.
\end{abstract}
\end{frontmatter}
\newpage
\section{Introduction}\label{sec: Matchingintro}
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 \cite{becker1991treatise}, \cite{grossbard1993theory}, \cite{pencavel1998assortative}, \cite{choo2006marries}, \cite{chiappori2016matching, legros2007beauty} and \cite{chiappori2008birth} among many others. \cite{becker1973} and \cite{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 \cite{willis1979education} and \cite{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 \cite{becker1991treatise}, \cite{weiss1997match}, \cite{qian1998changes}, \cite{silventoinen2003assortative}, \cite{hitsch2010matching}, and \cite{oreffice2010anthropometry}.}
A seminal paper, \cite{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 \cite{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 \citep{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., \cite{ekeland2010notes}, \cite{chiappori2010hedonic}, \cite{chiong2016duality}, \cite{lindenlaub2017}, \cite{galichon2022cupid}, and many more). \cite{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 \cite{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 \cite{ai2003} and \cite{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 \cite{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 `\textit{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{sec: Matchingmodel} proposes the optimal transport approach for multidimensional matching. Section \ref{sec: MatchingModelID} proposes the empirical matching models and establishes the identification. Section \ref{sec:MatchingEst} presents the sieve estimators. Section \ref{sec: asymptotics} derives the asymptotic properties of our sieve GLS estimator. Section \ref{sec: Simul} conducts simulation experiments. Section \ref{sec: Matchingemp} revisits \citet{lindenlaub2017}'s empirical analysis. Section \ref{sec: conclusion} concludes. Technical proofs and additional theoretical details are provided in the appendix.
\subsection{Related literature}
\citet{choo2006marries} (CS henceforth) introduces an empirical transferable utility (TU) model considering discrete characteristics and multidimensional unobserved heterogeneity.\footnote{See \citet{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. \citet{dupuy2014} extend this framework to continuous types. \citet{galichon2022cupid} allow for non-logit parametric distributions of unobserved heterogeneity and \cite{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 \citet{lindenlaub2017} and \citet{bojilov2016}. Unlike the CS framework, both papers focus on models where agents form matches given their multidimensional \textit{continuous} attributes that are assumed to be joint normally distributed.\footnote{Alternatively, \citet{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{\citet{shen1997} establishes asymptotic properties of smooth functionals of sieve MLE. \citet{newey2003}, \citet{ai2003, ai2007}, and \citet{blundell2007} propose efficient sieve IV and sieve minimum distance (SMD) estimators. \citet{chen2009} further show that the SMD estimator under proper penalization is consistent and efficient when residuals are potentially nonsmooth. \citet{chen2007} provides an extensive overview of sieve estimation of semi-nonparametric models.} In particular, we employ the sieve GLS estimator proposed in \citet{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. \citet{caffarelli1992CPAM,caffarelli1992JAMS,caffarelli1996} show the smoothness of transport maps when the distributions of characteristics on both sides are compactly supported. \citet{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.
\section{Optimal transport approach for multidimensional matching}\label{sec: Matchingmodel}
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 \cite{browning2014} and \cite{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:
\begin{equation}\label{eq: wage and profit}
w\left(x\right)=\sup_{y\in\mathcal{Y}}\left\{s\left(x,y\right)-v\left(y\right)\right\},\quad
v\left(y\right)=\sup_{x\in\mathcal{X}}\left\{s\left(x,y\right)-w\left(x\right)\right\}.
\end{equation}
$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 \eqref{eq: wage and profit} 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:
\begin{equation}\label{eq: dual-MK problem}
\inf_{w\in\mathcal{W}, v\in\mathcal{V}}
\left\{\mathbb{E}_{P}\left[w\left(X\right)\right]
+\mathbb{E}_{Q}\left[v(Y)\right]\right\}, \textit{ s.t. } w\left(x\right)+v\left(y\right)\geq s\left(x,y\right),\
\forall\left(x,y\right)\in\mathcal{X}\times\mathcal{Y}.
\end{equation}
Here, $\mathcal{W}$ and $\mathcal{V}$ are function spaces that contain all integrable functions with respect to $P$ and $Q$, respectively. A solution to \eqref{eq: dual-MK problem}, $\left(w^{*},v^{*}\right)$, also satisfies \eqref{eq: wage and profit}, 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 \citep{galichon2017survey}. Using the expression of $v(Y)$, we reformulate the optimization problem as:
\begin{equation}\label{eq: dual-matching}
\inf_{w\in\mathcal{W}}
\left(\mathbb{E}_{P}\left[w\left(X\right)\right]
+\mathbb{E}_{Q}\left[\sup_{x\in\mathcal{X}}
\left\{ s\left(x,Y\right)-w\left(x\right)\right\} \right]\right).
\end{equation}
If $w^{*}\left(x\right)$ is a solution to \eqref{eq: dual-matching}, 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:
\begin{equation}\label{eq: feasibility}
\int_{\mathcal{Y}}d\pi\left(x,y\right)=P\left(x\right),\ \forall x\in\mathcal{X},\quad
\int_{\mathcal{X}}d\pi\left(x,y\right)=Q\left(y\right),\ \forall y\in\mathcal{Y}.
\end{equation}
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 \citet{villani2008} or Theorem 1 in \cite{chiappori2010hedonic} for details.} Indeed, \eqref{eq: dual-MK problem} is the dual problem of the well-known Monge-Kantorovich optimal transport problem:
\begin{equation}
\sup_{\pi\in\mathcal{M}\left(P,Q\right)}
\mathbb{E}_{\pi}\left[s\left(X,Y\right)\right]
:=\int_{\mathcal{X}\times\mathcal{Y}}
s\left(x,y\right)d\pi\left(x,y\right),\label{eq: Kantorovich-matching}
\end{equation}
where $\mathcal{M}\left(P,Q\right)$ is the set of all probability measures on $\mathcal{X}\times\mathcal{Y}$ satisfying feasibility constraints \eqref{eq: feasibility}. We can see $w$ and $v$ as dual variables for feasibility constraints \eqref{eq: feasibility}. 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 \eqref{eq: Kantorovich-matching} 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 \eqref{eq: Kantorovich-matching} using a deterministic matching function, $T:\mathcal{X}\to\mathcal{Y}$ as follows:
\begin{equation}\label{eq: monge problem}
\max_{T(\cdot)} \mathbb{E}_{P}[s(X,T(X))], \quad {\rm s.t. } \ T(P)=Q.
\end{equation}
This is called the Monge problem proposed by \cite{monge1781histoire}. The solution to this problem, $T^{*}$, that maximizes the average overall surplus is called an optimal transport map. The primal problem \eqref{eq: Kantorovich-matching} always has a solution that is not necessarily deterministic, while the Monge problem \cite{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 \eqref{eq: Kantorovich-matching} is computationally very challenging except in a few cases where analytically tractable solutions exist \citep{peyre2019computational}. In contrast, the dual problem \eqref{eq: dual-MK problem} 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^{*}$.
\begin{assumption}
\label{assu:den1}
(i) $s\left(x,y\right)$ is differentiable as a function of $x$, for all $y$; (ii) For any fixed $x\in\mathcal{X}$ and $y_{1}\neq y_{2}\in\mathcal{Y}$,
$\nabla_{x}s\left(x,y_{1}\right)\neq\nabla_{x}s\left(x,y_{2}\right)$.
\end{assumption}
\begin{assumption}
\label{assu:Phi1} For any $v:\mathcal{Y}\mapsto\mathbb{R}\cup\left\{\pm\infty\right\}$, $w\left(x\right)=\sup_{y\in\mathcal{Y}}\left\{s(x,y)-v(y)\right\}$ is differentiable
almost surely on
\[
\left\{x|\exists y\in\mathcal{Y},\
s\left(x,y\right)-w\left(x\right)\geq s\left(z,y\right)-w\left(z\right),\
\forall z\in\mathcal{X}\right\}.
\]
\end{assumption}
\noindent Assumption \ref{assu:den1}(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 \cite{carlier2003duality} for more discussions). Assumption \ref{assu:Phi1}
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 \citet{villani2008}). Under the above assumptions, the following proposition holds.
\begin{prop}\label{prop:MatchingEq}Let Assumptions \ref{assu:den1}--\ref{assu:Phi1} hold. Then, there exists a unique (up to a constant) equilibrium wage function, $w^{*}\left(x\right)$, solving
the dual problem \eqref{eq: dual-matching}. Furthermore, the function, $T^{*}\left(x\right)$,
satisfying $\nabla w^{*}\left(x\right):=\left.\nabla_{x}s\left(x,y\right)\right|_{y=T^{*}\left(x\right)}$,
is the unique equilibrium assignment for the Monge problem \eqref{eq: monge problem}.
\end{prop}
This proposition is a direct application of Theorem 10.28 in \citet{villani2008}. Given the surplus function that satisfies Assumption \ref{assu:Phi1}, 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 \cite{lindenlaub2017}:
\begin{equation}\label{eq: bi-linear production function}
s\left(x,y\right):=s\left(x,y;A,b\right)
=x'Ay+x'b,
\end{equation}
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{prop:MatchingEq}, consider \citet{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 \citet{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 \eqref{eq: Kantorovich-matching} and its dual problem \eqref{eq: dual-matching} with $s\left(x,y\right)$ can be rewritten in terms of $s^{o}\left(x,y\right)=x'Ay$ as follows
\begin{equation}\label{eq: dual-matching2}
\begin{split}
&\inf_{w\in\mathcal{W}}\left\{\mathbb{E}_{P}\left[w\left(X\right)\right]
+\mathbb{E}_{Q}\left[\sup_{x\in\mathcal{X}}
\left\{s\left(x,Y\right)-w\left(x\right)\right\}\right]\right\}\\
&=\sup_{\pi\in\mathcal{M}\left(P,Q\right)}
\mathbb{E}_{\pi}\left[s\left(X,Y\right)\right]
=\sup_{\pi\in\mathcal{M}\left(P,Q\right)}
\mathbb{E}_{\pi}\left[s^{o}\left(X,Y\right)\right]+\mathbb{E}_{P}\left[X\right]'b\\
&=\inf_{w^{o}\in\mathcal{W}}\left\{\mathbb{E}_{P}\left[w^{o}\left(X\right)\right]
+\mathbb{E}_{Q}\left[\sup_{x\in\mathcal{X}}
\left\{s^{o}\left(x,Y\right)-w^{o}\left(x\right)\right\}\right]\right\}
+\mathbb{E}_{P}\left[X\right]'b.
\end{split}\end{equation}
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$.
\begin{assumption}\label{assu:Px}
(i) $P$ and $Q$ have finite second moments, and (ii) $P$ is absolutely continuous with respect to the Lebesgue measure.
\end{assumption}
\begin{assumption}\label{assu:bilinearID}
The matrix $A$ in the production technology \eqref{eq: bi-linear production function} is invertible.
\end{assumption}
\noindent Assumptions \ref{assu:Px}--\ref{assu:bilinearID} serve as primitive conditions to satisfy Assumptions \ref{assu:den1}--\ref{assu:Phi1}, given our production technology \eqref{eq: bi-linear production function}. The following statement derives the equilibrium assignment and wage in terms of the solution to the dual Monge-Kantorovich problem \eqref{eq: dual-matching2}.
\begin{prop}
\label{Prop:Matchingiden}Let Assumption \ref{assu:Px} holds. Then, there exists the unique (up to constant) convex solution, $w^{o*}\left(x\right)$, to the second dual problem in \eqref{eq: dual-matching2}, and the equilibrium wage ($w^{*}$) and assignment ($y^{*}=T^{*}(x)$ where $y^{*}=(y_{1}^{*},\ldots,y_{d}^{*})'$) are given by
\[
w^{*}\left(x\right)=w^{o*}\left(x\right)+x'b+c,\quad
A\begin{pmatrix}y_{1}^{*}\left(x\right) & \cdots & y_{d}^{*}\left(x\right)\end{pmatrix}'
=\nabla w^{o*}\left(x\right),
\]
where $c$ is the constant of integration. In addition, if Assumption \ref{assu:bilinearID} holds,
\[
\begin{pmatrix}y_{1}^{*}\left(x\right) & \cdots & y_{d}^{*}\left(x\right)\end{pmatrix}'
=A^{-1}\nabla w^{o*}\left(x\right).
\]
\end{prop}
\noindent We can interpret this problem as assigning from $\mathcal{X}$ to
$A\mathcal{Y}:=\left\{AY:Y\in\mathcal{Y}\right\}$. Assumption \ref{assu:Px} guarantees the existence of the convex solution to the dual problem \eqref{eq: dual-matching2}.\footnote{Proposition \ref{prop:MatchingEq} 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{assu:Px}(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$. \citet{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, \cite{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{fig: lindenlaub} 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.
\begin{figure}[tbh]
\begin{center}
\begin{tikzpicture}
\node at (0,0) (X) {$x\sim P$};
\node at (0,-3) (Y) {$y\sim Q$};
\node at (5,0) (X1) {$\begin{array}{l}\tilde{x}_{C}=\Phi^{-1}\left(F_{x_{C}}\left(x_{C}\right)\right)\sim{\rm N}\left(0,1\right)\\ \tilde{x}_{M}=\Phi^{-1}\left(F_{x_{M}}\left(x_{M}\right)\right)\sim{\rm N}\left(0,1\right)\end{array}$};
\node at (5,-3) (Y1) {$\begin{array}{l}\tilde{y}_{C}=\Phi^{-1}\left(F_{y_{C}}\left(y_{C}\right)\right)\sim{\rm N}\left(0,1\right)\\ \tilde{y}_{M}=\Phi^{-1}\left(F_{y_{M}}\left(y_{M}\right)\right)\sim{\rm N}\left(0,1\right)\end{array}$};
\node at (11.5,0) (X2) {$\begin{pmatrix}\tilde{x}_{C}\\ \tilde{x}_{M}\end{pmatrix}\sim{\rm N}\left(0,\begin{pmatrix}1 & \rho_{\tilde{x}}\\ \rho_{\tilde{x}} & 1\end{pmatrix}\right)$};
\node at (11.5,-3) (Y2) {$\begin{pmatrix}\tilde{y}_{C}\\ \tilde{y}_{M}\end{pmatrix}\sim{\rm N}\left(0,\begin{pmatrix}1 & \rho_{\tilde{y}}\\ \rho_{\tilde{y}} & 1\end{pmatrix}\right)$};
\draw [->] (X) -- (Y);
\draw [->] (X) -- (X1);
\draw [->] (X1) -- (X2);
\draw [->] (Y) -- (Y1);
\draw [->] (Y1) -- (Y2);
\draw [->] (X2) -- (Y2);
\node [above] at (8.5,0) {?};
\node [above] at (8.5,-3) {?};
\node [right] at (0,-1.5) {$T^{*}\left(x\right)$};
\node [left] at (11.5,-1.5) {$\tilde{T}^{*}\left(\tilde{x}\right)$};
\end{tikzpicture}
\end{center}
\caption{\label{fig: lindenlaub}\citet{lindenlaub2017}'s transformation. $\Phi$ denotes the standard normal c.d.f. $F_{x_{C}}$ and $F_{x_{M}}$ denote the c.d.f. for $x_{C}$ and $x_{M}$, respectively.}
\end{figure}
The actual impact of technological shifts on wage distribution may differ from the prediction based on the Gaussian model. \citet{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{appen: skewdisp}), the obtained wage distribution's skewness does not reach its minimum when $\alpha_{CC}=\alpha_{MM}$.
\section{\label{sec: MatchingModelID}Empirical model and identification}
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.
\cite{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, \citet{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, \citet{sanders2014}, \citet{lindenlaub2017}, and \citet{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:
\begin{equation}\label{eq:OTmodel}\begin{split}
w_{i}&=w^{*}\left(x_{i}\right)+\varepsilon_{wi}
=w^{o*}\left(x_{i}\right)+x_{i}'b+c+\varepsilon_{wi},\\
y_{i}&=y_{i}^{*}+\varepsilon_{yi}
=A^{-1}\nabla w^{o*}\left(x_{i}\right)+\varepsilon_{yi}.
\end{split}\end{equation}
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 \cite{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$:
\begin{equation}\label{eq:OTexoX}
\mathbb{E}\left[\varepsilon_{wi}|x_{i}\right]=0,\quad
\mathbb{E}\left[\varepsilon_{yi}|x_{i}\right]=0.
\end{equation}
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
\begin{equation*}
\rho_{w}\left(w_{i},x_{i};\theta,w\right)=w_{i}-\left(w\left(x_{i}\right)+x_{i}'b\right),\quad
\rho_{y}\left(y_{i},x_{i};\theta,w\right)=y_{i}-A^{-1}\nabla w\left(x_{i}\right).
\end{equation*}
For each observation $i$, the model \eqref{eq:OTmodel} satisfies the moment conditions \eqref{eq:OTexoX}. This implies that the following conditional moments hold:
\begin{equation}\label{eq:CondMom}
\mathbb{E}\left[\rho\left(z_{i};\theta,w\right)|x_{i}\right]=0,
\end{equation}
at a true parameter $\left(\theta_{0},w_{0}\right)$. Then $\left(\theta_{0},w_{0}\right)$ are identified via the model \eqref{eq:CondMom} by Proposition \ref{Prop:Matchingiden} and the exogeneity of $x_i$ as well as the following assumption on $\mathcal{Y}$.
\begin{assumption}\label{assu:bilinearID2}
There exist $y_{1},\ldots,y_{d},y_{d+1}\in\mathcal{Y}$ such that $\left\{y_{1}-y_{2},\ldots,y_{d}-y_{d+1}\right\}$ is linearly independent.
\end{assumption}
Assumptions \ref{assu:Px} and \ref{assu:bilinearID}, combined with Proposition \ref{Prop:Matchingiden}, 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{assu:bilinearID2} is sufficient to identify the nonzero vector $b_{0}$, as stated in the following theorem.
\begin{theorem}
\label{thm:MatchingID} Let Assumptions \ref{assu:Px}-\ref{assu:bilinearID2} hold and the moment conditions \eqref{eq:CondMom} be satisfied. Then, $\theta_{0}$ and $w_{0}=w_{0}^{o*}+c_{0}$ are identified.
\end{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$.
\section{\label{sec:MatchingEst}Sieve-based semiparametric estimation}
The model parameters are identified by the semiparametric conditional moment restrictions \eqref{eq:CondMom}. 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,
{\small\begin{equation}\label{eq: funcspace}
\mathcal{W}_{n}
=\left\{w_{n}:\mathcal{X}\to\mathbb{R},
w_{n}\left(x;\gamma\right)=\sum_{j_{C}=0}^{k_{Cn}}\sum_{j_{M}=0}^{k_{Mn}}
\gamma_{j_{C}j_{M}}p_{j_{C}}\left(x_{C}\right)p_{j_{M}}\left(x_{M}\right),
\gamma_{j_{C}j_{M}}\in\mathbb{R}\right\},
\end{equation}}
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}$:
\begin{equation*}\begin{split}
\partial w_{n}\left(x;\gamma\right)/\partial x_{C}
&=\sum_{j_{C}=0}^{k_{Cn}}\sum_{j_{M}=0}^{k_{Mn}}
\gamma_{j_{C}j_{M}}\left(\partial p_{j_{C}}\left(x_{C}\right)/\partial x_{C}\right)
p_{j_{M}}\left(x_{M}\right),\\
\partial w_{n}\left(x;\gamma\right)/\partial x_{M}
&=\sum_{j_{C}=0}^{k_{Cn}}\sum_{j_{M}=0}^{k_{Mn}}
\gamma_{j_{C}j_{M}}p_{j_{C}}\left(x_{C}\right)
\left(\partial p_{j_{M}}\left(x_{M}\right)/\partial x_{M}\right).
\end{split}\end{equation*}
We first consider the model \eqref{eq:OTmodel} 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 \cite{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 \eqref{eq:OTmodel} 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
\begin{equation}\label{eq: SieveNLFI}
L\left(\theta,w\left(\cdot\right)\right)
=-\frac{n}{2}\log\det\left(\frac{1}{n}\sum_{i=1}^{n}\rho_{i}\rho_{i}'\right).
\end{equation}
The value of $\left(\theta,w\right)$ maximizing \eqref{eq: SieveNLFI} 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{assu:bilinearID}, we can apply the sieve generalized least squares (GLS) procedure \citep{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{tab: SGLS algorithm} outlines the three-step procedure to compute the SGLS estimator.
\begin{table}[ht!]
\caption{Three-step procedure for Sieve GLS estimation \citep{chen2007}}\label{tab: SGLS algorithm}
\noindent \centering{}
\begin{tabular}{l}
\hline
\textbf{\textsc{Algorithm:}}\textbf{ }Computing the Sieve GLS Estimator
of $\theta$ and $w$\tabularnewline
\hline
1. Obtain an initial consistent sieve LS estimator $\left(\tilde{\theta}_{n},\tilde{w}_{n}\right)$
by\tabularnewline
$\quad \ \min_{\left(\theta,w\right)}\sum_{i=1}^{n}\rho\left(z_{i};\theta,w\right)^{\prime}\rho\left(z_{i};\theta,w\right),$\tabularnewline
2. Obtain a consistent estimator $\hat{\Sigma}_{0}\left(x\right)$ of
$\Sigma_{0}\left(x\right)={\rm Var}\left[\left.\rho\left(z_{i};\theta,w\right)\right|x_{i}=x\right]$\tabularnewline
$\quad \ $using $\left(\tilde{\theta}_{n},\tilde{w}_{n}\right)$ and
sieve LS estimation.\tabularnewline
3. Obtain the optimally weighted sieve GLS estimator $\left(\hat{\theta}_{n},\hat{w}_{n}\right)$
by\tabularnewline
$\quad \ \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)$.\tabularnewline
\hline
\end{tabular}
\end{table}
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{tab:comparison-estimators}.
\begin{table}[ht!]
\caption{Comparison of key assumptions in estimation procedures}
\label{tab:comparison-estimators}
\centering
\begin{tabular}{lcccc}
\hline
\multirow{2}{*}{Assumptions} & \cite{lindenlaub2017} & \multicolumn{3}{c}{Sieve Estimators} \\
\cmidrule(lr){2-2} \cmidrule(lr){3-5}
& ML & SML & SLS & SGLS \\
\hline
Joint normality of $x$ and $y$ & {\checkmark} & {} & {} & {} \\
Normality of measurement errors & {\checkmark} & {\checkmark} & {} & {} \\
Uncorrelated measurement errors & {\checkmark} & {\checkmark} & {} & {} \\
Homoskedasticity of measurement errors & {\checkmark} & {\checkmark} & {\checkmark} & {} \\
\hline
\end{tabular}
\end{table}
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., \citet{lorentz1986}).
This provides an approach to imposing shape restrictions on the sieve estimator with a linear constraint which can be solved easily.\footnote{\citet{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{appen: bernstein convexity}.
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.
\section{Asymptotic theory for the SGLS estimator}\label{sec: asymptotics}
We establish consistency, convergence rate, asymptotic normality, and semiparametric efficiency of our SGLS estimator using results in \cite{chen1998}, \cite{ai2003}, and \cite{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]
=\begin{pmatrix}
w_{1}\left(x\right)-w_{2}\left(x\right)+x'\left(b_{1}-b_{2}\right) \\
\nabla_{C}w_{0}\left(x\right)\left(\kappa_{C1}-\kappa_{C2}\right)
+\kappa_{C0}\left(\nabla_{C}w_{1}\left(x\right)-\nabla_{C}w_{2}\left(x\right)\right)\\
\nabla_{M}w_{0}\left(x\right)\left(\kappa_{M1}-\kappa_{M2}\right)
+\kappa_{M0}\left(\nabla_{M}w_{1}\left(x\right)-\nabla_{M}w_{2}\left(x\right)\right)
\end{pmatrix}.
\]
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. \cite{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.
\begin{assumption}\label{assu:IID}
(i) $\left\{w_{i},y_{i}',x_{i}'\right\}_{i=1}^{n}$ are i.i.d.; (ii) $\mathcal{X}$ is compact and a Cartesian product of compact intervals $\mathcal{X}_{C}$ and $\mathcal{X}_{M}$.
\end{assumption}
\begin{assumption}\label{assu:Sigma}
$\Sigma\left(x\right)$ and $\Sigma_{0}\left(x\right)\equiv\mathrm{Var}\left(\rho\left(z_{i},\lambda_{0}\right)|x_{i}=x\right)$ are positive definite and bounded uniform over $x\in\mathcal{X}$.
\end{assumption}
\begin{assumption}\label{assu:compact}
$\Lambda\equiv\Theta\times\mathcal{W}$ is compact under $\lVert\cdot\rVert_{s}$.
\end{assumption}
\begin{assumption}\label{assu:Holder}
(i) $w\in\Lambda^{m}\left(\mathcal{X}\right)$ with $m>2$; (ii) $\forall w\in\Lambda^{m}\left(\mathcal{X}\right), \exists w_{n}\left(x;\gamma\right)\in\mathcal{W}_{n}$ such that $\lVert w_{n}-w\rVert_{\infty}=O(\left(k_{Cn}k_{Mn}\right)^{-m/2})$ with $k_{Cn},k_{Mn}=O\left(n^{1/2\left(m+1\right)}\right)$.
\end{assumption}
Assumptions \ref{assu:IID}--\ref{assu:compact} 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{assu:Px}--\ref{assu:bilinearID} guarantee it by Theorem \ref{thm:MatchingID}. Assumption \ref{assu:Holder} 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 \citet{caffarelli1992CPAM,caffarelli1992JAMS,caffarelli1996} which covers the case of compactly supported $\mathcal{X}$ and $\mathcal{Y}^{*}.$ \citet{cordero-erausquin2019} provides an extended result for distributions with unbounded supports.} The following proposition establishes the convergence rate of $\hat{\lambda}_n.$
\begin{prop}\label{prop:convrate}
If Assumptions \ref{assu:Px}-\ref{assu:Holder} hold, then $\lVert\hat{\lambda}_{n}-\lambda_{0}\rVert=o_{p}\left(n^{-1/4}\right)$.
\end{prop}
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
\begin{align*}
D_{v_{1}}\left(x\right)
&=\begin{pmatrix}v_{1}\left(x\right)\\
\kappa_{C0}\nabla_{C}v_{1}\left(x\right)-\nabla_{C}w_{0}\left(x\right)\\
\kappa_{M0}\nabla_{M}v_{1}\left(x\right)\end{pmatrix},\
&D_{v_{2}}\left(x\right)
&=\begin{pmatrix}v_{2}\left(x\right)\\
\kappa_{C0}\nabla_{C}v_{2}\left(x\right)\\
\kappa_{M0}\nabla_{M}v_{2}\left(x\right)-\nabla_{M}w_{0}\left(x\right)\end{pmatrix},\\
D_{v_{3}}\left(x\right)
&=\begin{pmatrix}v_{3}\left(x\right)-x_{C}\\
\kappa_{C0}\nabla_{C}v_{3}\left(x\right)\\
\kappa_{M0}\nabla_{M}v_{3}\left(x\right)\end{pmatrix},\
&D_{v_{4}}\left(x\right)
&=\begin{pmatrix}v_{4}\left(x\right)-x_{M}\\
\kappa_{C0}\nabla_{C}v_{4}\left(x\right)\\
\kappa_{M0}\nabla_{M}v_{4}\left(x\right)\end{pmatrix}.
\end{align*}
Let $v^{*}=\left(v_{1}^{*},v_{2}^{*},v_{3}^{*},v_{4}^{*}\right)$, where $v_{j}^{*}$ solves
\begin{equation}\label{eq: v*}
\inf_{v_{j}}\mathbb{E}\left[D_{v_{j}}\left(x_{i}\right)'
\Sigma\left(x_{i}\right)^{-1}D_{v_{j}}\left(x_{i}\right)\right].
\end{equation}
\begin{assumption}\label{assu:Dv*}
(i) $\mathbb{E}\left[D_{v^{*}}\left(x_{i}\right)'D_{v^{*}}\left(x_{i}\right)\right]$ is positive definite; (ii) Each element of $v^{*}$ belongs to the H\"{o}lder space $\Lambda^{m}\left(\mathcal{X}\right)$ with $m>2$.
\end{assumption}
\begin{assumption}\label{assu:interior}
$\theta_{0}\in\mathrm{int}\left(\Theta\right)$.
\end{assumption}
Under Assumptions \ref{assu:Px}--\ref{assu:Dv*}, it is clear to see from Lemma B.1 in \citet{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.$
\begin{theorem}\label{thm:AsympSieveGLSE}
Let Assumptions \ref{assu:Px}--\ref{assu:interior} hold. Then, $\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\rightarrow_{d}N\left(0,V_{1}^{-1}V_{2}V_{1}^{-1}\right)$, where
\begin{equation}\label{eq: V1V2}\begin{split}
V_{1}
&=\mathbb{E}\left[D_{v^{*}}\left(x_{i}\right)'\Sigma\left(x_{i}\right)^{-1}D_{v^{*}}\left(x_{i}\right)\right],\\
V_{2}
&=\mathbb{E}\left[D_{v^{*}}\left(x_{i}\right)'
\Sigma\left(x_{i}\right)^{-1}\Sigma_{0}\left(x_{i}\right)\Sigma\left(x_{i}\right)^{-1}
D_{v^{*}}\left(x_{i}\right)\right].
\end{split}\end{equation}
\end{theorem}
The asymptotic variance $V_{1}^{-1}V_{2}V_{1}^{-1}$ can be consistently estimated (see, Remark 4.2 in \citet{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{thm:AsympSieveGLSE} 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.
\begin{assumption}\label{assu:Sigma2}
$\hat{\Sigma}\left(x\right)=\Sigma_{0}\left(x\right)+o_{p}\left(n^{-1/4}\right)$ uniformly over $x\in\mathcal{X}$.
\end{assumption}
Let $v_{0}=\left(v_{01},v_{02},v_{03},v_{04}\right)$, where $v_{0j}$ solves \eqref{eq: v*} with $\Sigma\left(x\right)$ replaced by $\Sigma_{0}\left(x\right)$. Now the following theorem establishes the semiparametric efficiency of $\hat{\theta}_n$.
\begin{theorem}\label{thm:EfficSieveGLSE}
Suppose that all conditions of Theorem \ref{thm:AsympSieveGLSE} with $\Sigma\left(x\right)=\Sigma_{0}\left(x\right)$ and $v^{*}=v_{0}$ hold, and Assumption \ref{assu:Sigma2} is satisfied. Then, $\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\rightarrow_{d}N\left(0,V_{0}^{-1}\right)$, with $V_{0}=\mathbb{E}\left[D_{v_{0}}\left(x_{i}\right)'\Sigma_{0}\left(x_{i}\right)^{-1}D_{v_{0}}\left(x_{i}\right)\right]$.
\end{theorem}
\section{\label{sec: Simul}Monte Carlo simulations}
This section evaluates the finite sample performances of our sieve estimators using known data-generating processes (DGPs). We first generate Monte Carlo samples from \cite{lindenlaub2017}'s quadratic-Gaussian model. Workers' skill bundle, $x$, and occupations' skill requirements, $y$, follow joint Gaussian distributions:
$$ \begin{pmatrix} x_C\\x_M \end{pmatrix}
\sim N\left(\begin{pmatrix} 0\\0 \end{pmatrix},
\begin{pmatrix} 1&\rho_x \\ \rho_x&1 \end{pmatrix}\right),\quad
\begin{pmatrix} y_C\\y_M \end{pmatrix}
\sim N\left(\begin{pmatrix} 0\\0 \end{pmatrix},
\begin{pmatrix} 1&\rho_y \\ \rho_y&1 \end{pmatrix}\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^*=\begin{pmatrix} y_C^*\\y_M^* \end{pmatrix}
=\underbrace{\begin{pmatrix} J_{11}&J_{12} \\ J_{21}&J_{22} \end{pmatrix}}_{:=J}
\begin{pmatrix} x_C\\x_M \end{pmatrix},$$
where
$$ J=\frac{1}{\sqrt{1+2\delta(\rho_x\rho_y+\sqrt{1-\rho_y^2}\sqrt{1-\rho_x^2})+\delta^2}}
\begin{pmatrix}
1+\delta\frac{\sqrt{1-\rho_y^2}}{\sqrt{1-\rho_x^2}} &
\delta\left(\rho_y-\rho_x\frac{\sqrt{1-\rho_y^2}}{\sqrt{1-\rho_x^2}}\right) \\
\rho_y-\rho_x\frac{\sqrt{1-\rho_y^2}}{\sqrt{1-\rho_x^2}} &
\delta+\frac{\sqrt{1-\rho_y^2}}{\sqrt{1-\rho_x^2}}
\end{pmatrix}.$$
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 \cite{lindenlaub2017}.
We estimate the production technology parameters using \cite{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 \cite{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.
\begin{figure}[ht!]
\begin{centering}
\includegraphics[width=\textwidth]{figures/simul_performance_gaussian.pdf}
\par\end{centering}
\caption{\label{fig:simul performance}Box plots of parameter estimates (Gaussian DGP)}
\end{figure}
The box plots in Figure \ref{fig:simul performance} 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{tab:simul performance}. 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.
\begin{table}[ht!]
\centering\footnotesize
\caption{\label{tab:simul performance}Finite sample performances of the estimators (Gaussian DGP)}
\begin{tabular}{ccccccc}
\toprule
&& ML & ML$^*$ & SML & SLS & SGLS \\
\midrule
$\alpha_{CC}$ & Bias& -0.0536 & -0.0041 & -0.0007 & -0.0018 & -0.0027\\
& RMSE& 0.0775 & 0.0513 & 0.0520 & 0.0523 & 0.0523\\
\hline
$\alpha_{MM}$ & Bias& 0.0992 & 0.0044 & 0.0065 & 0.0031 & 0.0022\\
& RMSE& 0.1077 & 0.0473 & 0.0491 & 0.0502 & 0.0489\\
\hline
$\beta_{C}$ & Bias& -0.0006 & -0.0007 & -0.0009 & -0.0009 & -0.0009\\
& RMSE& 0.0416 & 0.0414 & 0.0427 & 0.0427 & 0.0427\\
\hline
$\beta_{M}$ & Bias& -0.0000 & -0.0001 & -0.0000 & -0.0000 & -0.0000\\
& RMSE& 0.0398 & 0.0395 & 0.0397 & 0.0397 & 0.0397\\
\bottomrule
\end{tabular}
\end{table}
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:=\begin{bmatrix} x_{C1}&x_{M1} \\ \vdots&\vdots \\ x_{Cn}&x_{Mn} \end{bmatrix},\quad
y:=\begin{bmatrix} y_{C1}&y_{M1} \\ \vdots&\vdots \\ y_{Cn}&y_{Mn} \end{bmatrix}.$$
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:
\begin{equation}
\max_{x} f'x,
\quad{\rm s.t.}\
\begin{bmatrix}
\underbrace{A_1}_{n \times n^2} \\ \underbrace{A_2}_{n \times n^2}
\end{bmatrix}x \le \underbrace{b}_{2n \times 1},\label{eq: linear programming}
\end{equation}
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 \eqref{eq: linear programming}.\footnote{Even with the moderate sample size $n=3000$, the constraint matrix is enormous ($6000 \times 9,000,000).$ Solving the linear program \eqref{eq: linear programming} over many Monte Carlo samples is computationally demanding. We employ \textsc{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:
$$ \begin{pmatrix} \varepsilon_w \\ \varepsilon_C \\ \varepsilon_M \end{pmatrix}
\sim N\left(\begin{pmatrix} 0\\0\\0 \end{pmatrix},
\begin{pmatrix} 2&1&1 \\ 1&1&0.5 \\ 1&0.5&1 \end{pmatrix}\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.
\begin{table}[ht!]
\centering\footnotesize
\caption{\label{tab:simul performance gumbel}Finite sample performances of the estimators (Gumbel DGP)}
\begin{tabular}{cccccccccc}
\toprule
& &\multicolumn{4}{c}{Gamma errors} & \multicolumn{4}{c}{Joint Gaussian errors} \\
\cmidrule(lr){3-6} \cmidrule(lr){7-10}
&& ML$^*$ & SML & SLS & SGLS & ML$^*$ & SML & SLS & SGLS \\
\midrule
$\alpha_{CC}$ &Bias& -0.0608 & -0.0021 & -0.0015 & -0.0020 & -0.2032 & -0.0023 & 0.0008 & 0.0112 \\
&RMSE& 0.3780 & 0.0538 & 0.0535 & 0.0538 & 0.2971 & 0.0914 & 0.0925 & 0.0809 \\
\hline
$\alpha_{MM}$ &Bias& -0.0550 & 0.0019 & 0.0024 & 0.0020 & -0.3339 & 0.0018 & 0.0034 & -0.0123 \\
&RMSE& 0.4209 & 0.0512 & 0.0527 & 0.0513 & 1.6626 & 0.0886 & 0.0898 & 0.0827 \\
\hline
$\beta_{C}$ &Bias& 0.0055 & 0.0003 & 0.0003 & 0.0003 & -0.0015 & -0.0054 & -0.0054 & -0.0024 \\
&RMSE& 0.1138 & 0.0406 & 0.0405 & 0.0406 & 0.1013 & 0.0762 & 0.0757 & 0.0760 \\
\hline
$\beta_{M}$ &Bias& -0.0056 & 0.0011 & 0.0011 & 0.0011 & -0.0327 & -0.0058 & -0.0056 & -0.0032 \\
&RMSE& 0.1033 & 0.0398 & 0.0398 & 0.0399 & 0.1137 & 0.0735 & 0.0728 & 0.0677 \\
\bottomrule
\end{tabular}
\end{table}
The estimation results are provided in Table \ref{tab:simul performance gumbel}. 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 \cite{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(\begin{pmatrix} 1\\1 \end{pmatrix},
\begin{pmatrix} 1&\rho \\ \rho&1\end{pmatrix}\right),\quad
K_2 \sim N\left(\begin{pmatrix} -1\\-1 \end{pmatrix},
\begin{pmatrix} 1&-\rho \\ -\rho&1 \end{pmatrix}\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(\begin{pmatrix} 1\\1\\1 \end{pmatrix},
\begin{pmatrix} 1&0.7&0.7 \\ 0.7&1&0.3 \\ 0.7&0.3&1 \end{pmatrix}\right),\quad
M_2 \sim N\left(\begin{pmatrix} -3\\-3\\-3 \end{pmatrix},
\begin{pmatrix} 1&0.7&0.7 \\ 0.7&1&0.3 \\ 0.7&0.3&1 \end{pmatrix}\right).$$
In this case, both $(x,y)$ and measurement errors have bi-modal distributions that are far from a normal distribution.
\begin{table}[ht!]
\centering\footnotesize
\caption{\label{tab:simul performance mixture}Finite sample performances of the estimators (Gaussian mixture DGP)}
\begin{tabular}{cccccc}
\toprule
&& ML$^*$ & SML & SLS & SGLS \\
\midrule
$\alpha_{CC}$ & Bias& 0.2896 & -0.0014 & -0.0016 & -0.0017 \\
& RMSE& 0.3653 & 0.0429 & 0.0425 & 0.0385 \\
\hline
$\alpha_{MM}$ & Bias& 0.2446 & -0.0017 & 0.0006 & -0.0027 \\
& RMSE& 0.3584 & 0.0454 & 0.0431 & 0.0337 \\
\hline
$\beta_{C}$ & Bias& 0.7123 & -0.0007 & -0.0013 & 0.0008 \\
& RMSE& 0.7204 & 0.0428 & 0.0426 & 0.0364 \\
\hline
$\beta_{M}$ & Bias& -0.1288 & -0.0001 & 0.0002 & -0.0011 \\
& RMSE& 0.1592 & 0.0421 & 0.0419 & 0.0305 \\
\bottomrule
\end{tabular}
\end{table}
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{tab:simul performance mixture}. 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.
\section{\label{sec: Matchingemp}Empirical application to U.S. worker-job matching}
We revisit the dataset constructed by \cite{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. \citet{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 \cite{lindenlaub2017}. The original supports of worker skills and firms' skill demands are provided in Table \ref{tab: statistics-matching}.} 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. \cite{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{tab: statistics-matching} 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.
\begin{table}[t!]
\caption{\label{tab: statistics-matching}Summary statistics of worker skills ($x$) and firms' skill demands ($y$)}
\centering{}
\begin{tabular}{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\hline
& \multicolumn{8}{c}{1990/91 ($n=2984$)} & \multicolumn{8}{c}{2009/10 ($n=4495$)}\tabularnewline
\cmidrule(lr){2-9} \cmidrule(lr){10-17}
& \multicolumn{2}{c}{$x_{C}$} & \multicolumn{2}{c}{$x_{M}$} & \multicolumn{2}{c}{$y_{C}$} & \multicolumn{2}{c}{$y_{M}$} & \multicolumn{2}{c}{$x_{C}$} & \multicolumn{2}{c}{$x_{M}$} & \multicolumn{2}{c}{$y_{C}$} & \multicolumn{2}{c}{$y_{M}$}\tabularnewline
\hline
{Mean} & 0&3596 & -0&2912 & 0&0135 & -0&1189 & 0&5667 & -0&6601 & 0&0468 & -0&2509\tabularnewline
{SD} & 0&7423 & 0&9923 & 0&8490 & 1&0240 & 0&7556 & 0&8358 & 0&9280 & 0&9656\tabularnewline
{Min} & -2&0595 & -1&7004 & -2&0622 & -1&6949 & -2&3019 & -1&8116 & -2&5200 & -1&6597\tabularnewline
{Max} & 2&1649 & 2&1855 & 2&0925 & 2&1895 & 1&9160 & 2&1838 & 3&0504 & 2&1351\tabularnewline
\hline
\end{tabular}
\end{table}
\cite{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{fig:dist-transformed}, 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.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/2dxy.png}
\par\end{centering}
\caption{\label{fig:dist-transformed}Joint densities of $\tilde{x}$ and $\tilde{y}$ (transformed data)}
\end{figure}
We employ Mardia's test \citep{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{tab: mardia-transformed} 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.
\begin{table}[tbh]
\caption{\label{tab: mardia-transformed}Mardia's multivariate normality test statistics (p-values in parentheses)}
\centering{}
\begin{tabular}{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\hline
& \multicolumn{4}{c}{1990/91 ($n=2984$)} & \multicolumn{4}{c}{2009/10 ($n=4495$)}\tabularnewline
\cmidrule(lr){2-5} \cmidrule(lr){6-9}
& \multicolumn{2}{c}{$\tilde{x}$} & \multicolumn{2}{c}{$\tilde{y}$} & \multicolumn{2}{c}{$\tilde{x}$} & \multicolumn{2}{c}{$\tilde{y}$}\tabularnewline
\hline
Skewness & 4&58 (0.333) & 100&09 (0.000) & 16&34 (0.003) & 145&14 (0.000)\tabularnewline
Kurtosis & 4&44 (0.000) & 0&29 (0.774) & 14&42 (0.000) & 1&98 (0.048)\tabularnewline
\hline
\end{tabular}
\end{table}
We estimate the production technology in each period separately. First, we modify \cite{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.
\begin{table}[t!]
\caption{\label{tab: est-matching-trans}Estimates of production technology parameters on transformed data}
\centering{}
{\footnotesize
\begin{tabular}
{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\hline
& \multicolumn{10}{c}{1990/91} & \multicolumn{10}{c}{2009/10}\tabularnewline
\cmidrule(lr){2-11} \cmidrule(lr){12-21}
& \multicolumn{2}{c}{ML} &\multicolumn{2}{c}{ML*} & \multicolumn{2}{c}{SML} & \multicolumn{2}{c}{SLS} & \multicolumn{2}{c}{SGLS} &\multicolumn{2}{c}{ML} &\multicolumn{2}{c}{ML*} & \multicolumn{2}{c}{SML} & \multicolumn{2}{c}{SLS} & \multicolumn{2}{c}{SGLS}\tabularnewline
\hline
$\alpha_{CC}$ & 0&203 & 0&765 & 0&454 & 0&000 & 0&000 & 0&739 & 1&119 & 2&048 & 2&293 & 2&290\tabularnewline
& (0&342)&(0&574)&(0&009)& (0&000)&(0&000) &(0&198)& (0&408)& (0&036)&(0&256)&(0&265)\tabularnewline
$\alpha_{MM}$ & 0&479 & 1&270 & 1&422 & 1&084 & 0&856 & 0&055 & 0&486 & 0&237 & 1&033 & 0&291\tabularnewline
& (0&175)&(0&148)&(0&031)& (0&113)&(0&098) &(0&154)& (0&633)& (0&006)&(0&215)&(0&059)\tabularnewline
$\beta_{C}$ & 1&686 & 1&711 & 1&692 & 1&719 & 1&585 & 2&203 & 2&208 & 2&115 & 2&063 & 2&198\tabularnewline
& (0&143)&(0&589)&(0&068)& (0&434)&(0&416) &(0&152)& (0&540)& (0&068)&(0&589)&(0&534)\tabularnewline
$\beta_{M}$ & -0&421 & -0&392 &-0&374 & -0&382 &-0&388 & 0&210 & 0&243 & 0&198 & 0&180 & 0&327\tabularnewline
& (0&141)&(0&406)&(0&068)& (0&329)&(0&309) &(0&152)& (0&731)& (0&076)&(0&545)&(0&532)\tabularnewline
\hline
\multicolumn{21}{l}{{\footnotesize Standard errors in the parentheses. ML indicates the original estimates in \cite{lindenlaub2017}.}}\tabularnewline
\end{tabular}}
\end{table}
The estimation results are provided in Table \ref{tab: est-matching-trans}. 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.
\begin{figure}[b!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/wage_polarization.pdf}
\par\end{centering}
\caption{\label{fig:wagepol}Actual and model predicted wage polarization (transformed data)}
\end{figure}
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 \textit{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{fig:wagepol} 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{fig:wagepol} (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.
\begin{figure}[b!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/wage_curvature.pdf}
\par\end{centering}
\caption{\label{fig:wage3d}Estimated wage functions (transformed data)}
\end{figure}
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{fig:wage3d}--\ref{fig:wage curvature}. 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.
\begin{figure}[t!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/wage_function.pdf}
\par\end{centering}
\caption{\label{fig:wage curvature}Predicted wage with respect to cognitive skill (transformed data)}
\end{figure}
To understand the driving forces behind wage polarization, we isolate the effects of technological and distributional changes in Figure \ref{fig:wagepol}. 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.
\begin{figure}[b!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/skill_dist_UT1.pdf}
\medskip
\includegraphics[width=1\textwidth]{figures/skill_dist_UT2.pdf}
\par\end{centering}
\caption{\label{fig:marginal dist}Marginal distributions of skill supply and demand}
\end{figure}
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{fig: lindenlaub}, 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{fig:marginal dist} 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{tab: est-matching-original}, 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.
\begin{table}[tbh]
\caption{\label{tab: est-matching-original}Estimates of production technology parameters on original data}
\centering{}
\footnotesize
\begin{tabular}
{cr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\hline
& \multicolumn{6}{c}{1990/91} & \multicolumn{6}{c}{2009/10}\tabularnewline
\cmidrule(lr){2-7} \cmidrule(lr){8-13}
& \multicolumn{2}{c}{SML} & \multicolumn{2}{c}{SLS} & \multicolumn{2}{c}{SGLS} & \multicolumn{2}{c}{SML} & \multicolumn{2}{c}{SLS} & \multicolumn{2}{c}{SGLS}\tabularnewline
\hline
$\alpha_{CC}$ & 0&000 & 0&000 & 0&000 & 6&032 & 5&784 & 6&471\tabularnewline
&(0&000) &(0&000)& (0&000)& (0&110)&(0&184)&(0&196)\tabularnewline
$\alpha_{MM}$ & 1&947 & 2&292 & 2&234 & 0&000 & 1&020 & 0&003\tabularnewline
&(0&060) &(0&081)& (0&079)& (0&000)&(0&024)&(0&000)\tabularnewline
$\beta_{C}$ & 2&238 & 2&246 & 2&108 & 3&395 & 3&299 & 3&707\tabularnewline
&(0&083) &(0&225)& (0&220)& (0&101)&(0&396)&(0&292)\tabularnewline
$\beta_{M}$ &-0&341 & -0&441 & -0&317 & 0&254 & 0&384 & 0&440\tabularnewline
&(0&073) &(0&211)& (0&213)& (0&097)&(0&333)&(0&239)\tabularnewline
\hline
\multicolumn{13}{l}{{\footnotesize{}Standard errors in the parentheses}}\tabularnewline
\end{tabular}
\end{table}
\begin{figure}[h!]
\begin{centering}
\includegraphics[width=1\textwidth]{figures/wage_polarization_UT.pdf}
\par\end{centering}
\caption{\label{fig:wagepol-UT}Actual and model predicted wage polarization (original data)}
\end{figure}
The estimated models on the original data effectively capture the patterns of wage polarization, particularly in the upper tail as in Figure \ref{fig:wagepol-UT}. 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.
\section{Conclusion}\label{sec: conclusion}
Theoretical matching models often face empirical challenges due to discrepancies between model assumptions and real-world data. \cite{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 \citep{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 \citet{eeckhout2011}, to multidimensional settings.
\bibliographystyle{ecta-fullname}
\bibliography{main}
\newpage