Estimation of matching models

Alfred Galichon (NYU & Sciences Po)

'math+econ+code' masterclass on optimal transport and economic applications

With python code examples

© 2018-2022 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.

If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass on optimal transport and economic applications, January 2022. https://github.com/math-econ-code/mec_optim

Learning objectives

• Matching with unobserved heterogeneities

• Estimation of matching models

• Generalized linear models

References

[B] Becker (1973). 'A Theory of Marriage: Part 1.' Journal of Political Economy.

[CS] Choo and Siow (2006). 'Who Marries Whom and Why'. Journal of Political Economy.

[MN] McCullagh and Nelder (1989). Generalized Linear Models, Second Edition. Chapman and Hall/CRC.

[COQ] Chiappori, Oreffice and Quintana-Domeque (2012). 'Fatter Attraction: Anthropometric and Socioeconomic Matching on the Marriage Market'. Journal of Political Economy.

[CSW] Chiappori, Salanié, and Weiss (2017). 'Partner Choice and the Marital College Premium'. American Economic Review.

[DG] Dupuy and Galichon (2014). 'Personality traits and the marriage market'. Journal of Political Economy.

[GS] Galichon and Salanié (2020). 'Cupid's Invisible Hand: Social Surplus and Identification in Matching Models'. Preprint (first version 2011).

Motivation: models of matching since Gary Becker

• In the footsteps of Becker, empirical studies on the marriage market had long been focused on one-dimensional models, which assumes that a single index is enough to capture the interactions on the marriage market, and positive assortative matching (PAM), which predicts that the matching equilibrium will tend to match the agents with higher indices with each other.

• However, it is desirable to move beyond PAM:

  • PAM is always loosely true, never precisely
  • there are often many observed characteristics, and it is not always the case that the sorting can be captured by a single-dimensional model
  • PAM is a theoretical prediction stemming from assumptions of supermodularity of the surplus function which do not necessarly hold
  • optimal transport provide tools to study multidimensional models

• However, any model of matching based on (unregularized) optimal transport will not be exploitable because it will generate far too strong predictions, namely that some matchings will never hold. This is rather counterfactual: in the data, one observes virtually any combination of type.

• Hence, need to regularize the matching model, and we shall do so by introducing unobserved heterogeneity. The model so obtained will be exploitable for estimation and identification purposes. The first such model (with transfers) is the model by [CS]. We shall see a generalization of this model by [GS] (2015).

Loading our libraries

We start with loading the libraries we will need. They are rather standard.

# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here

import pandas as pd
import numpy as np
from scipy import optimize
import scipy.sparse as spr
import gurobipy as grb
from sklearn import linear_model

A look at our data

• Our data are Choo and Siow's original data. Choo and Siow wanted to study the impact of the legalization of abortion by the Roe vs. Wade decision by the supreme court on the 'value of marriage'. Roe vs. Wade decreased the role of marriage in covering out-of-the-wedlocks pregnancies ('shotgun weddings').

• The decision did, however, not make a change uniformly in the United States as a number of states had already legalized abortion (reform states). Choo and Siow thus offer a diffs-in-diffs approach in order to compute the change in the value of marriage.

• Choo and Siow's data are thus made of the marriages between men and women in reformed states (R) vs nonreformed states (NR), in 1972 and in 1982. One should expect to see a higher drop in marriage value in NR states.

thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/marriage-ChooSiow/'

n_singles = pd.read_csv(thepath+'n_singles.txt', sep='\t', header = None)
marr = pd.read_csv(thepath+'marr.txt', sep='\t', header = None)
navail = pd.read_csv(thepath+'n_avail.txt', sep='\t', header = None)

The data used by Choo and Siow is census data on marriages between age categories, from age 16 (row/column 0) to age 75 (row/age 59). It is thus 60x60 tables:

marr.head()

	0	1	2	3	4	5	6	7	8	9	...	50	51	52	53	54	55	56	57	58	59
0	22704	10954	3932	1550	672	414	190	114	78	64	...	0	0	0	0	0	0	0	0	0	0
1	40266	38980	16368	5714	2116	1101	691	427	260	127	...	0	0	0	0	0	0	0	0	0	0
2	39219	49753	40315	18001	5401	2291	1429	856	611	338	...	0	0	0	0	0	0	0	0	0	0
3	29500	44788	47212	39429	16730	6270	3211	1910	762	438	...	0	0	0	0	0	0	0	0	0	0
4	18952	32860	40103	43488	36265	15318	6076	2872	1552	1246	...	0	0	0	0	0	0	0	0	0	0

5 rows × 60 columns

The data also includes the number of single individuals per age category:

n_singles.transpose().head()

	0	1	2	3	4	5	6	7	8	9	...	50	51	52	53	54	55	56	57	58	59
0	1010132	907226	772448	597919	454792	341370	315797	253618	171614	152228	...	68383	69177	61358	70752	69112	60909	60749	65302	62922	61117
1	790793	671332	602900	492020	404882	313755	255414	195418	141283	128787	...	197456	201183	191519	218031	219157	199926	200052	202146	210315	202775

2 rows × 60 columns

The following loads the data and rescales them appropriately:

def prepareChooSiow_data(n_singles,marr,nbCateg = 25):
    μhat_x0 = np.array(n_singles[0].iloc[0:nbCateg])
    μhat_0y = np.array(n_singles[1].iloc[0:nbCateg])
    μhat_xy = np.array(marr.iloc[0:nbCateg:,0:nbCateg])
    Nhat = 2 * μhat_xy.sum() + μhat_x0.sum() + μhat_0y.sum()    
    μhat_a = np.concatenate([μhat_xy.flatten(),μhat_x0,μhat_0y]) / Nhat # rescale the data so that the total number of individual is one

    return μhat_a, Nhat,nbCateg,nbCateg

cs_μhat_a, cs_Nhat,cs_nbx,cs_nby = prepareChooSiow_data(n_singles,marr)

print('Total number of individuals in the sample = '+str(cs_Nhat)+'.')

Total number of individuals in the sample = 14885023.

Building the model

The analysis here follows [GS], who build on the logit model by [CS].

• Consider a heterosexual marriage matching market. The set of types (observable characteristics) is $\mathcal{X}$ for men, and $\mathcal{Y}$ for women. There are $n_{x}$ men of type $x$, and $m_{y}$ women of type $y$.

• Assume that if a man $i\in\mathcal{I}$ of type $x_{i}$ and a woman $j\in\mathcal{J}$ of type $y_{j}$ match, they get respective utilities \begin{align*} & \alpha_{x_{i}y_{j}}+w_{ij}+\varepsilon_{iy_{j}}\\ & \gamma_{x_{i}y_{j}}-w_{ij}+\eta_{x_{i}j} \end{align*} where $w_{ij}$ is the transfer from $i$ to $j$. If they remain single $i$ and $j$ get respectively $\varepsilon_{i0}$ and $\eta_{0j}$.

• The random utility vectors $\left( \varepsilon_{y}\right) $ and $\left( \eta_{x}\right) $ are drawn from probability distributions $\mathbf{P}_{x}$ and $\mathbf{Q}_{y}$, respectively. In the sequel we shall work with a finite number of agents of each type, and then we'll investigate the limit of these results.

• The fact that the preferences for the other side of the market terms $\varepsilon_{iy_{j}}$ and $\eta_{x_{i}}$ do not vary within observed types is a very important implicit assumption called separability.

Optimal matching

The matching surplus between $i$ and $j$ is therefore $$\tilde{\Phi}_{ij}=\Phi_{x_{i}y_{j}}+\varepsilon_{iy_{j}}+\eta_{x_{i}j}$$ where $\Phi_{xy}=\alpha_{xy}+\gamma_{xy}$. The value of optimal matching is thus, under its dual form, \begin{align*}\min_{u_{i},v_{j}} & \sum_{i\in\mathcal{I}}u_{i}+\sum_{j\in\mathcal{J}} v_{j}\\ s.t.~ & u_{i}+v_{j}\geq\Phi_{x_{i}y_{j}}+\varepsilon_{iy_{j}}+\eta_{x_{i}j}\\ & u_{i}\geq\varepsilon_{i0}\\ & v_{j}\geq\eta_{j0} \end{align*}

Written like this, the lp has $\left\vert \mathcal{I}\right\vert +\left\vert \mathcal{J}\right\vert $ variables and $\left\vert \mathcal{I} \right\vert \times\left\vert \mathcal{J}\right\vert +\left\vert \mathcal{I} \right\vert +\left\vert \mathcal{J}\right\vert $ constraints. Assuming that there are $K$ individuals per type for each type, this is $K\left( \left\vert \mathcal{X}\right\vert +\left\vert \mathcal{Y}\right\vert \right) $ variables and $K^{2}\left( \left\vert \mathcal{X}\right\vert \times\left\vert \mathcal{Y}\right\vert \right) +K\left( \left\vert \mathcal{X}\right\vert +\left\vert \mathcal{Y}\right\vert \right) $ constraints.

The number of constraints is quadratic with respect to $K$. Fortunately, a little thinking about the implications of separability will help us reduce this complexity.

A property of equilibrium

We have:

---

Lemma. Consider the set $\mathcal{I}_{xy}$ of men of type $x$ matched to women of type $y$ at equilibrium. If $\mathcal{I}_{xy}$ is nonempty, then $u_{i}-\varepsilon_{iy}$ is a constant across $\mathcal{I}_{xy}$.

---

Proof. For $i\in\mathcal{I}$ such that $x_{i}=x$, \begin{align*} u_{i} & =\max_{j\in\mathcal{J}}\left\{ \tilde{\Phi}_{ij}-v_{j}% ,\varepsilon_{i0}\right\} \\ & =\max_{y\in\mathcal{Y}}\left\{ U_{xy}+\varepsilon_{iy},\varepsilon _{i0}\right\} \end{align*} where $U_{xy}=\max_{j:y_{j}=y}\left\{ \Phi_{xy}+\eta_{x_{i}j}-v_{j}\right\} $, thus $u_{i}\geq U_{xy}+\varepsilon_{iy}$ with equality on $\mathcal{I}% _{xy}$. With similar notations, $v_{j}\geq V_{xy}+\eta_{xj}$ with equality on $\mathcal{J}_{xy}$. As a result, if $\mathcal{I}_{xy}$ is nonempty, then $U_{xy}+V_{xy}=\Phi_{xy}$ and $\forall i\in\mathcal{I}_{xy},$ $u_{i}% =U_{xy}+\varepsilon_{iy}$.

A simplification

In the sequel, we shall see that adding an auxiliary variable to the previous lp will lead to decreasing the computational complexity of the problem.

Observe that the first set of constraints is reexpressed by saying that, for every $x\in\mathcal{X}$, $y\in\mathcal{Y}$, $$ \min_{i:x_{i}=x}\left\{ u_{i}-\varepsilon_{iy}\right\} +\min_{j:y_{j} =y}\left\{ v_{j}-\eta_{xj}\right\} \geq\Phi_{xy}. $$

Hence, letting $U_{xy}=\min_{i:x_{i}=x}\left\{ u_{i}-\varepsilon _{iy}\right\} $ and $V_{xy}=\min_{j:y_{j}=y}\left\{ v_{j}-\eta_{xj}\right\} $, a solution of the previous lp should satisfy $$ u_{i}=\max_{y\in\mathcal{Y}}\left\{ U_{xy}+\varepsilon_{iy},\varepsilon _{i0}\right\} \text{ and }v_{j}=\max_{x\in\mathcal{X}}\left\{ V_{xy} +\varepsilon_{xj},\varepsilon_{0j}\right\} . $$

The problem rewrites as \begin{align} \min_{u_{i},v_{j},U_{xy},V_{xy}} & \sum_{i\in\mathcal{I}}u_{i}+\sum _{j\in\mathcal{J}}v_{j}\label{simplifiedDual}\\ s.t.~ & U_{xy}+V_{xy}\geq\Phi_{xy}~\left[ \mu_{xy}\geq0\right] \nonumber\\ & u_{i}\geq U_{x_{i}y}+\varepsilon_{iy_{j}}~\left[ \mu_{iy}\right] \nonumber\\ & v_{j}\geq V_{xy_{j}}+\eta_{x_{i}j}~\left[ \mu_{xi}\right] \nonumber\\ & u_{i}\geq\varepsilon_{i0}~\left[ \mu_{i0}\right] \nonumber\\ & v_{j}\geq\eta_{j0}~\left[ \mu_{0x}\right] \nonumber \end{align}

This problem has $K\left( \left\vert \mathcal{X}\right\vert +\left\vert \mathcal{Y}\right\vert \right) +\left\vert \mathcal{X}\right\vert \times\left\vert \mathcal{Y}\right\vert $ variables and $\left( \left\vert \mathcal{X}\right\vert \times\left\vert \mathcal{Y}\right\vert \right) +K\left( 2\left\vert \mathcal{X}\right\vert \times\left\vert \mathcal{Y} \right\vert +\left\vert \mathcal{X}\right\vert +\left\vert \mathcal{Y} \right\vert \right) $ constraints.

The number of constraint is now linear with respect to $K$.

Consequences

1. Lagrange multipliers:

• The Lagrange multiplier $\mu_{xy}$ is interpreted as the number of matchings between types $x$ and $y$.

• The Lagrange multiplier $\mu_{iy}$ ($y\in\mathcal{Y}_{0}$) is interpreted as a 0-1 indicator that man $i$ chooses a type $y$

• The Lagrange multiplier $\mu_{xj}$ ($x\in\mathcal{X}_{0}$) is interpreted as a 0-1 indicator that woman $j$ chooses a\ type $x$

2. Utilities:

• Man $i$ solves a discrete choice problem $u_{i}=\max_{y\in\mathcal{Y} }\left\{ U_{xy}+\varepsilon_{iy},\varepsilon_{i0}\right\} $

• Woman $j$ solve a discrete choice problem $v_{j}=\max_{x\in\mathcal{X} }\left\{ V_{xy}+\eta_{xj},\eta_{0j}\right\} .$ \end{itemize}

• $U_{xy}$ and $V_{xy}$ are related by $U_{xy}+V_{xy}\geq\Phi_{xy}$ with equality if $\mu_{xy}>0$.

Large market limit

Now look at the limit of previous markets when the number of market participants gets large, holding fixed the frequency of each types.

In the large population limit $n_{x}$ and $m_{y}$ are now interpreted as the mass distribution of respective types $x$ and $y$.

We shall from now on assume that $\mathbf{P}_{x}$ and $\mathbf{Q}_{y}$, the distributions of random utility vectors $\left( \varepsilon_{y}\right) $ and $\left( \eta_{x}\right) $, have a density with full support. This will ensure that the Emax operators associated with the choice problems of the men and the women respectively \begin{align*} & G_x(U_{x.}) = \mathbb{E}_\mathbf{P} \left[\max_{y\in\mathcal{Y} }\left\{ U_{xy}+\varepsilon_{iy},\varepsilon_{i0}\right\} \right]\text{, and }\\ & H_y(V_{.y}) = \mathbb{E}_\mathbf{Q} \left[\max_{x\in\mathcal{X} }\left\{ V_{xy}+\eta_{xj},\eta_{0j}\right\} \right],\end{align*}as well as the corresponding entropies of choice $G_x^{\ast}$ and $H_y^{\ast}$ are continuously differentiable.

Under these assumptions, the problem becomes \begin{align*} \min_{U,V} ~& G\left( U\right) +H\left( V\right) \\ s.t.~ & U_{xy}+V_{xy}\geq\Phi_{xy}~\left[ \mu_{xy}\right] \end{align*} where \begin{align*} G\left( U\right) & =\sum_{x\in\mathcal{X}}n_{x}\mathbb{E}_{\mathbf{P}% }\left[ \max_{y\in\mathcal{Y}}\left\{ U_{xy}+\varepsilon_{iy},\varepsilon _{i0}\right\} \right] \\ H\left( V\right) & =\sum_{y\in\mathcal{Y}}m_{y}\mathbb{E}_{\mathbf{Q}% }\left[ \max_{x\in\mathcal{X}}\left\{ V_{xy}+\eta_{xj},\eta_{0j}\right\} \right] \end{align*}

By first order conditions, $$ \frac{\partial G\left( U\right) }{\partial U_{xy}}=\mu_{xy}=\frac{\partial H\left( V\right) }{\partial V_{xy}}. $$ and $\mu_{xy}>0$ for every $x\in\mathcal{X}$ and $y\in\mathcal{Y}$.

Social planner's problem

The primal problem corresponding the problem above is $$ \max_{\mu_{xy}\geq0}\sum_{\substack{x\in\mathcal{X}\\y\in\mathcal{Y}}}\mu _{xy}\Phi_{xy}-\mathcal{E}\left( \mu\right) $$ where $$ \mathcal{E}\left( \mu\right) =G^{\ast}\left( \mu\right) +H^{\ast}\left( \mu\right) $$

Recall $G^{\ast}\left( \mu\right) =\max\left\{ \sum_{xy}\mu _{xy}U_{xy}-G\left( U\right) \right\} $ is the Legendre transform of $G$, and similarly for $H^{\ast}$.

Constructing the surplus:

def build_surplus(nbx,nby):
    xs, ys = np.repeat(range(1,nbx+1),nby).reshape(nbx,nby)/nbx, np.repeat(range(1,nby+1),nbx).reshape(nbx,nby).T/nby
    phi1_xy = -((xs-ys)**2).flatten()
    phimat = np.vstack([phi1_xy,
                        -((xs-ys)**2 * ((xs+ys)/2)**2).flatten(),
                         -( (xs-ys)**2 * ((xs+ys-2)/2)**2).flatten(),
                        -( (xs-ys)**2 *(xs+ys-1)**2).flatten()
                       ]).T
    phimat_mean = phimat.mean(axis = 0 , keepdims=True)
    phimat_stdev = phimat.std(axis = 0 , keepdims=True)
    φ_xy_k = np.hstack([np.ones((nbx*nby,1)) ,(phimat - phimat_mean)/phimat_stdev]  )
    _,nbk = φ_xy_k.shape
    return np.vstack( [φ_xy_k ,np.zeros( (nbx,nbk) ),np.zeros( (nby,nbk) )])

Denoting $a$ a generic element of

$\mathcal{A}=\mathcal{X}\times \mathcal{Y}\cup \mathcal{X}\times \left\{ 0\right\} \cup \left\{ 0\right\} \times \mathcal{Y}$ the set of matches, and
$\mathcal{K}=\left\{ 1,...,K\right\} $ the index set of the parameter,
we may view

• $\Phi =\left( \phi _{a}^{k}\right) $ as a $\mathcal{A}\times \mathcal{K}$ matrix, and
• $\hat{\mu}$ as an element of $\mathbb{R}^\mathcal{A}$,

and we define $d_{xy}=2$, $d_{x0}=d_{0y}=1$.

Let's create a class TUlogit that handles this:

class TUlogit:
    def __init__(self,nbx,nby,φ_a_k , μhat_a ):
        self.nba, self.nbk = φ_a_k.shape
        self.nbx, self.nby  = nbx, nby
        self.φ_a_k, self.μhat_a = φ_a_k, μhat_a
        self.d_a = np.append(2*np.ones(nbx*nby),np.ones(nbx+nby))
        μ_x_y = self.μhat_a[0:(nbx*nby)].reshape((nbx,nby))
        self.n_x  = μ_x_y.sum(axis = 1) + self.μhat_a[(nbx*nby):(nbx*nby+nbx)]
        self.m_y  = μ_x_y.sum(axis = 0) + self.μhat_a[-nby:]

        
choo_siow_mkt = TUlogit(cs_nbx,cs_nby,build_surplus(cs_nbx,cs_nby),cs_μhat_a)  

Identification of the matching surplus

---

Theorem. By first order conditions, we get the identifcation formula of $\Phi$: $$ \Phi_{xy}=\frac{\partial G^{\ast}\left( \mu\right) }{\partial\mu_{xy}} +\frac{\partial H^{\ast}\left( \mu\right) }{\partial\mu_{xy}}. $$

---

This means that the surplus function is identified nonparametrically given the matching patterns $\mu$ and assuming a fixed distribution of unobserved heterogeneity.

Hence only the joint surplus $\Phi _{xy}=\alpha_{xy}+\gamma_{xy}$ is identified. However, if the transfers $\hat{w}_{xy}$ are observed too (e.g. wages in labour market), then $U_{xy}=\alpha_{xy}+w_{xy}$ and $V_{xy}=\gamma_{xy}-w_{xy}$, so that $\alpha$ and $\gamma$ are separately identified by $$ \left\{ \begin{array} [c]{c}% \hat{\alpha}_{xy}=\frac{\partial G^{\ast}\left( \mu\right) }{\partial\mu_{xy}}-\hat{w}_{xy}\\ \hat{\gamma}_{xy}=\frac{\partial H^{\ast}\left( \mu\right) }{\partial\mu_{xy}}+\hat{w}_{xy}% \end{array} \right. $$

Choo and Siow's logit model

In Choo and Siow's model [CS], the heterogeneities in tastes are Gubmel, we have $$ \mathcal{E}\left( \mu\right) =2\sum_{\substack{x\in\mathcal{X}% \\y\in\mathcal{Y}}}\mu_{xy}\log\mu_{xy}+\sum_{x\in\mathcal{X}}\mu_{x0}\log \mu_{x0}+\sum_{y\in\mathcal{Y}}\mu_{0y}\log\mu_{0y}. $$ Note that $\mathcal{E}\left( \mu\right) < + \infty$ if and only if $\mu \in\mathcal{M}\left( n,m\right) $.

By first order conditions above, Choo-Siow's TU-logit model implies the following matching function: $$ \mu_{xy}=M_{xy}\left( \mu_{x0},\mu_{0y}\right) :=\sqrt{\mu_{x0}}\sqrt {\mu_{0y}}\exp\left( \frac{\Phi_{xy}}{2}\right) $$

This is a gravity equation of sorts. The full link with gravity equations is explored in the next lecture.

As a result, $\partial\mathcal{E}\left( \mu\right) /\partial\mu _{xy}=2\log\mu_{xy}-\log\mu_{x0}-\log\mu_{0y}$, which implies that $\Phi_{xy}$ is estimated by Choo and Siow's identification formula $$ \hat{\Phi}_{xy}=\log\frac{\hat{\mu}_{xy}^{2}}{\hat{\mu}_{x0}\hat{\mu}_{0y}}. $$

Solving equilibrium in the Choo-Siow model

Write down the equilibrium equations in the TU-logit model: $$ \left\{ \begin{array} [c]{c} \sum_{y\in\mathcal{Y}}\sqrt{\mu_{x0}}\sqrt{\mu_{0y}}\exp\left( \frac {\Phi_{xy}}{2}\right) +\mu_{x0}=n_{x}\\ \sum_{x\in\mathcal{X}}\sqrt{\mu_{x0}}\sqrt{\mu_{0y}}\exp\left( \frac {\Phi_{xy}}{2}\right) +\mu_{0y}=m_{y}% \end{array} \right. $$

Setting $a_{x}=\sqrt{\mu_{x0}}$, $b_{y}=\sqrt{\mu_{0y}}$, and $K_{xy}=\exp\left( \Phi_{xy}/2\right) $, this rewrites as $$ \left\{ \begin{array} [c]{c} \sum_{y\in\mathcal{Y}}K_{xy}a_{x}b_{y}+a_{x}^{2}=n_{x}\\ \sum_{x\in\mathcal{X}}K_{xy}a_{x}b_{y}+b_{y}^{2}=m_{y}% \end{array} \right.$$

which is a variant of the equations previously seen to accomodate unmatched agents.

We can easily adapt the IPFP to this setting. The IPFP will consists in iteratively solving quadratic equations: $$ \left\{ \begin{array} [c]{l} a_{x}^{2t+1}=\sqrt{n_{x}+\left( \sum_{y\in\mathcal{Y}}b_{y}^{2t}% K_{xy}/2\right) ^{2}}-\sum_{y\in\mathcal{Y}}b_{y}^{2t}K_{xy}/2\\ b_{y}^{2t+2}=\sqrt{m_{y}+\left( \sum_{x\in\mathcal{X}}a_{x}^{2t+1}% K_{xy}/2\right) ^{2}}-\sum_{x\in\mathcal{X}}a_{x}^{2t+1}K_{xy}/2 \end{array} \right. $$

Dual problem

The dual problem is given by $$ \min_{u,v}\left\{ \begin{array} [c]{c}% \sum_{x}n_{x}u_{x}+\sum_{y}m_{y}v_{y}\\ +2\sum_{xy}\sqrt{n_{x}m_{y}}\exp\left( \frac{\Phi_{xy}-u_{x}-v_{y}}{2}\right) \\ +\sum_{x}n_{x}\exp\left( -u_{x}\right) +\sum_{y}m_{y}\exp\left( -v_{y}\right) \end{array} \right\} $$

Remarks.

• This problem is an unconstrained convex optimization problem, so this formulation will be quite useful.

• If $\left( u,v\right) $ is solution, $u_{x}=-\log\mu_{0|x}=-\log\left( \mu_{x0}/n_{x}\right) $ and $v_{y}=-\log\left( \mu _{0|y}\right) $.

• Note that the IPFP algorithm just seen interprets as (blockwise) coordinate descent method in the dual problem.

DIY approach

We first try the "do-it-yourself" approach, following [GS]. Construct the objective function:

def ObjFunc(uvlambda,scale,mkt):
    nbx,nby,nbk=mkt.nbx,mkt.nby,mkt.nbk
    u_x,v_y,λ_k = uvlambda[0:nbx] / scale,uvlambda[nbx:(nbx+nby)]/scale,uvlambda[(nbx+nby):(nbx+nby+nbk)] / scale
    
    Phi_a = mkt.φ_a_k @ λ_k
    mu_xy = np.exp((Phi_a[0:(nbx*nby)].reshape((nbx,nby)) - u_x.reshape((-1,1))-v_y.reshape((1,-1)))/2).flatten()
    mu_x0 = np.exp(Phi_a[(nbx*nby):(nbx*nby+nbx)]-u_x)
    mu_0y = np.exp(Phi_a[-nby:]-v_y)
    
    return (mkt.n_x.dot(u_x)+mkt.m_y.dot(v_y) - (Phi_a * mkt.μhat_a).sum() + 2 * mu_xy.sum() + mu_x0.sum() + mu_0y.sum())

Construct its gradient:

def grad_ObjFunc(uvlambda,scale,mkt):
    nbx,nby,nbk=mkt.nbx,mkt.nby,mkt.nbk
    u_x,v_y,λ_k = uvlambda[0:nbx]/scale,uvlambda[nbx:(nbx+nby)]/scale,uvlambda[(nbx+nby):(nbx+nby+nbk)]/scale
    
    Phi_a = mkt.φ_a_k @ λ_k
    mu_x_y = np.exp((Phi_a[0:(nbx*nby)].reshape((nbx,nby)) - u_x.reshape((-1,1))-v_y.reshape((1,-1)))/2).reshape((nbx,nby))
    mu_x0 = np.exp(Phi_a[(nbx*nby):(nbx*nby+nbx)]-u_x)
    mu_0y = np.exp(Phi_a[-nby:]-v_y)
    mu_a = np.concatenate((mu_x_y.flatten(),mu_x0,mu_0y))
    
    grad_u = mkt.n_x - mu_x_y.sum(axis = 1) - mu_x0
    grad_v = mkt.m_y - mu_x_y.sum(axis = 0) - mu_0y
    grad_lambda = (mu_a-mkt.μhat_a ) @ mkt.φ_a_k
    
    return np.concatenate((grad_u,grad_v,grad_lambda.flatten()))/scale

uvl_0 = np.repeat(0,choo_siow_mkt.nbx+choo_siow_mkt.nby+choo_siow_mkt.nbk)
grad_ObjFunc(uvl_0,1,choo_siow_mkt)

array([-2.59293961e+01, -2.59318975e+01, -2.59374157e+01, -2.59469373e+01,
       -2.59559106e+01, -2.59643718e+01, -2.59682064e+01, -2.59738609e+01,
       -2.59826665e+01, -2.59852835e+01, -2.59877521e+01, -2.59881513e+01,
       -2.59907795e+01, -2.59924915e+01, -2.59927869e+01, -2.59936593e+01,
       -2.59943827e+01, -2.59944988e+01, -2.59951420e+01, -2.59949721e+01,
       -2.59952618e+01, -2.59952660e+01, -2.59954071e+01, -2.59948270e+01,
       -2.59945361e+01, -2.59343738e+01, -2.59385664e+01, -2.59436264e+01,
       -2.59523431e+01, -2.59602771e+01, -2.59693410e+01, -2.59758246e+01,
       -2.59812204e+01, -2.59870677e+01, -2.59887056e+01, -2.59896728e+01,
       -2.59898079e+01, -2.59917911e+01, -2.59924536e+01, -2.59920896e+01,
       -2.59934063e+01, -2.59935473e+01, -2.59938525e+01, -2.59937524e+01,
       -2.59938244e+01, -2.59942826e+01, -2.59937386e+01, -2.59939844e+01,
       -2.59934018e+01, -2.59931882e+01,  6.24885633e+02, -8.02789702e-02,
       -8.05367354e-02, -6.73882331e-02, -5.27153039e-02])

Optimize using scipy.optimize:

def fit_diy(self,scale=1,tol=1e-9):
    uvl_0 = np.repeat(0,self.nbx+self.nby+self.nbk)
    outcome = optimize.minimize(ObjFunc,jac = grad_ObjFunc,args=(scale,self), method = 'CG', x0=uvl_0,options={'gtol':tol,'disp': True})
    return  outcome['x'][0:-self.nbk]/scale,outcome['x'][-self.nbk:]/scale

TUlogit.fit_diy = fit_diy
uv_GS_diy,λ_GS_diy_k = choo_siow_mkt.fit_diy(10000,tol = 1e-8)
print('λ_GS_diy_k = '+str(λ_GS_diy_k))

Optimization terminated successfully.
         Current function value: 5.594425
         Iterations: 37
         Function evaluations: 589
         Gradient evaluations: 589
λ_GS_diy_k = [-10.33712436   1.39174524   1.04488304   1.53388576   0.62883546]

GLM computation

We can reformulate the model as a generalized linear model, see [MN] for a reference on the latter.

Setting
$B= \begin{pmatrix} I_{X}\otimes 1_{Y} & 1_{X}\otimes I_{Y} \\ I_{X} & 0_{X\times Y} \\ 0_{Y\times X} & I_{Y}% \end{pmatrix}% $
and
$C=\left( -B, \Phi \right) ,$
which we code into:

def TUlogit_C(self):
    nbx, nby, nbk = self.nbx, self.nby, self.nbk
    B = spr.vstack([spr.hstack([spr.kron(spr.identity(nbx), np.ones((nby,1))),spr.kron(np.ones((nbx,1)),spr.identity(nby))]),
                    spr.hstack([spr.identity(nbx) , spr.csr_matrix((nbx,nby))]),
                    spr.hstack([spr.csr_matrix((nby,nbx)),spr.identity(nby) ])])
    return spr.hstack([-B,self.φ_a_k])

TUlogit.C = TUlogit_C

Letting $\theta =\left( u,v,\lambda \right) $,
one has that the estimation problem above reformulates as

$\min_{\theta }\left\{ 1^{\top }\delta \exp \left( \delta ^{-1}\left( C\theta \right) \right) -\left( \delta \hat{\mu}\right) ^{\top }\delta ^{-1}C\theta \right\} $

and thus it results from the Poisson regression of
$Y=\hat{\mu}$ on $X=\delta ^{-1}C$ with weights vector $d$.

Poisson coerciveness

Coercivity detector:

$\max_{\mu \geq 0,u\geq 0}u$

s.t. $X^{\top }\mu +0u=b$

$I\mu -1u\geq 0$

def poisson_coerciveness(X,y):
    n,m = X.shape
    b = X.T @ y
    m=grb.Model()
    u = m.addMVar(shape = 1)
    mu = m.addMVar(shape=n)
    m.setObjective(1 * u, grb.GRB.MAXIMIZE)
    m.addConstr(X.T @ mu == b)
    m.addConstr( mu - np.ones((n,1)) @ u >= 0)
    m.optimize()
    if m.status == grb.GRB.Status.OPTIMAL:
        solution = np.array(m.getAttr('x'))
        u = solution[0]
    else:
        u = np.NaN
    return(u)

We use the Poisson regression included in the scikit-learn library, see
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.PoissonRegressor.html

def TUlogit_fit_glm(self,tol = 1e-12,pretest = False):
    clf = linear_model.PoissonRegressor(fit_intercept=False,tol=tol ,verbose=3,alpha=0)
    clf.fit( spr.diags(1/ self.d_a) @ self.C(), self.μhat_a,sample_weight = self.d_a)
    return clf.coef_[0:-self.nbk], clf.coef_[-self.nbk:]

TUlogit.fit_glm = TUlogit_fit_glm

def TUlogit_isCoercive(self):
    return poisson_coerciveness( self.C(), self.μhat_a)

TUlogit.isCoercive = TUlogit_isCoercive


def TUlogit_assess(self,θ):
    C= self.C()
    return np.abs((np.exp(C @ θ) - self.μhat_a).reshape((1,-1)) @ C).max()

TUlogit.assess = TUlogit_assess

uv_GS_GLM, λ_GS_GLM_k = choo_siow_mkt.fit_glm()
print('λ_GS_GLM_k = '+str(λ_GS_GLM_k))
print('λ_GS_diy_k = '+str(λ_GS_diy_k))
print(choo_siow_mkt.assess(np.append(uv_GS_GLM, λ_GS_GLM_k)))
print(choo_siow_mkt.assess(np.append(uv_GS_diy, λ_GS_diy_k)))

λ_GS_GLM_k = [-10.32125662   1.38355341   1.02107659   1.54579903   0.62931836]
λ_GS_diy_k = [-10.33712436   1.39174524   1.04488304   1.53388576   0.62883546]
0.11428192310395208
0.11428184295353003

choo_siow_mkt.isCoercive()

Using license file C:\Users\Alfred\gurobi.lic
Academic license - for non-commercial use only
Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (win64)
Optimize a model with 730 rows, 676 columns and 5775 nonzeros
Model fingerprint: 0x5bce378c
Coefficient statistics:
  Matrix range     [9e-04, 4e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e-03, 1e-01]
Presolve removed 51 rows and 50 columns
Presolve time: 0.01s
Presolved: 679 rows, 626 columns, 5050 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    5.6173242e-03   3.680846e+00   0.000000e+00      0s
     673    2.0886948e-05   0.000000e+00   0.000000e+00      0s

Solved in 673 iterations and 0.08 seconds
Optimal objective  2.088694766e-05

2.0886947663455243e-05