Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
50,256 characters · 15 sections · 6 citation commands
Estimating Separable Matching Models
{ Keywords: matching, marriage, assignment, estimations comparison.} \newline { JEL codes: C78, C13, C15.}
The estimation of models of two-sided matching has made considerable progress in the past decade. While some of this work has used matching under non-transferable utility, many applications have focused on markets where utility is transferable. The pioneering contribution of \citecomplete{choo-siow:06} introduced a simple and highly tractable specification. They used their model to estimate the effect of the 1973 liberalization of abortion in the US on marriage outcomes. In doing so, they used a nonparametric estimator of the matching patterns. Their specification is a natural extension of the multinomial logit model, and it has become quite popular.
The choo-siow:06 specification rests on three main assumptions that will be defined later in the paper: separability; large market; and standard type I extreme value random utility. In \citecomplete{cupid:restud}, we showed that the third, distributional assumption is not necessary: for any (separable) distribution of the errors, the joint surplus is nonparametrically identified. The nonparametric estimator of choo-siow:06 was feasible in their case as they only conditioned on the ages of the partners in a couple. It breaks down, however, when more covariates are considered as matching cells become too small; and by construction, it does not allow for parameterized error distributions. Structural models of household behavior also naturally introduce parameters.
In all of these cases, the analyst must resort to parametric models. This note shows two very simple methods to estimate parametric versions of separable matching models with perfectly transferable utility, with special emphasis on the choo-siow:06 model and more generally on “semilinear” models, where the joint surplus is linear in the parameters (again, to be formally defined later in the paper).
Our first method applies a minimum-distance estimator to the identification equation derived in \citecomplete{cupid:restud}, which relates the joint surplus to the derivatives of a generalized entropy function evaluated at the observed matching patterns. For any fixed distribution of the error terms, the generalized entropy can be evaluated and differentiated, numerically if needed. The estimator selects parameter values and also provides a simple specification test. In semilinear models, the estimator can be obtained in closed form.
The second method we present applies more specifically to the semilinear Choo and Siow model. We show that the moment-matching estimator we described in \citecomplete{cupid:restud} can be reframed as a generalized linear model, more specifically as the pseudo-maximum likelihood estimator of a Poisson regression with two-sided fixed effects. This is available as linear_model in the scikit-learn library in Python, as fepois in the R package fixest and as \texttt{ ppmlhdfe} in Stata, among other common statistical packages.
We conclude with a brief discussion of the pros and cons of these two methods. Both are coded in a Python package called cupid_matching that is available on the standard repositories\footnote{ See \url{http://bsalanie.github.io} for more information, and \url{https://share.streamlit.io/bsalanie/cupid_matching_st/main/cupid_streamlit.py} for an interactive Streamlit app that demonstrates solving and estimating a \citecomplete{choo-siow:06} model.}.
This paper applies to a bipartite matching market with perfectly transferable utility. For simplicity, we refer to potential partners as “men” and “women”. We use the same notation as in \citecomplete{cupid:restud}. We assume that the analyst can only observe which of a finite set of types each individual belongs to. Men and women of a given type differ along some dimensions that they all observe, while the analyst does not. Each man $i\in \mathcal{I}$ belongs to one group of (observable) type $x_{i}\in \mathcal{X}$; and, similarly, each woman $j\in \mathcal{J }$ belongs to one (observable) type $y_{j}\in \mathcal{Y}$. We will say that “man $i$ is of type $x$” and “woman $j$ is of type $y$.” We denote $\mu_{ij}$ the indicator function for a matching between man $i$ and woman $j$, which is equal to $1$ if $i$ and $j$ are matched and to $0$ otherwise. Similarly, $ \mu_{i0}$ and $\mu_{0j}$ are the indicator of $i$ or $j$ to remain unmatched, respectively. Without loss of generality, we assimilate $\mathcal{ X}$ to $\{1,\ldots,X\}$ and $\mathcal{Y}$ to $\{1,\ldots,Y\}$. As $\mathcal X$ and $\mathcal Y$ will later serve as choice sets of partners types for women and men, respectively, and as the marital options needs to include remaining unmatched, we shall add the option to remain unmatched $0$ to these sets and denote $\mathcal X_{0} = \mathcal X \cup \{0\}$ and $\mathcal Y_{0} = \mathcal Y \cup \{0\}$ the respective sets of marital options of women and men.
We denote $n_{x}$ the mass of men of type $x\in \mathcal X$, and $m_{y}$ the mass of women of type $y \in \mathcal Y$. We denote $\bm{q}=(\bm{n}, \bm{m})$ the vector that collects the margins $\bm{n} $ and $\bm{m}$ of the problem
In addition to the margins $\bm{q}$, the analyst observes matchings at the type level. We denote $\mu_{xy}$ the mass of the couples where the man belongs to type $x$, and where the woman belongs to type $y$, which is formally defined as $\mu_{xy} = \sum_{ i\in \mathcal{I}, j\in\mathcal{J}} \mu_{ij}1\{ x_i = x\}1\{ y_j = y\}$. We also denote $\mu_{x0} = \sum_{ i\in \mathcal{I}} \mu_{i0}1\{ x_i = x\}$, and $\mu_{0y} = \sum_{ j\in\mathcal{J}} \mu_{0j}1\{ y_j = y\}$ the mass of single individuals who are respectively men of type $x$ and women of type $y$. We will be interested in the limiting market with a large number of men in any type $x$, and of women in any type $y$. Since the problem is homogeneous, we shall normalize the total mass $N$ of households to one; that is, we rescale $\bm{\mu}$ and $\bm{q}$ by a multiplicative factor $N$ such that $ \sum_{x,y}\mu_{xy}+\sum_x \mu_{x0}+\sum_y \mu_{0y} =1$. Again, we use the boldface notation $\bm{\mu}$ to denote the vector of matching numbers. We denote $\mathcal A = \mathcal X \times \mathcal Y \cup \mathcal X \times \{ 0 \} \cup \{0 \} \times \mathcal Y$ the set of possible marital arrangements (matched household of type $xy$, or single households of type $x0$ or $0y$), so that $\bm{\mu}$ is a vector of $\mathbb R^{\mathcal A}$.
A matching is the specification of who matches with whom. It is feasible if each individual is matched to at most one partner. It is stable if no individual who has a partner would prefer to be single, and if no two individuals would prefer forming a couple over their current situation.
We model the joint surplus $\tilde{\Phi}_{ij}$, which is the sum of the cardinal utilities that both a man $i$ and a woman $j$ jointly obtain by being matched together, and we assume a separable matching surplus:
Separability allows for a restricted form of “matching on unobservables”; it rules out interaction terms on characteristics that are unobserved on both sides of the market, e.g.\ some unobserved preference of man $i$ for some unobserved characteristics of woman $j$.
\citecomplete{csw:17} and \citecomplete{cupid:restud} showed that under separability, at any stable matching $\bm{\mu}$ there exist two matrices $\bm{U}$ and $\bm{V}$ such that for all $(x,y)$, $U_{xy}+V_{xy} = \Phi_{xy}$, and $U_{x0}=V_{0y}=0$ , and such that man $i$ of type $x$ is assigned option $y=0,1,\ldots,Y$ which maximizes $U_{xy}+\varepsilon_{iy}$ (where option $0$ means remaining unmatched, and option $y\neq 0$ means being matched with a woman of type $y$); similarly woman $j$ of type $y$ is assigned option $x=0,1,\ldots,X$ which maximizes $V_{xy}+\eta_{xj}$.
Consider the classic “Emax” function $G_x$ defined as follows. In this paragraph we let $\bm{U} =(U_1,\ldots,U_{Y})$ be a $Y$-dimensional vector. Then we define
As a maximum of linear functions, $G_x$ is a convex function. We denote $ \partial G_x(\bm{U})$ its subgradient; because of the assumptions made on $ \mathbb{P}_{x}$, it is a singleton almost everywhere.
Now take the Legendre-Fenchel transform of $G_x$: for any $(\nu_1,\ldots, \nu_{Y})$ such that $\sum_{y\in \mathcal{Y} }\nu_y\leq 1$, we define
It is another convex function; and since $G_x$ is convex, $G_x$ is the Legendre-Fenchel transform of $G^\ast_x$. As a consequence,
This convex duality is at the core of the identification and inference results in \citecomplete{cupid:restud}.
Defining $H_y$ and $H^\ast_y$ in the same way, we get the generalized entropy: for any feasible matching $\bm{\mu}$,
The function $\mathcal{E}$ only depends on the matching patterns $\bm{\mu}$ and the margins $\bm{q}= (\bm{n},\bm{m})$. It is concave; its shape depends on the distributions $(\mathbb{P}_{x})$ and $(\mathbb{Q}_{y})$ of the unobserved heterogeneity terms $\bm{\varepsilon}$ and $\bm{\eta}$.
We assume that the analyst observes a random sample of size $N$ from a large population of households. By simple counting (possibly using sampling weights), she obtains estimators of the matching patterns $\hat{\mu}_{xy}, \hat{\mu}_{x0}$, and $\hat{\mu}_{0y}$, as well as the margins:
and a consistent estimator $\Sigma_{\hat{\bm{\mu}}}$ of their asymptotic variance-covariance matrix, given by $$ \Sigma_{\hat{\bm{\mu}}} = \text{diag}(\hat{\bm{\mu}} )-\hat{\bm{\mu}} \hat{\bm{\mu}}^\top. $$
Recall that we have assumed that each $\mathbb{P}_x$ (resp.\ each $\mathbb{Q} _y$) has full support on $\mathrm{I\kern-.17em R}^{Y+1}$ (resp.\ $\mathrm{I \kern-.17em R}^{X+1}$). Then all $\mu_{xy}, \mu_{x0}, \mu_{0y}$ must be positive; as a consequence, the $G_x, H_y, G^\ast_x, H^\ast_y$ functions are continuously differentiable everywhere, as is the generalized entropy function $\mathcal{E}$.
\citecomplete{cupid:restud} showed that at the stable matching $\bm{\mu}$, the joint surplus matrix $\bm{\Phi}$ can be obtained by the following simple formula:
These are the first-order conditions of the maximization of the total joint surplus
Suppose that the distributions $\mathbb{P}_x$ and $\mathbb{Q}_y$ are specified up to a parameter vector $\bm{\alpha} \in \mathbb{R}^{d_{\bm{\alpha}}}$, while the joint surplus matrix $\bm{\Phi}$ is specified up to a parameter vector $\bm{\beta} \in \mathbb{R}^{d_{\bm{\beta}}}$. We write the generalized entropy function $\mathcal{E}^{\bm{\alpha}}$ and the parameterized surplus vector $\bm \Phi^{\bm \beta}$. Then one can use (ref) as the basis for a minimum distance estimator\footnote{Note that in general, one should choose $d_{\bm{\alpha}}+d_{\bm{\beta}} \leq X \times Y$ to ensure identification.}. That is, we write a mixed hypothesis as
stacking all $X\times Y$ conditions in (ref) in a vector $\bm{D}^{ \bm{\lambda}}$.
We choose $\hat{\bm{\lambda}}$ to minimize $\lVert\bm{D}^{\bm{\lambda}}(\bm{ \hat{\mu}}, \bm{\hat{q}})\rVert^2_{\bm{S}}$ for some positive definite $(X\times Y, X\times Y)$ matrix $\bm{S}$. By the general theory of minimum distance estimators, we know that this yields a consistent estimator of $\bm{\lambda}$ if the model is well specified, and that if we choose $\bm{S}= \bm{\hat{\Omega}}^{-1}$ where $\bm{\hat{\Omega}}$ consistently estimates $V \bm{D}^{\bm{\lambda}}(\bm{\hat{\mu}},\bm{\hat{q}})$ (and can be obtained by the delta method), the minimum distance estimator will reach its efficiency bound. Further, if the model is well specified and the choice of $\bm{S}$ is the efficient one, the minimized value of the squared norm follows a $\chi^2$ of degree $X\times Y-d_{\bm{\alpha}}-d_{\bm{\beta}}$. Note that this optimization problem is not a convex optimization problem in general.
Minimum-distance estimation is a particularly appealing strategy if both the derivatives of the generalized entropy function $\mathcal{E}^{\bm{\alpha}}$ and the surplus matrix $\bm{\Phi}^{\bm{\beta}}$ are linear in the parameters:
and
for some vectors of basis functions $\bm{e}(\bm{\mu}, \bm{q})$ and $\bm{\phi} $. Then
is linear in the parameters $\bm{\lambda}$. (Recall that, for every $(x,y) \in \mathcal X \times \mathcal Y$, the vector $ e_{xy}(\bm{\mu}, \bm{q})$ is of size $d_{\bm{\alpha}}$, and $\bm {\phi}_{xy}$ is of size $d_{\bm \beta}$.)
These two conditions call for several remarks. Condition (ref) is a natural choice for a flexible specification. Condition (ref) trivially holds in models where the $\mathbb{P}_x$ and $ \mathbb{Q}_y$ are parameter-free, like the ubiquitous \citecomplete{choo-siow:06} specification. As we will see, it holds in several other leading examples. Note also that the parameter-free part $\bm{e}^0$ is necessary in order to normalize the scale of the error terms, which is otherwise not identified in this discrete-choice model.
Under conditions (ref) and (ref), the minimum distance estimator can be implemented by linear least-squares. Let $ \bm{\hat{F}}$ denote the $(X\times Y, d_{\bm{\alpha}}+d_{\bm{\beta}})$ matrix that stacks $ \bm{e}(\bm{\hat{\mu}},\bm{\hat{q}})$ and $\bm{\phi}$ vertically, so that $ \bm{D}^{\lambda}(\bm{\hat{\mu}},\bm{\hat{r}})=\bm{\hat{e}}^0+\bm{\hat{F}} \bm{\lambda}$, where $\bm{\hat{e}}^0 = \bm{e}^0(\bm {\hat{\mu}},\bm{\hat{q}} ) $ Then for any choice of $\bm{S}$, the minimum distance estimator $\bm{\hat{\lambda}}$ solves the linear system
Since $\bm{\hat{e}}^0$ and $\bm{\hat{F}}$ are functions of $(\bm{\hat{\mu}}, \bm{\hat{q}})$, the variance $\bm{\hat{\Omega}}(\bm{\lambda})$ of $ \bm{\hat{D}}^{\bm{\lambda}}$ can be computed from $\bm{\hat{V}}$ using the delta method. Again, taking $\bm{S}$ to be the inverse of $\bm{\hat{\Omega}}( \bm{\hat{\lambda}})$ is the efficient choice. This procedure is summarized in Box \hyperref[min-distance-estimation-full]{1}.
If the distributions $(\mathbb{P}_x)$ and $(\mathbb{Q}_y)$ are parameter-free, the matrix $\bm{\Omega}^\ast$ does not depend on $ \bm{\lambda}$ any more, and $\bm{\hat{F}}$ is simply the matrix $\bm{\phi}$. The estimators of $\bm{\lambda}=\bm{\beta}$ can be obtained following the procedure described in Box \hyperref[min-distance-estimator]{2}.
Note that since $\bm{\phi}$ is non-random, the variance of $\bm{\hat{D}}^{ \bm{\lambda}}$ is the variance of the derivative of the generalized entropy. Step 2 therefore requires evaluating the second derivatives of the generalized entropy $\mathcal{E}$: by the delta method,
It is easy to see from the definition in (ref) that the first derivative of $\mathcal{E}$ with respect to $\mu_{xy}$ only depends on the conditional matching patterns $\mu_{\cdot\vert x}=(\mu_{x1}/n_x, \ldots, \mu_{xY}/n_x)$ of men of type $x$, and on those of women of type $y$. As a consequence, the Hessians of $\mathcal{E}$ are very sparse and are often easy to evaluate.
We start with two examples for which the generalized entropy and its derivatives are available in closed form; in both cases, the derivatives are linear in the parameters $\bm{\alpha}$. In our third example, the calculation requires finding the fixed point of a contraction, in a way that is familiar from empirical industrial organization.
Let us start with an easy extension of the \citecomplete{choo-siow:06} logit model: the distributions $\mathbb{P}_x$ and $\mathbb{Q}_y$ are type I-EV iid vectors with unknown scale factors $\sigma_x$ and $\tau_y$ respectively. Then $\bm{\alpha}=(\bm{\sigma},\bm{\tau})$ and the derivatives of the generalized entropy function are linear in $\bm{\alpha}$:
where $\mu_{x0}=n_x-\sum_{y\in \mathcal{Y}} \mu_{xy}$ and $ \mu_{0y}=m_y-\sum_{x\in \mathcal{X}} \mu_{xy}$. The second derivatives of the generalized entropy take a very simple form:
and
Scale normalization is done by fixing the value of one of the parameters in $ \bm{\alpha}$. The choo-siow:06 homoskedastic model obtains when all $ \sigma_x$ and $\tau_y$ equal one; a gender-heteroskedastic model would have all $\sigma_x$ equal to one and all $\tau_y$ equal to an unknown $\tau$. \citecomplete{csw:17} applied a minimum distance estimator to the homoskedastic and heteroskedastic logit models.
Consider a two-layer nested logit model. Take men of type $x$ first. Alternative $0$ (singlehood) is obviously special; we put it alone in its nest. Each other nest $n\in \mathcal{N}_x$ contains alternatives $y\in \mathcal{Y}_n$. The correlation of alternatives within nest $n$ is proxied by $1-\left(\rho^x_{n}\right)^2$ (with $\rho^x_{0}=1$ for the nest made of alternative $0$). Similarly, for women of type $y$, alternative $0$ is in a nest by itself with parameter $\delta^y_0=1$ and alternatives $x\in \mathcal{ X}_{n^\prime}$ are in a nest $n\in \mathcal{N}^\prime_y$ with parameter $ \delta^y_{n^\prime}$. We collect the parameters $\bm{\rho}$ and $\bm{\delta}$ into $\bm{\alpha}$.
The formul\ae\ in Example 2.1 of \citecomplete{cupid:restud} imply that if $y$ is in nest $n \in \mathcal{N}_x$ and $x$ is in nest $n^\prime \in \mathcal{N}_y$, then
where we defined $\mu_{xn}=\sum_{t\in \mathcal{Y}_n} \mu_{xt}$ and $ \mu_{n^\prime y}=\sum_{z\in \mathcal{X}_{n^\prime}} \mu_{zy}$. Once again, this is linear in the parameters $\bm{\alpha}$; it remains linear if we impose constraints on the nests (for instance, that $\mathcal{N}_x$ is the same for all types $x$) and/or linear constraints on the $\bm{\rho}$ parameters (for instance, that $\rho_{xn}$ only depends on $n$).
Let us now describe a random coefficient logit model. Consider a man $i$ of type $x$, endowed with preferences $\bm{e}_i$ over a set of $d$ observable characteristics $\bm{Z}$ of potential partners. We add an idiosyncratic shock $\bm{\zeta}_i$ that is distributed as a standard iid type I extreme value vector over $\mathrm{I\kern-.17em R}^{Y+1}$, independently of $\bm{e}_i $, and a scale factor $s>0$:
or in matrix form: $\bm{\varepsilon}=\bm{Z} \bm{e}+s\bm{\zeta}$. This specification is standard in empirical IO. In \citecomplete{blp:1995}: the covariates in $\bm{Z}$ stand for the observed characteristics of the products; the $ \bm{e}$ are individual valuations of these characteristics, and the $\bm{\zeta}$ are idiosyncratic shocks.
Let individual preferences $\bm{e}$ of men of type $x$ have distribution $ \mathbb{P}^e_{x}$. We will seek to estimate the parameters $\bm{\beta}$ of the joint surplus, the scale factor $s$, and the parameters of the distributions $\mathbb{P}^e_{x}$. We collect $s$ and the parameters of $ \mathbb{P}^e_{x}$ in a vector $\bm{\alpha}$.
To compute the derivative of the generalized entropy function, we recall from \citecomplete{cupid:restud} that
By the envelope theorem, the derivative of $G^{\ast}_x(\bm{\nu};\bm{\alpha})$ with respect to $\bm{\nu}$ is the vector $\bm{U}$ that solves the system
This is exactly isomorphic to the inversion problem in \citecomplete{blp:1995}, with the unknown $\bm{U}$ standing for the product effects and $\bm{\nu}$ playing the role of the product market shares. After replacing $\bm{\nu}$ with the observed $\bm{\mu}_{x\cdot}/n_x$, the system can be solved by any of the algorithms that are standard in this literature. The solution gives row $x$ of the matrix $\bm{U}$. Proceeding in the same way for other types of men, and solving for $\bm{V}$ for women, gives the derivatives of the generalized entropy function:
The limit case $s=0$ yields the pure characteristics model of \citecomplete{berrypakes:07}. Then the system to be solved for row $x$ of $\bm{U}$ is
If each $\bm{Z}_t$ is a scalar, the inequalities boil down to
and the system of equations to be solved for $\bm{U}$ is
Now take the generalized entropy function $\mathcal{E}$ as known/assumed; and assume that the joint surplus vector $\bm{\Phi} \in \mathbb{R}^ { \mathcal{X} \times \mathcal{Y} } $ is semilinear: $\bm{\Phi}^{\bm{\beta}}= \bm{\phi} \bm{\beta}$, where $\bm{\beta}$ is a vector of dimension $d_{\bm{\beta}}$ and $ \phi$ is a $| \mathcal{X}| |\mathcal{Y}| \times d_{\bm{\beta}}$ matrix. \citecomplete{cupid:restud} introduced a moment-matching procedure that gives a consistent estimator of the parameter vector $\bm{\beta}$ if the model is well-specified. The moment matching estimator\/ equalizes the observed and simulated comoments, that is the expectations of the basis functions $\bm{\phi}$ under the observed and simulated matching patterns:
where $\bm{\mu}^{\bm{\beta}}$ denotes the stable matching patterns for the parameter vector $\bm{\beta}$. As explained in \citecomplete{cupid:restud}, these are the first-order conditions of the following maximization problem:
where $\mathcal{W}(\bm{\beta},\bm{q})=\max_{\bm{\mu}}\left(\bm{\mu} \bm{\Phi} ^{\bm{\beta}} +\mathcal{E}(\bm{\mu},\bm{q})\right)$ is the value of the total joint surplus. With a semilinear specification for $\bm{\Phi}^{ \bm{\beta}}$, both of these problems are globally convex.
We now show that in the specific (but popular) case of the \citecomplete{choo-siow:06} model, moment matching can be reformulated as a generalized linear model, and estimated by a Poisson regression with two-sided fixed effects.
Define $\mathcal{A=X\times Y\cup X\times }\left\{ 0\right\} \cup \left\{ 0\right\} \times \mathcal{Y}$ the set of possible marital arrangements. Define a vector $w \in \mathbb{R}^\mathcal{A}$ by $w_{xy}=2$ if $x\in \mathcal{X} \text{ and } y\in \mathcal{Y}$ and $w_{xy} = 1$ if $x = 0$ or if $y=0$, so that $w_{xy}$ is the size of household $xy$, namely 2 if matched, 1 if single. The following theorem summarizes our results.
The proof of Theorem (ref) is given in Appendix (ref). The result is very useful in that it allows for inference on $\bm{\beta},\bm{u}$ and $\bm{v}$ in semilinear logit models with standard statistical packages such as glm in R, or scikit-learn in Python. Note that like \citecomplete{sstgravity:06} in the international trade literature, we end up fitting a Poisson regression to a model that is definitely not generated by a Poisson count process. The motivation is different, however. They start from a semiparametric model of the gravity equation and use the robustness of the Poisson pseudo-maximum likelihood estimator. We start from a more complex, fully specified structural model and we show that a semiparametric estimator (moment-matching) is numerically equivalent to the maximum likelihood estimator of a Poisson model.
In the sequel we will denote $\bm{I}_m$ the $(m,m)$ identity matrix; $\bm{p} _{(m,n)}$ the $(m,n)$ matrix whose elements all equal $p$; and $\bm{p} _m\equiv \bm{p}_{(m,1)}$. Also, we say that we stack an $(X,Y)$ matrix in “row-major order” when we create a vector of $X\times Y$ elements whose first $Y$ elements are the first row of the matrix, etc.
As a result, we get that:
We coded these two estimation methods in a Python package called cupid_matching that is available from the standard repositories\footnote{ See \url{https://pypi.org/project/cupid-matching/}.}. To test the quality of the estimators, we generated data both from a choo-siow:06 model and from a semilinear nested logit model. We use both the Poisson estimator and the minimum-distance estimator on the former model, and only the minimum-distance estimator of course on the latter.
In both cases, we take $X=Y=20$ and we use $K=8$ basis functions: $1, x, y, x^2, xy, y^2, \mathrm{1\kern-.40em 1}(x\geq y), \max(x-y, 0)$. The true data-generating process has
so that the true $\bm{\beta}$ is $(1.0,0.0,0.0,-0.01,0.02,-0.01,0.5, 0.0)$. This could be interpreted as the joint surplus from marriage as a function of the ages of the husband $x$ and of the wife $y$. It is highest when the partners have the same age; if they don't, it is larger when the husband is the older partner. We use equal numbers of men and women; and we choose vectors $\bm{n}=\bm{m}$ whose elements form a decreasing geometric sequence with rate $0.8$ (there are fewer individuals available for marriages at higher ages).
The semilinear logit model is entirely described above. We use the IPFP algorithm described in Section 4.2 of \citecomplete{cupid:restud} to solve for the stable matching patterns $\bm{\mu}$ for the margins $\bm{n}$ and $\bm{m}$. To generate a sample, we draw randomly $N=10,000$ households from the multinomial probability distribution generated by $\bm{\mu}$. We generated $S=1,000$ such samples. We used minimum distance estimation and Poisson GLM on each sample. While the minimum distance estimator only uses a linear regression, the Poisson GLM method uses numerical optimization under the hood. In our simulations using the sklearn Python package, the algorithm went astray on 50 of our 1,000 samples, mostly because of overflow errors. We discarded these samples from our analysis.
As Figure (ref) shows, on the remaining 950 samples the two estimators perform about equally well. Both estimators exploit the same $X\times Y$ moment conditions
and both minimize a quadratic form of these conditions. The difference is in the weighting matrix. We saw in Section (ref) that the minimum-distance estimator uses the variance of the derivative of the entropy at the observed matching. On the other hand, the Poisson estimator uses the Hessian of the entropy at the current parameter values. While the two estimators are quite close in our simulations, one can imagine situations in which the divergence would be larger.
For both men and women, we defined three nests that consist of $\{0\}$, $\{1,\ldots,10\}$, and $\{11,\ldots,20\}$. We take the true nest parameters to be all equal to $0.5$ (that is, $\rho_n^x=\delta_{n^ \prime}^y=0.5$ for all $n,n^\prime: x\in n$ and $y\in n^\prime$).
To generate samples from the nested logit model, we proceed as with the logit model. The only difference is that setting up the system to be solved for equilibrium requires a bit more work. We describe our IPFP algorithm in Appendix (ref).
The minimum distance estimator converges fast on all samples. However, we found that a sample size of 10,000 households was much too small to get reliable estimates of the parameters. Figure (ref) gives the distribution of the estimates of the four nest parameters (first two rows) and the eight coefficients of the bases for larger sample sizes: respectively $N=100,000$ and $N=1,000,000$. There is a clear downwards bias on the nest parameters $\rho$ and $\delta$ when $N=100,000$, to the point that some estimates are negative. Some of the coefficients of the bases are also badly estimated. With $N=1,000,000$, the minimum distance estimator performs much better.
Each of the two methods we presented here has its pros and cons.
The minimum-distance estimator applies to all separable models; it is most convenient in semilinear models. To achieve maximum efficiency, and to test the specification, one needs to evaluate the second derivatives of the entropy with respect to the matching patterns. This may be difficult. In addition, the data often contains zero cells---some $\hat{\mu}_{xy}$ may be zero. Then the corresponding equation in (ref) is only an inequality and it must be dropped from the system of estimating equations. An alternative is to add a small positive number $\delta$ to each $\hat{\mu} _{xy}$, to increase the margins $\hat{n}_x$ and $\hat{m}_y$ accordingly, and to estimate on this adjusted data.
The Poisson regression estimator only applies to semilinear \citecomplete{choo-siow:06} models. It is appealing in its simplicity of use, as one can rely on standard statistical packages. It is also more robust to zero cells: nothing in Section (ref) relied on taking derivatives with respect to $\bm{\mu}$ at the observed matching patterns.
Our simulations suggest that it takes large sample sizes to get reliable estimates of distributional parameters (our $\bm{\alpha}$). In labor markets or in marriage markets, large samples are readily available. When they are not (as with matching between firms), it may be better to stick to the \citecomplete{choo-siow:06} specification. Fortunately, the simulations reported in \citecomplete {cns:19} are encouraging as to its robustness.
\printbibliography