Recovering Latent Variables by Matching

Recovering Latent Variables by Matching

\baselineskip21pt

\setcounter{page}{0}\thispagestyle{empty}

Introduction

In this paper we propose a method to nonparametrically estimate a class of models with latent variables. We focus on linear factor models whose latent factors are mutually independent. These models have a wide array of economic applications, including measurement error models, fixed-effects models, and error components models. We briefly review the existing literature in Section (ref). In many empirical settings, such as in our application to the study of the cyclical behavior of income shocks, it is appealing not to restrict the functional form of the distributions of latent variables and adopt a nonparametric approach.

Nonparametric estimation based on empirical characteristic functions has been extensively studied in the literature (e.g., Carroll and Hall, 1988; Stefanski and Carroll, 1990). However, while such Fourier-based methods apply to general mutivariate linear factor models with independent components (Horowitz and Markatou, 1996; Li and Vuong, 1998; Bonhomme and Robin, 2010), they tend to be sensitive to the choice of regularization parameters, and they do not guarantee that the estimated densities be non-negative and integrate to one. Recently, Efron (2016) motivated his “parametric g-modeling” approach by the difficulties of nonparametric estimation in this context; see also Efron and Hastie (2016, Chapter 21) and Koenker and Gu (2019).

In this paper we propose a novel nonparametric estimator, and provide evidence that it performs well even in relatively small samples. Our approach differs from the literature in two main aspects. First, we generate a sample of pseudo-observations that may be interpreted as the order statistics of the latent variables. Moments, densities, or other functionals can then be estimated based on them. In particular, densities will be non-negative and integrate to one by construction. Means or other features of the distribution of the latent variables conditional on the data, such as optimal predictors, can also be directly estimated.

The second main feature of our approach is that it is based on matching. Specifically, we generate 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. The model predictions are computed as independent combinations of the pseudo latent observations. This “observation matching” estimation approach can be interpreted as a nonparametric counterpart to (simulated) method-of-moments estimators, which are commonly used in parametric econometric models. Our nonparametric approach, which amounts to minimizing a quadratic Wasserstein distance between empirical distribution functions, exploits linearity and independence to provide a computationally convenient estimator.

As an illustration, in Figure (ref) we show the results of several iterations of our algorithm, in a fixed-effects model with two observation periods and 100 individuals. We start the algorithm from parameter values that are far from the true ones (in the left column). As shown on the top panel, the outcome observations in the data (in crosses) are first matched to model-based predictions (in circles). Pseudo-observations of the latent variables are then updated based on the matched outcome values. The objective function we aim to minimize is the sum of squares of the segments shown on the top panel. The bottom panel shows the estimates of the latent individual-specific effect sorted in ascending order (on the y-axis), against the true values (on the x-axis). We see that, within a few iterations, the model's predictions and the empirical observations tend to agree with each other (in the top panel), and that the distribution of the pseudo latent observations gets close to the population distribution (in the bottom panel).

Our approach builds on and generalizes an important idea due to Colin Mallows (2007), who proposed a “deconvolution by simulation” method based on iterating between sorts of the data and random permutations of pseudo-observations of a latent variable. Mallows (2007) focused on the classical deconvolution model with scalar outcome and known error distribution. Our main goal in this paper is to extend Mallows' insight by proposing a framework to analyze estimators based on matching predicted values from the model to data observations.

In particular, as an extension of Mallows' (2007) original idea, we show how our method can handle multivariate outcomes, hence extending the scope of application to fixed-effects models and multi-factor models. While a number of estimation methods are available for scalar nonparametric deconvolution with known error distribution, the multivariate case -- which is of interest in many economic applications -- remains challenging. Our estimator exploits that the multi-factor models we consider are linear in the independent latent variables, even though they imply nonlinear restrictions on density functions.

A key step in our analysis is to relate the estimation problem to optimal transport theory. Optimal transport is the subject of active research in mathematics, see for example Villani (2003, 2008). Economic applications of optimal transport are many fold, as documented in Galichon (2016). In our context, optimal transport provides a natural way to estimate models with multivariate outcomes via “generalized sorting” algorithms (i.e., matching algorithms) based on linear programming.

To establish the consistency of our estimator we use that, in large samples, our estimator minimizes the Wasserstein distance between the population distribution of the data and the one implied by the model. This problem has a unique solution under suitable conditions on the characteristic functions of the factors (Sz\'ekely and Rao, 2000). Consistency then follows from verifying the conditions for the consistency of sieve extremum estimators (e.g., Chen, 2007) in this setting. When analyzing the multivariate case, our arguments rely on properties of Wasserstein distances established in the optimal transport literature.

We illustrate the performance of our estimator on simulated data. Under various specifications of a nonparametric fixed-effects model, we find that our estimator recovers accurately the true underlying quantile functions and densities, even for samples with only 100 individual observations. This finite-sample performance is remarkable in a fully nonparametric model with multiple latent variables. In addition, we find that our estimator outperforms characteristic-function based estimators, particularly due to improved estimation of the tails of the distributions. In contrast with Fourier methods, our estimator imposes that quantile functions be monotone, and that densities be non-negative. In the related problem of nonparametric instrumental variables estimation, Chetverikov and Wilhelm (2017) show that imposing monotonicity in estimation can help alleviating ill-posedness issues. We conjecture that this feature contributes to explain the finite-sample performance of our estimator.

We then apply our method to study the cyclicality of permanent and transitory income shocks in the US. Answering this question is important, since a well-calibrated cyclical income process is a key input to many models of business cycle dynamics. Storesletten et al. (2004) estimate on the Panel Study of Income Dynamics (PSID) that the dispersion of persistent shocks is countercyclical. However, using a nonparametric descriptive analysis, Guvenen et al. (2014) find using administrative data that the dispersion of log-income growth is acyclical, whereas skewness is procyclical, and Busch et al. (2018) find similar results using the PSID.

We revisit this debate by working with a permanent-transitory model of log-income dynamics, and estimating the annual densities of permanent and transitory shocks nonparametrically. Using the PSID, we estimate that income shocks are not normally distributed, confirming previous evidence using other nonparametric methods. Our main finding is that the dispersion of income shocks is approximately acyclical, whereas the skewness of permanent shocks is procyclical. By comparison, our nonparametric estimates suggest that the dispersion and skewness of shocks to hourly wages vary little with the business cycle.

Our matching-based, minimum Wasserstein distance estimator is related to recent work in machine learning and statistics on the estimation of parametric generative models (see Bernton et al., 2017; Genevay et al., 2017; Bousquet et al., 2017). In contrast with this emerging literature, the models we consider here are nonparametric. In an early theoretical contribution, Bassetti et al. (2006) study consistency in minimum Wasserstein distance estimation. Recently, Rigollet and Weed (2019) develop a minimum Wasserstein deconvolution approach for uncoupled isotonic regression, and Rigollet and Weed (2018) relate maximum-likelihood scalar deconvolution under Gaussian noise to entropic regularized optimal transport. Lastly, our general estimation strategy is also related to Galichon and Henry's (2011) analysis of partially identified models.

As we show at the end of the paper, our matching approach can be generalized to nonparametric estimation of other latent variables models. We briefly describe such generalizations to cross-sectional random coefficients models with exogenous covariates (Beran and Hall, 1992; Ichimura and Thompson, 1998), panel data random coefficients models (Arellano and Bonhomme, 2012), nonparametric deconvolution under heteroskedasticity (Delaigle and Meister, 2008), and nonparametric finite mixture models (Hall and Zhou, 2003).

The outline of the paper is as follows. In Section (ref) we describe linear independent factor models, and we briefly review applications and existing estimation approaches. In Section (ref) we introduce our matching estimator. In Sections (ref) and (ref) we study computation and consistency, respectively. In Sections (ref) and (ref) we present the simulation exercises and empirical application. In Section (ref) we outline several extensions. Lastly, we conclude in Section (ref). Proofs and additional material are collected in the appendix.

Independent factor models

We focus on linear independent factor models of the form $Y=AX$, where $Y=(Y_1,...,Y_T)'$, $X=(X_{1},...,X_{K})'$, $A$ is a known or consistently estimable $T\times K$ matrix, and the components $X_1,...,X_K$ are mutually independent. In this section we review several examples of models and applications that have such a structure. We focus on the case $K>T$, so the system is singular and the realizations of the latent variables are not identifiable, although under suitable conditions their distributions will be.

\paragraph{Nonparametric deconvolution.}

When $T=1$, $Y=X_1+X_2$, and $X_2$ has a known or consistently estimable distribution, one obtains the scalar nonparametric deconvolution model. This model has been extensively studied in statistics and econometrics. Nonparametric deconvolution is often used to deal with the presence of measurement error. In such settings, $Y$ is an error-ridden variable, $X_1$ is the true value of the variable, and $X_2$ is an independent, classical measurement error (e.g., Carroll et al., 2006; Chen et al., 2011; Schennach, 2013a). Other economic applications of nonparametric deconvolution are 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 et al., 1997; Wu and Perloff, 2006), and the estimation of the distribution of time-invariant random coefficients of binary treatments in panel data models (Arellano and Bonhomme, 2012).

The statistical literature on nonparametric deconvolution provides conditions under which the distribution of $X_1$ is nonparametrically identified, alongside numerous estimation approaches such as kernel deconvolution estimators (Carroll and Hall, 1988; Delaigle and Gijbels, 2002; Fan, 1991), wavelet methods (Pensky and Vidakovic, 1999; Fan and Koo, 2002), regularization techniques (Carrasco and Florens, 2011), and nonparametric maximum likelihood methods (Kiefer and Wolfowitz, 1956; Gu and Koenker, 2017).

\paragraph{Nonparametric distribution of fixed effects.}

A leading example of a linear independent factor model is the fixed-effects model:

where $Y_1,...,Y_T$ are observed outcomes and $\alpha,\varepsilon_1,...,\varepsilon_T$ are latent and mutually independent. Working with $T=2$, Kotlarski (1967) provided simple conditions under which the density functions of the latent factors are nonparametrically identified in model ((ref)).

This fixed-effects structure arises frequently in economic applications. As an example, $\alpha$ can be a latent skill of an individual, measured with error (as in Cunha et al., 2010). In other applications, researchers may be interested in estimating the distribution of worker, teacher, firm, school, hospital, neighborhood, or bank-specific fixed-effects, for example. Compared to commonly used Gaussian specifications (e.g., Kane and Staiger, 2008; Angrist et al., 2017; Chetty and Hendren, 2018), a nonparametric estimator of the distribution of $\alpha$ in ((ref)) will be robust to functional form assumptions under the maintained assumption of mutual independence. Non-Gaussianity, such as skewness or fat tail behavior, is relevant in many empirical settings. The fixed-effects model ((ref)) and its generalizations are sometimes estimated using flexible parametric specifications such as finite Gaussian mixtures (e.g., Carneiro et al., 2003). Alternatively, nonparametric estimators based on empirical characteristic functions can be constructed, by mimicking and extending Kotlarski's proof (Li and Vuong, 1998; Li, 2002; Horowitz and Markatou; 1996).

\paragraph{Error components: generalized nonparametric deconvolution.}

A prominent error component model is the permanent-transitory model for the dynamics of log-income: $Y_t=\eta_t+\varepsilon_t$, where $\eta_t=\eta_{t-1}+v_t$ is a random walk with independent innovations, and all $\varepsilon_t$'s and $v_t$'s are independent over time and independent of each other and of the initial $\eta_0$ (e.g., Hall and Mishkin, 1982; Blundell et al., 2008). This model is a special case of a linear independent factor model $Y=AX$, where $Y=(Y_1,...,Y_T)'$ are observed outcomes, $X=(X_1,...,X_K)'$ are mutually independent latent factors, and $A$ is a known $T\times K$ matrix. Identification of such generalized deconvolution models is established in Sz\'ekely and Rao (2000). Bonhomme and Robin (2010) propose nonparametric characteristic-function based estimators of factor densities, and apply them to study income dynamics; see also Botosaru and Sasaki (2015).\footnote{Quantile-based estimation in linear and nonlinear factor models is introduced in Arellano and Bonhomme (2016), and applied in Arellano et al. (2017) to document income dynamics in the PSID.} In such settings, a nonparametric approach is able to capture the skewness and kurtosis of income shocks. In addition, an important application of error components models is to relax independence in fixed-effects models such as ((ref)). This can be done provided $T$ is large enough.\footnote{Modeling $\varepsilon_t$ in ((ref)) 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 et al., 2019). Ben Moshe (2017) shows how to allow for arbitrary subsets of dependent factors, and proposes characteristic-function based estimators. In addition, in model ((ref)) Schennach (2013b) points out that full independence between the factors is not necessary, and that sub-independence suffices to establish identification.} Such specifications can be estimated using the methods we introduce in this paper.

Latent variable estimation by matching

In this section, to introduce the main ideas we start by describing our estimator in the scalar nonparametric deconvolution model. We then show how the same approach can be used to estimate linear multi-factor models with independent factors.

Nonparametric deconvolution

Let $Y=X_1+X_2$ be a scalar outcome, where $X_1$ and $X_2$ are independent, $X_1$ is unobserved to the econometrician, and its distribution is unspecified. We assume that $Y$, $X_1$ and $X_2$ are continuously distributed, and postpone more specific assumptions until Section (ref). Let $F_Z$ denote the cumulative distribution function (cdf) of any random variable $Z$. We assume that two random samples, $Y_1,...,Y_N$ and $X_{12},...,X_{N2}$, drawn from $F_Y$ and $F_{X_2}$, respectively, are available.\footnote{The sample sizes being the same for $Y$ and $X_2$ 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.}

Our goal is to estimate a sample of pseudo-observations $\widehat{X}_{11},...,\widehat{X}_{N1}$, whose empirical cdf is asymptotically distributed as $F_{X_1}$ as $N$ tends to infinity. To do so, we minimize a distance between the sample of observed $Y$'s and a sample of $Y$'s predicted by the model. We rely on the quadratic Wasserstein distance (see, e.g., Chapter 7 in Villani, 2003), which is the minimum Euclidean distance between observed $Y$'s and predicted $Y$'s with respect to all possible reorderings of the observations.

Formally, assume without loss of generality that $Y_i\leq Y_{i+1}$ and $X_{i2}\leq X_{i+1,2}$ for all $i$. Let $\Pi_N$ denote the set of permutations $\pi$ of $\{1,...,N\}$. Moreover, let $\overline{C}_N>0$ and $\underline{C}_N>0$ be two constants, and let ${\cal{X}}_{N}$ be the set of parameter vectors $X_1=(X_{11},...,X_{N1})\in\mathbb{R}^N$ such that $|X_{i1}|\leq \overline{C}_N$ and $\underline{C}_N\leq (N+1)(X_{i+1,1}-X_{i1})\leq \overline{C}_N$ for all $i$. The constants $\underline{C}_N$ and $\overline{C}_N$ play a role in our consistency argument below, and we will study how their choice affects our estimator in simulations. We propose to compute:

where $\sigma$ is a random permutation in $\Pi_N$ (i.e., a uniform draw on $\Pi_N$), independent of $Y_1,...,Y_N,X_{12},...,X_{N2}$.

To interpret the objective function on the right-hand side of ((ref)), note that, for any random permutation $\sigma$, $Z_i\equiv X_{\sigma(i),1}+X_{i,2}$, $i=1,...,N$, are $N$ draws from the model. Predicted values from the model could be generated in other ways. For example, one could instead compute $X_{i1}+\widetilde{X}_{i2}$, where $\widetilde{X}_{i2}$ are i.i.d. draws from the empirical distribution of $X_{i2}$. Alternatively, one could generate $R>1$ predictions per observation $i$, although here we take $R=1$ to minimize computation cost.\footnote{Specifically, one could compute $X_{\sigma(i,r),1}+X_{i2}$, with $\sigma(\cdot,1),...\sigma(\cdot,R)$ being $R$ independent permutations. In that case, $\pi$ would be a generalized permutation (or “pure matching”), mapping $\{1,...,N\}^R$ to $\{1,...,N\}$.}

A simple way to reduce the dependence of the estimator on the random $\sigma$ draw is to compute $\widehat{X}^{(m)}_{i1}$, for $i=1,...,N$ and $m=1,...,M$, where $\sigma^{(1)},...,\sigma^{(M)}$ are independent random permutations drawn from $\Pi_N$, and to report the averages: $\widehat{X}_{i1}=\frac{1}{M}\sum_{m=1}^M\widehat{X}^{(m)}_{i1}$, for $i=1,...,N$. For fixed $M$, such averages will be consistent as $N$ tends to infinity under similar conditions as our baseline estimator.

The estimator $\widehat{X}_1$ in ((ref)) minimizes the Wasserstein distance between the empirical distributions of the model predictions $Z_i=X_{\sigma(i),1}+X_{i2}$ and the outcome observations $Y_i$. The Wasserstein distance is defined as:

Since $Y_i$ and $Z_i$ are scalar, the Hardy-Littlewood-P\'olya rearrangement inequality implies that the solution to ((ref)) is to sort $Y_i$'s and $Z_i$'s in the same order. That is, letting $\widehat{\pi}$ denote the minimum argument in ((ref)), $\widehat{\pi}(i)=\limfunc{Rank}(Z_i)\equiv N\widehat{F}_Z(Z_i)$ is the rank of $Z_i$.

Nonparametric factor models

We now apply the same idea to a general linear independent multi-factor model $Y=AX$, where $A$ is a $T\times K$ matrix with generic element $a_{tk}$, and $X=(X_1,...,X_K)'$ with $X_1,...,X_K$ mutually independent. For simplicity we assume that $X$ and $Y$ have zero mean.\footnote{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 fixed-effects model, assuming that $\mathbb{E}(X_1)=0$ suffices for identification. Our algorithm can easily be adapted to such cases.} We seek to compute pseudo-observations $\widehat{X}_{11},...,\widehat{X}_{N1}$, ..., $\widehat{X}_{1K},...,\widehat{X}_{NK}$, which minimize the Wasserstein distance between the sample of observed $Y$'s, which here are $T\times 1$ vectors, and the sample of $Y$'s predicted by the factor model.

As before, let $\overline{C}_N>0$ and $\underline{C}_N>0$ be two constants, and let ${\cal{X}}_N$ be the set of $(X_1,...,X_N)\in\mathbb{R}^{NK}$ such that $|X_{i,k}|\leq \overline{C}_N$ and $\underline{C}_N\leq (N+1)(X_{i+1,k}-X_{ik})\leq \overline{C}_N$ for all $i$ and $k$, and $\sum_{i=1}^N X_{ik}=0$ for all $k$. We define:

where $\sigma_1,...,\sigma_K$ are independent random permutations in $\Pi_N$, independent of $Y_{11},...,Y_{NT}$.

As in the scalar case, $Z_{it}\equiv\sum_{k=1}^Ka_{tk}X_{\sigma_k(i),k}$, $i=1,...,N$, $t=1,...,T$, are $NT$ predicted values from the factor model. Hence, as before, the vector $\widehat{X}$ minimizes the Wasserstein distance between the empirical distributions of the data $(Y_{i1},...,Y_{iT})$ and of the model predictions $(Z_{i1},...,Z_{iT})$. A difference with the scalar deconvolution model is that, when $Y_i$ are multivariate, the minimization with respect to $\pi$ inside the brackets in ((ref)) does not have an explicit form in general. However, from optimal transport theory it is well-known that the solution can be obtained as the solution to a linear program. We will exploit this feature in our estimation algorithm.

\paragraph{Densities and expectations.}

In Section (ref) we will provide conditions under which $\widehat{X}_{ik}$, $i=1,...,N$, consistently estimate the quantile function of $X_k$. More precisely, we will show that $\max_{i=1,...,N}\, |\widehat{X}_{ik}-F_{X_k}^{-1}(\frac{i}{N+1})|$ tends to zero in probability asymptotically. This provides uniformly consistent estimators of the quantile functions of the latent variables, which can in turn be used for density estimation under a slight modification of the parameter space ${\cal{X}}_N$. Indeed, let us restrict the parameter space to elements $X=(X_1,...,X_N)$ in ${\cal{X}}_N$ which satisfy the following additional restrictions on second-order differences: $(N+1)^2\left|X_{i+2,k}-2X_{i+1,k}+X_{i,k}\right|\leq \overline{C}_N$, for all $i$ and $k$. Let us then define, for a bandwidth parameter $b>0$ and a kernel function $\kappa\geq 0$ that integrates to one:

We will show that $\widehat{f}_{X_k}$ is uniformly consistent for the density of $X_k$, under standard conditions on the kernel $\kappa$ and bandwidth $b$.

In addition, our estimator delivers simple consistent estimators of unconditional and conditional expectations, as we show in Appendix (ref). As an example of practical interest, in the fixed-effects model ((ref)) the best predictor of $X_1$ under squared loss can be estimated as:

where the weights $\widehat{\omega}_{i}$ are given by: $$\widehat{\omega}_{i}=\frac{\prod_{t=1}^T\widehat{f}_{X_{t+1}}(Y_{it}-\widehat{X}_{i1})}{\sum_{j=1}^N\prod_{t=1}^T\widehat{f}_{X_{t+1}}(Y_{jt}-\widehat{X}_{j1})},\quad i=1,...,N.$$

Computation

The optimization problems in ((ref)) and ((ref)) are mixed integer quadratic programs. Although the literature on mixed integer programming has recently made substantial progress (e.g., Bliek et al., 2014), exact algorithms are currently limited in the dimensions they can allow for. Here we describe a simple, practical method to minimize ((ref)) and ((ref)).

Algorithm

The algorithm we propose is based on the observation that, for given $X_{11},...,X_{NK}$ values, ((ref)) is a linear assignment (or discrete optimal transport) problem, hence it can be solved by any linear programming routine. In turn, given $\pi$, ((ref)) is a monotone least squares problem. Our estimation algorithm is as follows. Here we focus on the general form ((ref)), since the estimator for the scalar deconvolution model ((ref)) is a special case of it.

Both steps in the algorithm are straightforward to implement. The matching step ((ref)) can be computed by a linear programming routine, due to the fact that the linear programming relaxation of a discrete optimal transport problem has integer-valued solutions.\footnote{See for example Chapter 3 in Galichon (2016) on discrete Monge-Kantorovitch problems, and Conforti et al. (2014) on integer programming problems and perfect formulations.} Formally, $\widehat{\pi}^{(s+1)}$ in ((ref)) is a solution to the following linear program:

where ${\cal{P}}_{N}$ denotes the set of $N\times N$ matrices with non-negative elements, whose rows and columns all sum to one. In the scalar nonparametric deconvolution case ((ref)), this gives $\widehat{\pi}^{(s+1)}(i)=\widehat{\limfunc{Rank}}\left(\widehat{X}^{(s)}_{\sigma(i),1}+X_{i2}\right)$ for all $i$.

In fact, it is possible to write $\widehat{X}=(\widehat{X}_{1},...,\widehat{X}_{N})$ in ((ref)) as the solution to a quadratic program:

which is not convex in general. Our estimation algorithm is a method to solve this non-convex quadratic program. However, the algorithm is not guaranteed to reach a global minimum in ((ref)). Our implementation is based on starting the algorithm from multiple random values. We will assess the impact of starting values on simulated data.

Comparison to Mallows (2007)

Our algorithm may be seen as a generalization of Mallows' (2007) “deconvolution by simulation” method. To highlight the connection, consider the scalar nonparametric deconvolution model. The two steps in our algorithm take the following form:

The Mallows (2007) algorithm is closely related to this algorithm. The main difference is that, instead of minimizing an objective function for fixed values of the random permutation $\sigma$, random permutations are re-drawn in each step of the algorithm. In addition, the ordering of the $X_{i1}$'s is not restricted, and neither are the values and increments of the $X_{i1}$'s. Formally, the sub-steps of the Mallows algorithm are the following:

To provide intuition about this algorithm, Mallows (2007) observes that, starting with draws from the true latent $X_1$, one expects the iteration to continue to draw from that distribution. However, starting from different values, the $\widehat{X}_1$ vectors implied by the algorithm will follow a complex, $N$-dimensional Markov Chain. Moreover, the consistency properties of the Mallows estimator are currently unknown. Lastly, note that the methods introduced in this paper naturally deliver counterparts to the Mallows algorithm for other models beyond deconvolution, such as general linear independent factor models.

Consistency analysis

In this section we provide conditions under which the estimators introduced in Section (ref) are consistent.

For $k\in\{1,...,K\}$, let us denote the quantile function of $X_k$ as: $$F_{X_k}^{-1}(\tau)=\limfunc{inf}\, \{x\in\limfunc{Supp}(X_k)\, :\, F_{X_k}(x)\geq \tau\},\text{ for all }\tau\in(0,1).$$ In addition, for any candidate quantile function $H_k$ that maps the unit interval to the real line, let us define the following Sobolev sup-norms: $$\|H_k\|_{\infty} =\sup_{\tau\in(0,1)}\, |H_k(\tau)|, \,\,\,\text{ and }\,\,\, \|H_k\| =\max_{m\in\{0,1\}}\sup_{\tau\in(0,1)}\, |\nabla^{m}H_k(\tau)|,$$ where $\nabla^{m}H_k$ denotes the $m$-th derivative of $H_k$ (when it exists). We will simply denote $\nabla=\nabla^1$ for the first derivative.

To a solution $\widehat{X}_k$ to ((ref)),\footnote{It is not necessary for $\widehat{X}_k$ to be an exact minimizer of ((ref)). As we show in the proof, it suffices that the value of the objective function at $(\widehat{X}_1,...,\widehat{X}_K)$ be in an $\epsilon_N$-neighborhood of the global minimum, for $\epsilon_N$ tending to zero as $N$ tends to infinity.} we will associate an interpolating quantile function $\widehat{H}_k$ such that $\widehat{H}_k\left(\frac{i}{N+1}\right)=\widehat{X}_{ik}$ for all $i$. We will then show that $ \|\widehat{H}_k-F_{X_k}^{-1}\|_{\infty}=o_p(1)$. This result will be obtained as an application of the consistency theorem for sieve extremum estimators in Chen (2007).

We make the following assumptions.

Though convenient for the derivations, the compact supports assumption in part $(i)$ is strong. This could be relaxed by working with weighted norms, at the cost of achieving a weaker consistency result. The simulation experiments we report below suggest that the estimator continues to perform well when supports are unbounded. Part $(ii)$ is a sufficient condition for the distributions of latent variables $X_k$ to be nonparametrically identified (Sz\'ekely and Rao, 2000). The constants $\underline{C}_N$ and $\overline{C}_N$ appearing in part $(iii)$ ensure that the $\widehat{X}_{ik}$ values are bounded and of bounded variation. Lastly, the independent random permutations $\sigma_1,...,\sigma_K$ in ((ref)) depend on $N$, although we have omitted this dependence for conciseness.

Consistency is established in the following theorem. Proofs are in Appendix (ref).

While Theorem (ref) does not formally cover the scalar deconvolution model, the same proof arguments can be used to show the following result, under similar assumptions to those of Theorem (ref).

An important step in the proof of Theorem (ref) is to define the population counterpart to the estimation problem ((ref)). Let $\mu_Y$ denote the population measure of $Y$. Moreover, for any candidate quantile functions $H=(H_1,...,H_K)$, let $\mu_{AH}$ denote the population measure of the random vector $Z\equiv \sum_{k=1}^KA_kH_k(V_k)$, where $V_1,...,V_K$ are independent standard uniform random variables on the unit interval. Finally, let ${\cal{M}}(\mu_Y,\,\mu_{AH})$ denote the set all possible joint distributions, or couplings, of the random vectors $Y$ and $\sum_{k=1}^KA_k H_k\left(V_k\right)$, with marginals $\mu_Y$ and $\mu_{AH}$. The population objective function is then: $$Q(H)\equiv\underset{\pi\in {\cal{M}}(\mu_Y,\,\mu_{AH})}{\limfunc{inf}}\,\mathbb{E}_{\pi}\left[\sum_{t=1}^T\left (Y_t-\sum_{k=1}^Ka_{tk} H_k\left(V_k\right)\right)^2\right],$$ which is the quadratic Wasserstein distance between the population distribution of the data and the one implied by the model. Under part $(ii)$ in Assumption (ref) that ensures identification, $Q(H)$ is minimized at the true quantile functions $H_k=F_{X_k}^{-1}$.

In the scalar deconvolution model, the population objective takes the explicit form: $$Q(H_1)\equiv\mathbb{E}\left[\left (F_Y^{-1}\left(\int_0^1 F_{X_2}\left(H_1(V_1)+F_{X_2}^{-1}(V_2)-H_1(\tau)\right)d\tau\right)-H_1(V_1)-F_{X_2}^{-1}(V_2)\right)^2\right],$$ where the expectation is taken with respect to independent standard uniform random variables $V_1$ and $V_2$. Note that the integral in this expression is simply the population rank of $H_1(V_1)+F_{X_2}^{-1}(V_2)$. When the characteristic function of $X_2$ has no real zeros, $Q(H_1)$ is minimized at $H_1=F_{X_1}^{-1}$.

\paragraph{Densities and expectations.}

Under slightly stronger assumptions, Theorem (ref) can be modified to obtain consistent estimators of both $F_{X_k}^{-1}$ and its derivative, which can then be used for density estimation. To see this, let us denote as ${\cal{X}}^{(2)}_N$ the set of $X$ in ${\cal{X}}_N$ which satisfy the restrictions on second-order differences: $(N+1)^2\left|X_{i+2,k}-2X_{i+1,k}+X_{ik}\right|\leq \overline{C}_N$, for all $i$ and $k$, and replace the minimization in ((ref)) by a minimization with respect to $X\in {\cal{X}}^{(2)}_N$. Imposing in Assumption (ref) that the densities of $X_k$ have bounded second-order derivatives, and modifying the proof of Theorem (ref) accordingly, we obtain that:

We then have the following result.

Lastly, given Corollary (ref) it can readily be checked that conditional expectations estimators, such as ((ref)) and those in Appendix (ref), are consistent in sup-norm for their population counterparts.

\paragraph{Remark: convergence rates and inference.}

It follows from existing convergence rates in nonparametric deconvolution models (e.g., Fan, 1991; Hall and Lahiri, 2008) that neither $\widehat{X}_{ik}$ (as an estimator of the quantile function of $X_k$) nor its functionals will converge at the root-$N$ rate in general. Bertail et al. (1999) propose an inference method under the condition that the estimator is $N^{\beta}$-consistent with a continuous asymptotic distribution, for some $\beta>0$. Their rate-adaptive method is attractive in our setting, although polynomial convergence rates may rule out cases of severe ill-posedness. Completing the characterization of the asymptotic behavior of our estimator is an important task for future work.

Performance on simulated data

In this section we illustrate the finite-sample performance of our estimator on data simulated from a nonparametric fixed-effects model. In Appendix (ref) we report additional results for a scalar nonparametric deconvolution model.

We focus on the model $Y_1=X_1+X_2$, $Y_2=X_1+X_3$, where $X_1,X_2,X_3$ are independent of each other and have identical distributions. We consider four specifications for the distribution of $X_k$ for all $k$: Beta$(2,2)$, Beta$(5,2)$, normal, and log-normal, all standardized so that $X_k$ has mean zero and variance one. To restrict the maximum values of $\widehat{X}_{ik}$, its increments, and its second-order differences, we consider two choices for the penalization constants: $(\underline{C}_N,\overline{C}_N)=(.1,10)$ (“strong constraint”), and $(\underline{C}_N,\overline{C}_N)=(0,10000)$ (“weak constraint”). To minimize the objective function in ((ref)) we start with $10$ randomly generated starting values, drawn from widely dispersed mixtures of five Gaussian distributions, and keep the solution corresponding to the minimum value of the objective. Lastly, we draw $M=10$ independent random permutations in $\Pi_N$, and average the resulting $M$ sets of estimates $\widehat{X}^{(m)}_{i1}$, for $i=1,...,N$.

In Appendix (ref) we study the sensitivity of the estimates to the penalization constants, the starting values, and the number $M$ of $\sigma$ draws, in a nonparametric deconvolution model. We find that the estimator is quite robust to these choices. In particular, we document that taking conservative choices for $\underline{C}_N$ and $\overline{C}_N$ (such as in the “weak constraint” case) results in a well-behaved estimator, suggesting that our matching procedure induces an implicit regularization, even in the absence of additional constraints on parameters. At the same time, we find that such a conservative choice may not be optimal in terms of mean squared errors of quantile estimates. The optimal choice of penalization constants is an interesting question for future work.\footnote{A simple recommendation for practice could be based on a truncated normal distribution. Let $\widehat{\sigma}_k$ denote a consistent estimate of the standard deviation of $X_k$, e.g. 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 $3275c\widehat{\sigma}_k$ (upper bound on second derivatives). When $c=1$, these constants are binding when $X_k$ follows a normal truncated at the 99th percentiles. A default choice could be $c=2$.}

In the first two columns in Figure (ref) we show the estimates of the quantile functions $\widehat{X}_{i1}=\widehat{F}_{X_1}^{-1}\left(\frac{i}{N+1}\right)$, for the four specifications and both penalization parameters. The results for the other two factors are similar and omitted for brevity. The solid and dashed lines correspond to the mean and 10 and 90 percentiles across 100 simulations, respectively, while the dashed-dotted line corresponds to the true quantile function. The sample size is $N=100$. Even for such a small sample size, our nonparametric estimator performs well, especially under a weaker constraint on the parameters (second column). In the last two columns of Figure (ref) we show density estimates for the same specifications. We take a Gaussian kernel and set the bandwidth based on Silverman's rule. Although there are some biases in the strong constraint case, our nonparametric estimator reproduces the shape of the unknown densities well.

In Table (ref) we report the mean integrated squared and absolute errors (MISE and MIAE, respectively) of our density estimators, for the four distributional specifications and $N=100$. We see that the estimator performs better under the weak constraint. Moreover, interestingly, as shown by the last two columns of Table (ref) our estimator outperforms characteristic-function based density estimators. Here the “Fourier” results are based on the estimator of Bonhomme and Robin (2010), and we use their recommended choice to set the regularization parameter in each replication. Inspection of the estimates suggests that the differences are mainly driven by estimates of the tails of the densities. Characteristic-function based estimators do not guarantee that densities be non-negative, and we find that their values tend to oscillate in the left and right tails. From results in Chetverikov and Wilhelm (2017), we conjecture that finite-sample performance benefits from the fact that our estimator enforces monotonicity of quantile functions and non-negativity of densities. However, proving this conjecture would require providing convergence rate results in addition to our consistency analysis.

Lastly, in Appendix (ref) we present numerical calculations of the rate of convergence of our estimator of latent quantiles, in data simulated from a nonparametric scalar deconvolution model. The results suggest the rate ranges between $N^{-\frac{3}{10}}$ and $N^{-\frac{7}{10}}$ in the data generating processes that we study. We also compare the performance of our method to Mallows' (2007) “deconvolution by simulation” estimator.

Empirical application: income risk over the business cycle in the PSID

In this section we use our nonparametric method to study the cyclical behavior of income risk in the US. In an influential contribution, Storesletten et al. (2004) report using the PSID that the dispersion of idiosyncratic income shocks increases substantially in recessions. Guvenen et al. (2014) re-examine this finding, using US administrative data and focusing on log-income growth. They find that the dispersion of log-income growth is acyclical, and that its skewness is procyclical. Recently, Busch et al. (2018) find similar results using the PSID and data from Sweden and Germany. Nakajima and Smyrnyagin (2019) use an approach similar to the one in Storesletten et al. (2004), making use of a larger PSID sample and different measures of income, and find that log-income shocks exhibit countercyclical dispersion and procyclical skewness. This literature is motivated by the key quantitative role of the cyclical behavior of the income process when calibrating models of business cycle dynamics.

Here we revisit this question, by estimating a nonparametric permanent-transitory model where log-income, net of the effect of some covariates, is the sum of a random walk $\eta_{it}=\eta_{i,t-1}+v_{it}$ and an independent innovation $\varepsilon_{it}$. In first-differences we have, denoting log-income growth as $\Delta Y_{it}=Y_{it}-Y_{i,t-1}$:

Model ((ref)) is a linear factor model with $2T-1$ independent factors. Indeed, we have: $$\underset{\equiv Y}{\underbrace{\left(

\right)}}=\underset{\equiv A}{\underbrace{\left(

\right)}}\,\underset{\equiv X}{\underbrace{\left(

\right)}}.$$ We leave the distributions of $v_{it}$ and $\varepsilon_{it}$ unrestricted. Our aim is to document the behavior of these distributions over the business cycle.

There are several differences between our model and estimation approach and the ones in Storesletten et al. (2004). A substantive difference is that we estimate the densities of the shocks nonparametrically, while they use a parametric model under Gaussian assumptions. This is important, since estimates of non-Gaussian models (e.g., Horowitz and Markatou, 1996; Geweke and Keane, 2000; Bonhomme and Robin, 2010; Arellano et al., 2017) and descriptive evidence (e.g., Guvenen et al., 2014; Guvenen et al., 2016) both suggest that income shocks are strongly non-Gaussian in the US. Another difference is that we rely on first-differences of log-income in estimation, while Storesletten et al. (2004) estimate the model in levels. This choice allows them to exploit a long past history of recessions, even before the PSID started to be collected, since past recessions and expansions affect the cross-sectional variance of log-income. On the other hand, income levels may also reflect other differences between cohorts, and our estimation in differences is robust to those. A last, less substantive difference is that we impose that the persistent component follows a unit root, while Storesletten et al. (2004) use an autoregressive process whose baseline value for the autoregressive coefficient is 0.96.

Studying aggregate dynamics using survey panel data like the PSID is complicated by attrition and confounding age effects. To minimize the impact of these factors, we follow the approach pioneered by Storesletten et al. (2004) and construct a sequence of balanced, four-year subpanels. In every subpanel, we require that households have non-missing data on income and demographics and comply with standard selection criteria: the household has positive annual labor income during the four years, the head is between 23 and 60 years old, and is not part of the SEO low-income sample or the immigrant sample. We estimate model ((ref)) on 21 subpanels, whose base years range between 1969 and 1989. Log-household income growth is net of indicators for age (of head), education, gender, race, marital status, state of residence, number of children, and family size. In estimation we set conservative values for the penalization constants (that is, we use the “weak constraint” values of the simulation section), we use a single starting value in the algorithm, and we average the results of $M=10$ draws.

Our first finding is that income shocks are strongly non-Gaussian. In Figure (ref) we report the estimated quantile functions and densities of permanent shocks $v_{it}$ and transitory shocks $\varepsilon_{it}$, averaged over years (in solid), together with normal fits (in dashed). The excess kurtosis of both shocks is in line with previous evidence reported in the literature (e.g., Geweke and Keane, 2000; Bonhomme and Robin, 2010).

We are interested in how features of these distributions vary with the business cycle. In the left column of Figure (ref) we plot the 90/10 percentile difference of log-income $P_{90}-P_{10}$ (a common measure of dispersion, in solid) together with log-GDP growth (in dashed), both of them net of a linear time trend. While permanent and transitory shocks tend to move countercyclically in the first part of the period, the relationship tends to become procyclical in the 1980's. As we report in Table (ref), the coefficient of log-GDP growth in a regression of the dispersion of permanent income shocks on log-GDP growth and a time trend is -0.25, with a Newey-West standard error of 0.30.\footnote{We compute the Newey-West formula with one lag. Using two or three lags instead has little impact. In the computation 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.} Hence, overall we do not find significant evidence that the dispersion of permanent shocks varies systematically with the business cycle. This result based on first-differenced estimation and a nonparametric approach contrasts with the main finding in Storesletten et al. (2004). In addition, we neither find that the dispersion of transitory shocks varies with the cycle.

Next, in the right column of Figure (ref) we plot the Bowley-Kelley quantile measure of skewness $[(P_{90}-P_{50})-(P_{50}-P_{10})]/(P_{90}-P_{10})$. The graphs of permanent and transitory income shocks suggest that skewness is procyclical. This is confirmed in Table (ref), which shows that the coefficient of log-GDP growth in the skewness regression is 3.07 for permanent shocks, and 2.36 for transitory shocks, significant at the 5% level in both cases. Our nonparametric estimates of a permanent-transitory model of income dynamics thus suggest that dispersion is approximately acyclical, and skewness is procyclical, in line with the conclusions of the descriptive evidence in Guvenen et al. (2014) and Busch et al. (2018).

As graphical way to illustrate the distributional dynamics of income over the business cycle, in Figure (ref) we plot the coefficients of log-GDP growth in regressions of the quantiles of permanent or transitory income shocks on log-GDP growth and a time trend. The estimates suggest a U-shape pattern along the distribution, both for permanent and transitory shocks. Expansions are associated with increases at the top and bottom of the distribution, while recessions are associated with the opposite pattern and a relative increase of the middle quantiles. In the upper panel of Figure (ref) we show how the model fits the distributions of log-income growth, suggesting that our model is able to reproduce the density and quantile cyclicality of log-income growth that we observe in the data.

We performed several exercises to probe the robustness of these findings, using the “strong constraint” penalization of Section (ref), measuring business cycle conditions using the unemployment rate instead of log-GDP growth, and varying the choice of starting values in the algorithm. While we found the year-to-year variation in Figure (ref) to depend on the chosen specification, in all our checks we found a lack of systematic cyclical variability of the dispersion of income shocks, and a significant procyclicality of the skewness of permanent shocks. Among the results reported in Table (ref), we found the procyclicality of the skewness of transitory shocks to be most sensitive to specification changes.

\paragraph{Hourly wages.}

We next use the information in the PSID about hours worked to compute similar measures of cyclicality based on hourly wages of household heads. Evidence from Italy and France (Hoffmann and Malacrino, 2019; Pora and Wilner, 2019) suggests that days and hours worked may contribute significantly to the observed cyclical patterns of skewness. For the US, Nakajima and Smyrnyagin (2019) obtain similar conclusions. In contrast, Busch et al. (2018) find a moderate role of hours worked in Germany. In the bottom panel of Table (ref) we see that the skewnesses of permanent and transitory shocks to hourly wages do not vary significantly with the cycle, and that the point estimates are greatly reduced compared to the case of total income. This suggests that hours worked largely contribute to the distributional income dynamics that we document. In the lower panels in Figure (ref) we show the model fit to log-hourly wage growth. The estimates show that quantiles of log-hourly wage growth vary little with the business cycle in our sample, and that our model is able to reproduce this pattern.

Extensions

In this section we briefly outline several extensions of our matching approach to random coefficients models, finite mixture models, and deconvolution models with heteroskedasticity. These extensions show that the idea of matching data observations to model predictions is applicable to a variety of settings with latent variables.

Random coefficients

Consider the linear cross-sectional random coefficients model:

where $(W_2,...,W_K)$ is independent of $(X_1,...,X_K)$, the scalar outcome $Y$ and the covariates $W_2,...,W_K$ are observed, and $(X_1,...,X_K)$ is a latent vector with an unrestricted joint distribution (e.g., Beran and Hall, 1992; Hoderlein et al., 2010). To construct a matching estimator in this case, we augment ((ref)) with: $W_k=V_k$, $k=2,...,K$, where the $V_k$'s are auxiliary latent variables independent of the $X_k$'s. In this augmented model, the joint distributions of $(X_1,...,X_K)$ and $(V_2,...,V_K)$ can be estimated by minimizing the Euclidean distance between the model's predictions of $Y,W$ observations, and their matched values in the data. A similar approach can be used in binary choice models with random coefficients (Ichimura and Thompson, 1998; Gautier and Kitamura, 2013)

To see how to adapt this idea to panel data random coefficients models, consider the random trends model:

where $(\alpha_i,\beta_i)$, $\varepsilon_{i1}$, ..., $\varepsilon_{iT}$ are mutually independent. Our matching approach applies directly to this case, by minimizing the following objective:

where $\sigma_1,...,\sigma_T$ are independent random permutations of $\{1,...,N\}$. In this case our algorithm consists in alternating optimal transport (matching) steps and least squares (update) steps. Note that in this case the algorithm delivers bivariate pseudo-observations $(\widehat{\alpha}_i,\widehat{\beta}_i)$, and that those can no longer be interpreted as estimates of order statistics; see Chernozhukov et al. (2017) for an optimal transport approach to multivariate quantiles.\footnote{When the trend $t$ in ((ref)) is replaced by a strictly exogenous regressor $X_{it}$, we can augment the model with auxiliary latent variables $V_{it}$ using the same strategy as in the cross-sectional case, and minimize the distance between the model's predictions of $Y,X$ and their matched values in the data. Note that in that case $X_{it}$ and $(\alpha_i,\beta_i)$ are allowed to be dependent.}

Finite mixtures

Consider next a finite mixture model with $G$ groups, for a $T$-dimensional outcome $Y$:

where $Z_1,...,Z_G$ and $X_{11},...,X_{GT}$ are unobserved, $Z_g\in\{0,1\}$ with $\sum_{g=1}^G Z_g=1$, and $(Z_1,...,Z_G)$ and all $X_{11}$, ..., $X_{GT}$ are mutually independent. The nonparametric model ((ref)) has been extensively analyzed in the literature (e.g., Hall and Zhou, 2003; Hu, 2008; Allmann et al., 2009; Bonhomme et al., 2016).

To construct a matching estimator in model ((ref)) we first note that, by the threshold crossing representation, there exist a parameter vector $\mu=(\mu_1,...,\mu_{G-1})$ and a standard uniform random variable $V$ such that $Z_g=Z_g(V,\mu)$, where $Z_1(V,\mu)=1$ if and only if $V\leq \mu_1$, $Z_g(V,\mu)=1$ if and only if $\mu_{g-1}<V\leq \mu_g$ for $g=2,...,G-1$, and $Z_{G}(V,\mu)=1$ if and only if $\mu_{G-1}<V$. We denote as ${\cal{M}}_{G-1}$ the set of vectors $\mu\in \mathbb{R}^{G-1}$ such that $0\leq \mu_1\leq \mu_2\leq ...\leq \mu_{G-1}\leq 1$. We then define the following estimator:

where $V_1,...,V_N$ are standard uniform draws, and $\sigma_{gt}$ are random permutations in $\Pi_N$ for all $g=1,...,G$, $t=1,...,T$, all independent of each other.

For given $\mu$, we use an algorithm analogous to the one described in Section (ref) to compute $\widehat{X}$. The outer minimization with respect to $\mu$ can be performed using simulated annealing or other methods to minimize non-differentiable objective functions. In Appendix (ref) we report simulation results for a nonparametric two-component mixture model. In that case grid search is a viable option. Moreover, a similar approach can be used to estimate finite mixtures of linear independent factor models (also known as “mixtures of factor analyzers”); see Ghahramani and Hinton (1997) and McLachlan et al. (2003).

Heteroskedastic deconvolution

Finally, consider the model

where $(X_1,S)$ is independent of $X_2$, and $X_2\sim F$, where $F$ is known and has zero mean. The econometrician observes a sample $Y_1,\widetilde{S}_1,...,Y_N,\widetilde{S}_N$ from $(Y,\widetilde{S})$, where $\widetilde{S}_i$ is a consistent estimator of $S_i$ for all $i$.

To motivate this setup, consider the estimation of neighborhood effects on income in Chetty and Hendren (2018), where $i$ is a commuting zone or county, and $Y_i$ is a neighborhood-specific estimate of the “causal effect” of place $i$. Within-$i$, a central limit theorem-type argument suggests that $Y_i$ is approximately normally distributed, with mean $X_{i1}$ and standard deviation $S_i$. Chetty and Hendren report, alongside $Y_i$ estimates, standard deviation estimates $\widetilde{S}_i$. In this example $F$ is the standard normal distribution. See Azevedo et al. (2019) for other examples and a parametric estimation approach.

To estimate the distribution of $X_{1}$ by matching, we minimize the following objective:

where $\sigma$ is a random permutation of $\{1,...,N\}$, ${\cal{S}}_N$ is the parameter space for $S$, and $\lambda>0$ is a constant. In this case our algorithm again consists in alternating optimal transport steps and least squares steps. As an illustration, in Appendix (ref) we estimate the density of neighborhood effects across US commuting zones using the data from Chetty and Hendren (2018).

Conclusion

In this paper we have proposed an approach to nonparametrically estimate models with latent variables. The method is based on matching predicted values from the model to the empirical observations. We have provided a simple algorithm for computation, and established consistency. We have also documented remarkable performance of our nonparametric estimator in small samples, and we have used it to shed new light on the cyclicality of permanent and transitory shocks to income and wages in the US. Progress on computation might be possible by leveraging recent advances on regularized optimal transport (Cuturi, 2013; Peyr� and Cuturi, 2019). Finally, an important question for future work will be to characterize rates of convergence and asymptotically valid confidence sets for our estimator.

\baselineskip14pt

\baselineskip21pt