Recovering Latent Variables by Matching

ABSTRACT

We propose an optimal-transport-based matching method to nonparametrically estimate linear models with independent latent variables. The method consists in generating pseudo-observations from the latent variables, so that the Euclidean distance between the model’s predictions and their matched counterparts in the data is minimized. We show that our nonparametric estimator is consistent, and we document that it performs well in simulated data. We apply this method to study the cyclicality of permanent and transitory income shocks in the Panel Study of Income Dynamics. We find that the dispersion of income shocks is approximately acyclical, whereas the skewness of permanent shocks is procyclical. By comparison, we find that the dispersion and skewness of shocks to hourly wages vary little with the business cycle. Supplementary materials for this article are available online.

KEYWORDS:

• Factor models
• Income dynamics
• Latent variables
• Matching
• Nonparametric estimation
• Optimal transport

Supplementary material

In addition to these extensions, in the online appendix we provide the proofs and present some additional simulations.

Acknowledgments

We thank to Colin Mallows, Tincho Almuzara, Alfred Galichon, Jiaying Gu, Kei Hirano, Pierre Jacob, Roger Koenker, Thibaut Lamadon, Guillaume Pouliot, Azeem Shaikh, Tim Vogelsang, Daniel Wilhelm, and audiences at various places for comments. Tincho Almuzara and Beatriz Zamorra provided excellent research assistance.

Notes

1 Codes to implement the estimator are available on the second author’s webpage.

2 When $A\prime A$ is nonsingular and A is known, $\widehat{X}={(A\prime A)}^{-1}A\prime Y$ recovers X exactly. We are interested in situations, such as deconvolution and filtering, where exact recovery of the latent variables is not possible.

3 An interesting possibility would be to jointly estimate A and the distribution of X. Although we do not study it formally, we comment on this possibility in Section 4.2.

4 Other applications in economics include the estimation of the heterogeneous effects of an exogenous binary treatment under the assumption that the potential outcome in the absence of treatment is independent of the gains from treatment (Heckman, Smith, and Clements 1997), and the estimation of the distribution of time-invariant random coefficients of binary treatments in panel data models (Arellano, and Bonhomme 2012).

5 See also Botosaru and Sasaki (2015). Our approach may be used to estimate linear autoregressive specifications of the form ${\eta }_t=\alpha +\rho {\eta }_{t-1}+v_t$, where we estimate $(\alpha ,\rho )$—that is, the matrix A—in a first step. An important application of error components models is to relax independence in repeated measurements models such as Equation (1). This can be done provided T is large enough. Modeling ${\epsilon }_t$ in Equation (1) as a finite-order moving average or autoregressive process with independent innovations preserves the linear independent factor structure of the model (Arellano, and Bonhomme 2012; see also Hu, Moffitt, and Sasaki 2019). In addition, in model (1) Schennach (2013b) pointed out that full independence between the factors is not necessary, and that sub-independence suffices to establish identification.

6 Ben Moshe (2017) showed how to allow for arbitrary subsets of dependent factors, and proposes characteristic-function based estimators.

7 The sample sizes being the same for Y and X₂ is not essential and can easily be relaxed. In a setting where the cdf $F_{X_2}$ is known, one can draw a sample from it, or alternatively work with an integral counterpart to our estimator.

8 Specifically, one could compute $X_{\sigma (i,r),1}+X_{i2}$, with $\sigma (·,1),\dots ,\sigma (·,R)$ being R independent permutations. In that case, π would be a generalized permutation, mapping ${\{1,\dots ,N\}}^R$ to $\{1,\dots ,N\}$.

9 It is common in applications to assume that some of the X_k’s have zero mean while leaving the remaining means unrestricted. For example, in the repeated measurements model, assuming that $\mathbb{E}(X_1)=0$ suffices for identification. Our algorithm can easily be adapted to such cases.

10 If $A=\left\{a_{tk}\right\}$ is unknown and $\widehat{A}=\left\{{\widehat{a}}_{tk}\right\}$ is a consistent estimate of it, then we replace a_tk by ${\widehat{a}}_{tk}$ in Equation (4). We proceed similarly in the algorithm we propose in the next section. Alternatively, one could jointly minimize the objective function on the right-hand side of Equation (4) with respect to both X and $\left\{a_{tk}\right\}$. Here, we do not study the formal properties of such a joint estimation method.

11 Notice that, since π is a permutation, ${\sum }_{i=1}^N{\sum }_{t=1}^T{Y}_{\pi (i),t}^2={\sum }_{i=1}^N{\sum }_{t=1}^T{Y}_{it}^2$ does not depend on π.

12 See, for example, Galichon (2016, chap. 3) on discrete Monge–Kantorovitch problems.

13 An entropic-regularized counterpart to Equation (4) is, for ${ϵ}_N>0$

$\widehat{X}=\underset{X\in {\mathcal{X}}_N}{\text{argmin}}\{\underset{P\in {\mathcal{P}}_N}{\text{min}}{\sum }_{i=1}^N{\sum }_{j=1}^N{P}_{ij}{\sum }_{t=1}^T{(Y_{jt}-{\sum }_{k=1}^K{a}_{tk}{X}_{{\sigma }_k(i),k})}^2\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+{ϵ}_N{\sum }_{i=1}^N{\sum }_{j=1}^N{P}_{ij}(\hspace{0.17em}\text{ln}\hspace{0.17em}(P_{ij})-1)\}.$

14 Strictly speaking, Mallows (2007) redefined ${\widehat{X}}_{i1}^{(s+1)}\equiv {\widehat{X}}_{{\sigma }^{(s)}(i),1}^{(s+1)}$ for all $i=1,\dots ,N$ at the end of step s, and then applies the random permutation ${\sigma }^{(s+1)}$ to the new ${\widehat{X}}^{(s+1)}$ values. This difference with the algorithm outlined here turns out to be immaterial, since the composition of ${\sigma }^{(s+1)}$ and ${\sigma }^{(s)}$ is also a random permutation of $\{1,\dots ,N\}$.

15 It is not necessary for ${\widehat{X}}_k$ to be an exact minimizer of Equation (4). As we show in the proof, it suffices that the value of the objective function at $({\widehat{X}}_1,\dots ,{\widehat{X}}_K)$ be in an ϵ_N-neighborhood of the global minimum, for ϵ_N tending to zero as N tends to infinity.

16 Part (i) ensures that $F_{X_k}^{-1}$ belongs to an $\Vert ·{\Vert }_{1,\infty }$ -ball, which is compact under $\Vert ·{\Vert }_{\infty }$ (Gallant and Nychka 1987). Compactness can be preserved when norms are replaced by weighted norms (e.g., using polynomial or exponential weights); see, for example, (Freyberger and Masten 2019, theor. 7), and the analysis in Newey and Powell (2003).

17 The alternative density estimator ${\tilde{f}}_{X_k}(x)\equiv 1/\nabla {\widehat{H}}_k({\widehat{H}}_k^{-1}(x))$ can be shown to be uniformly consistent for $f_{X_k}$ as N tends to infinity under the same conditions.

18 A simple recommendation for practice is based on a truncated normal distribution. Let ${\widehat{\sigma }}_k$ denote a consistent estimate of the standard deviation of X_k, for example, obtained by covariance-based minimum distance, and let c > 0 be a tuning parameter. Possible penalization constants are: $2.3c{\widehat{\sigma }}_k$ (upper bound on quantile values), $2.5c^{-1}{\widehat{\sigma }}_k$ and $37c{\widehat{\sigma }}_k$ (lower and upper bounds for first derivatives), and $3300c{\widehat{\sigma }}_k$ (upper bound on the second derivatives). When c = 1, these constants are binding when X_k follows a normal truncated at the 99th percentiles. As a default choice one may take c = 2.

19 When implementing the Fourier estimator we enforce the non-negativity and integral constraints ex-post. To select the tuning parameter, we minimize the Monte Carlo MISE of the estimator on a grid of values.

20 Indeed, we have $$\begin{array}{l}\underset{\equiv Y}{\underset{}{\underbrace{\left(\text{}\text{}\begin{array}{c}\Delta Y_1\\ \Delta Y_2\\ \Delta Y_3\\ \mathrm{...}\\ \Delta Y_T\end{array}\text{}\text{}\right)}}}=\underset{\equiv A}{\underset{}{\underbrace{\left(\text{}\begin{array}{cccccccc}1& 0& \mathrm{...}& 0& 1& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0& -1& 1& \mathrm{...}& 0\\ 0& 0& \mathrm{...}& 0& 0& -1& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& 1& 0& 0& \mathrm{...}& -1\end{array}\text{}\right)}}}\\ \underset{\equiv X}{\underset{}{\underbrace{\left(\begin{array}{c}{v}_1-{\epsilon }_0\\ v_2\\ \mathrm{...}\\ v_T+{\epsilon }_T\\ {\epsilon }_1\\ {\epsilon }_2\\ \mathrm{...}\\ {\epsilon }_{T-1}\end{array}\right)}}}.\end{array}$$

21 For example, Storesletten, Telmer, and Yaron (2004) estimated an AR(1) process for the persistent component, whose baseline value for the autoregressive coefficient is 0.96. While they estimated the model in levels, our motivation for estimating (11) in the first-differences is that differences are robust to heterogeneity between cohorts.

22 We compute the Newey-West formula with one lag. Using two or three lags instead has little impact. In this calculation we do not account for the fact that the quantiles are estimated, our rationale being that the cross-sectional sizes are large relative to the length of the time series.

23 Indeed, in the univariate case $X_{i+1,k}\ge X_{ik}$ for all i is equivalent to ${\sum }_{j=1}^m{X}_{{\sigma }_k(i_j),k}({\sigma }_k(i_{j+1})-{\sigma }_k(i_j))\le 0$ for all $m\le N$ and length-m cycle $i_1,\dots ,i_m,i_{m+1}=i_1$.

24 Note that (14) is linear in $X_{\mathcal{l},i,k}$’s. However, it may be impractical to enforce all restrictions in (14) in the update step. In applications, a possibility is to select S_N restrictions at random, where S_N depends on the sample size.

Additional information

Funding

Arellano acknowledges research funding from the Ministerio de Economía y Competitividad, Grant ECO2016-79848-P. Bonhomme acknowledges support from the NSF, Grant SES-1658920.