Bivariate Distribution Regression; Theory, Estimation and an Application to Intergenerational Mobility

Bivariate distribution regression; theory, estimation and an application to intergenerational mobility

\onehalfspacing

\thispagestyle{empty}

\sloppy

We employ distribution regression (DR) to estimate the joint distribution of two outcome variables conditional on chosen covariates. While Bivariate Distribution Regression (BDR) is useful in a variety of settings, it is particularly valuable when some dependence between the outcomes persists after accounting for the impact of the covariates. Our analysis relies on a result from Chernozhukov2018distribution which shows that any conditional joint distribution has a local Gaussian representation. We describe how BDR can be implemented and present some associated functionals of interest. As modeling the unexplained dependence is a key feature of BDR, we focus on functionals related to this dependence. We decompose the difference between the joint distributions for different groups into composition, marginal and sorting effects. We provide a similar decomposition for the transition matrices which describe how location in the distribution in one of the outcomes is associated with location in the other. Our theoretical contributions are the derivation of the properties of these estimated functionals and appropriate procedures for inference. Our empirical illustration focuses on intergenerational mobility. Using the Panel Survey of Income Dynamics data, we model the joint distribution of parents' and children's earnings. By comparing the observed distribution with constructed counterfactuals, we isolate the impact of observable and unobservable factors on the observed joint distribution. We also evaluate the forces responsible for the difference between the transition matrices of sons' and daughters'.

\noindentKeywords: Bivariate distribution regression; joint distribution; local Gaussian correlation, Decomposition, intergenerational mobility

\noindentJEL: C14, C21

\onehalfspacing

\onehalfspacing

Introduction

Although distribution regression (DR) was introduced in the 1970s by Williams1971, and revisited by Foresi1995 in the 1990s, it has only relatively recently been more actively employed in empirical economic investigations. For example, inference for DR with a continuous outcome variable was only developed in Chernozhukov2013inference and for discrete, mixed discrete or continuous outcomes in Chernozhukovs2020. These papers have resulted in an examination of the economically interesting functionals which can be obtained from an estimated distribution of an economic variable conditional on a vector of exogenous variables. This has resulted in the use of DR in a range of different models and settings. These include, for example, models with endogeneity Chernozhukov2020nonseparable, sample selection Chernozhukov2018distribution, FVV2024, dynamic panels fernandezval2023dynamic, and most recently difference-in-difference estimation fernandezval2024.

While these extensions cover a wide class of models they, with the exception of some studies we note below, are restricted to models with univariate outcomes. However, for many empirical questions it is necessary to examine how the joint distribution of two outcomes varies in response to changes in the conditioning variables. For example, in investigating household labor supply responses to changes in tax rates, one should quantify how the hours worked by each of the spouses respond and, perhaps more importantly, how changes in their individual labor supplies are correlated, thereby affecting overall household labor supply. Alternatively, in examining the impact of changes in the minimum wage it would be useful to investigate how both wages and hours of work individually and jointly respond. These examples are just illustrative of the large number of economic questions in which it is important to examine more than one outcome and, perhaps more importantly, how the joint distribution of the outcomes is affected.

The estimation of models with multiple outcomes by DR has been considered in Meier2020 and Wang2022. Meier2020 extends univariate DR to the multivariate setting by using a multivariate indicator function at every location of the distribution. Wang2022 propose a fast factorization method that also employs univariate DR. However, this transformation to a univariate approach is potentially restrictive as it requires a single outcome variable to be constructed from the multiple outcomes of interest. As DR then employs this constructed outcome, it does not enable each of the individual outcomes to be modeled separately. Moreover, it also bypasses the modeling of the correlation between the unobservables driving each of the outcome variables. In contrast, BDR enables each of the individual outcomes to be explicitly modeled and facilitates an analysis of the local dependence structure. Furthermore, the choice of the link function in Meier2020 and Wang2022 might provide a poor approximation to the underlying distribution. In contrast, a result from Chernozhukov2018distribution, guarantees that the joint distribution is locally Gaussian and BDR is able to approximate it.

We employ DR to investigate settings in which the outcome of interest is the joint distribution of two outcome variables conditional on a set of covariates. Moreover, we describe how bivariate distribution regression (BDR) can be implemented, as well as the potentially interesting functionals that can be derived from BDR estimates. We separately analyze the roles of the marginal distributions and the dependence structure. As each is modeled separately in BDR, a high degree of flexibility is incorporated between the role of both observable and unobservable factors. Alternative approaches to modeling conditional joint distributions exist. These include the use of copulas Klein2022 and non-parametric estimators Bouzebda2019. Due to its semi-parametric properties, BDR can be seen as a mixture between these approaches. In contrast to copula models, BDR allows for greater flexibility in incorporating the impact of the covariates and requires relatively weaker parametric assumptions. Non-parametric approaches are more likely to suffer from the curse of dimensionality. A detailed overview of modeling choices for multivariate distributions can be found in Meier2020.

Bivariate or multivariate distribution regression has been employed elsewhere. For example, fernandezval2023 model the joint distribution of spouses' wages in the presence of selection rules, fernandezval2024 employ it in implementing difference-in-difference procedures with bivariate outcomes, while fernandezval2022cpscouples employ it to describe the joint distribution of spouses' annual earnings. The first two papers provide a theoretical contribution related to the literature on which they are focused and the third is purely a descriptive analysis. This paper provides important original theoretical contributions related to inference, the validity of the bootstrap employed, and the decompositions of the joint distribution and transition matrices.

Our empirical illustration is related to intergenerational mobility. This is a research area with potentially important policy implications on a range of topics such as tax policy and inequality reduction. It has received significant attention from policymakers and has generated a substantial academic literature. In a seminal paper, Chetty2014 showed that schools, neighborhoods, and family stability are influential determinants of mobility. Others have shown that race Chetty2020, immigration status Abramitzky2021, human capital Adermon2021, or wealth Adermon2018 are correlated with various mobility measures (see Mogstad2021 for a recent overview). This suggests that one should include many covariates when modeling intergenerational mobility noting that the appropriate manner to do so is not straightforward given that the relationship between the covariates and earnings appears to vary over the earnings' distributions. Previous studies typically present estimates of relative mobility across subsamples of the data or employ a linear regression approach in which the child's outcomes (rank of income, human capital, or wealth) are regressed on the parents' analog and (sometimes) covariates. By employing BDR to model the joint distribution of fathers' and children's income we provide an additional methodological approach for evaluating intergenerational mobility. First, we allow for the covariates to have varying effects at different points of the earnings' distributions. Second, we can allow for different covariates to have an impact on each of the marginal distributions and the dependence structure. Third, as the joint distribution captures the entire dependence structure between the respective income measures, standard mobility measures, such as the rank-rank relationship and the conditional expected rank, can be directly derived from the BDR estimates. Fourth, we can explore how these mobility measures are influenced by the dependence in the unobservables and how this dependence may vary across the income distribution. Finally, BDR estimates can be employed to produce counterfactual bivariate distributions corresponding to scenarios capturing alternative values of the local correlation.

The remainder of this paper is organized as follows. Section (ref) formally presents our methodological contribution. In Section (ref) we present a number of functionals which are potentially of interest in empirical work. Section (ref) outlines the associated estimation procedure. Section (ref) presents the asymptotic theory. Section (ref) contains our empirical application. Section (ref) concludes.

Bivariate Distribution Regression

Let $F_{Y,W \mid X}$ be the joint distribution of $(Y,W)$ conditional on $X$. We consider the bivariate distribution regression (BDR) model:

where $\Phi_{2}(\cdot, \cdot; \rho)$ is the distribution of the standard bivariate normal with parameter $\rho$, and $u \mapsto g(u)$ is a known link function with range $[-1,1]$ such as the Fisher transformation $g(u) = \tanh(u)$. In (ref) the marginal distributions of $Y$ and $W$ conditional on $X$ follow distribution regression models: $$ F_{Y \mid X}(y \mid x) = \Phi(x'\mu_y), \quad F_{W\mid X}(w \mid x) = \Phi(x'\nu_w), $$ where $\Phi$ is the distribution of the standard normal. The function $(y,w,x) \mapsto g(x'\delta_{yw})$ measures the local dependence or sorting between $Y$ and $W$ at $(Y,W,X) = (y,w,x)$. This sorting can vary with respect to observed covariates $X$, and along the distribution as indexed by $(y,w)$. The simplest case is $g(x'\delta_{yw}) = g(\delta_{yw})$, where the sorting only varies with respect to unobservables.

The BDR model can be motivated by the local Gaussian representation (LGR) of Chernozhukov2018distribution, which establishes that for any conditional joint distribution,

for some functions $\mu$, $\nu$ and $\rho$. The LGR is the right-hand-side of the previous equation and is unique, that is there is a one-to-one mapping between a conditional joint distribution and its LGR. In the LGR the marginal distributions of $Y$ and $W$ conditional on $X$ are represented by nonparametric distribution regression models, that is $$ F_{Y \mid X}(y \mid x) = \Phi(\mu(y \mid x)), \quad F_{W \mid X}(w \mid x) = \Phi(\nu(w \mid x)). $$ The BDR model can be seen as a semiparametric specification for the LGR where the nonparametric functions $\mu$, $\nu$ and $\rho$ are replaced by (generalized) linear indexes with function-valued parameters. In particular, the parameters $y \mapsto \mu_y$ and $w \mapsto \nu_w$ measure the effect of the covariates on the marginal distributions of $Y$ and $W$, and $(y,w) \mapsto \rho_{yw}$ measures the effect of the covariance on the local dependence (copula) between $Y$ and $W$.

Let $\bar \mathcal{Y}$ and $\bar \mathcal{W}$ denote strict subsets of $\mathcal{Y}$ and $\mathcal{W}$, the supports of $Y$ and $W$, respectively. We impose the following restrictions on the coefficients at the tails to facilitate estimation and inference, $$ \mu_{y,1} = \mu_{\bar y_y, 1} + (y - \bar y_y)\alpha_{\bar y_y}, \quad \mu_{y,-1} = \mu_{\bar y_y,-1},\quad y \in \mathcal{Y} \setminus\bar{\mathcal{Y}}, $$ $$ \nu_{w,1} = \nu_{\bar w_w,1} + (w - \bar w_w)\alpha_{\bar w_w}, \quad \nu_{w,-1} = \nu_{\bar w_w,-1},\quad w \in \mathcal{W} \setminus\bar{\mathcal{W}}, $$ and $$ \delta_{y,w} = \delta_{\bar y_y \bar w_w},\quad y \in \mathcal{Y} \setminus\bar{\mathcal{Y}}, \quad w \in \mathcal{W} \setminus\bar{\mathcal{W}}, $$ where $\bar y_y := \arg \min_{y' \in \bar{\mathcal{Y}}} |y-y'|$, $\alpha_{\bar y_y} > 0$, $\bar w_w := \arg \min_{w' \in \bar{\mathcal{W}}} |w-w'|$, $\alpha_{\bar w_w} > 0$, and the subscript $1$ and $-1$ denote the first element of a vector and its complement. For example, $\mu_{y,1}$ is the first element of $\mu_y$ and $\mu_{y,-1}$ is a vector with the rest of the elements. Here, we follow chernozhukov2025 and postulate that the random variables $Y$ and $W$ behave in the tails like random variables with distribution $\Phi$, after subtracting the location shifts $x'\mu_{\bar y}$ and $x'\nu_{\bar w}$, and dividing by the scales $\alpha_{\bar y_y}$ and $\alpha_{\bar w_w}$, which are different at the upper and lower tails. For the local dependence parameter $\delta_{yw}$, we postulate that it is constant at the tails. We could allow for additional tail parameters for the intercept of $\delta_{yw}$, but we set them to zero for simplicity.\footnote{Note that, unlike the intercepts of $\mu_y$ and $\nu_w$, the intercept of $\delta_{yw}$ does not need to satisfy any monotonicity restriction at the tails.}

Functionals of Interest

Rather than examine statistics which can be obtained in the univariate setting, we exploit the rich information contained in the conditional bivariate distribution and focus on functionals which highlight the dependence structure in the data. Accordingly, the local dependence parameter $\delta_{yw}$ features in each of the objects that follow.

Joint Distribution and Decomposition

The joint distribution of $Y$ and $W$ can be represented in terms of the BDR model as $$ F_{Y,W}(y,w) = \int \Phi_2(x'\mu_y,x'\nu_w;g(x'\delta_{yw})) {\mathrm{d}} F_X(x), $$ where $F_X$ is the distribution of $X$. Using this representation, we can decompose the difference between the joint distribution of two groups, say $D=1$ and $D=0$. To do so, it is convenient to introduce counterfactual distributions that combine BDR parameters and distributions of $X$ from different groups. Let

where $\mu_y^i$ is the parameter $\mu_y$ in group $i$ and the rest of the terms are defined analogously. In this notation the actual distribution in group $d$ is $F^{(d,d,d,d)}_{Y,W}$.

The difference in the actual distributions between groups $1$ and $0$ can then be decomposed as

where the first term in brackets is a composition effect, the second term is a dependence or sorting effect, and the third term is the effect of the difference in the marginal distributions. The last term can be further decomposed into

where the first term is the marginal effect related to $W$ and the last term is the marginal effect related to $W$.

Transition Matrices and Decomposition

We can also construct transition matrices and decompose differences across those for different groups into composition, sorting and marginal distribution components using counterfactual distributions. Let $\{y_j: 0 \leqslant j \leqslant J\}$ and $\{w_k: 0 \leqslant k \leqslant K\}$ be grids of values covering the supports of $Y$ and $W$, where we set $y_0=w_0 = - \infty$ and $y_K = w_K = +\infty$. Then, the $(j,k)$ element of the counterfactual transition matrix $T^{(k,d,r,s)}$ is

for $j = 1,\ldots,J$ and $k = 1,\ldots,K$. These matrices provide a parsimonious representation of the joint distribution of $Y$ and $W$.

The difference in the actual transition matrices between groups $1$ and $0$ can then be decomposed as

We provide examples of this decomposition in Section (ref).

Estimation

Rather than employing probit, as is done in univariate DR, BDR estimates a sequence of bivariate probits. Let $I^y := 1(Y \leqslant y)$, $\bar I^y := 1 - I^y$, $J^w := (W \leqslant w)$ and $\bar J^w := 1 - J^w$, for $(y,w) \in \mathcal{Y}_n \times \mathcal{W}_n$, where $\mathcal{Y}_n$ and $\mathcal{W}_n$ are grids of points in the supports of $Y$ and $W$, respectively.

Under the BDR model, the joint probability function of $I^y$ and $J^w$ conditional on $X=x$ is

Given a random sample of size $n$, $\{(Y_i,W_i,X_i)\}_{=1}^n$, the average conditional log-likelihood at $(y,w)$ is

The conditional maximum likelihood estimator (CMLE) of $(\mu_y,\nu_w, \delta_{yw})$ is

Several remarks should be noted regarding computation. First, when $\rho(x'\delta_{yw}) = \rho(\delta_{yw})$, then the CMLE is the bivariate probit estimator of $(I^y,J^w)$ on $X$. Second, the CMLE can be computed in two steps. In the first step, we obtain the estimators of $\mu_y$ and $\nu_w$ by probit distribution regression of $I^y$ on $X$ and $J^w$ on $X$, respectively. In the second step, we plug in the estimators of $\mu_y$ and $\nu_w$ from the first step in the average log-likelihood and maximize with respect to $\delta$, that is

where $\widehat \mu_y$ and $\widehat \nu_w$ are the distribution regression estimators of $\mu_y$ and $\nu_w$. Finally, to trace the joint and marginal distributions of $Y$ and $W$, we need to obtain the CMLE for a grid of values of $(y,w)$ in the support of $(Y,W)$. For example, we can use the tensor product of two grids for $Y$ and $W$, including sampling percentiles from the $1$ or $2\%$ to the $99$ or $98\%$.

Let ${\mathbb{E}_n}$ denote the empirical expectation, that is ${\mathbb{E}_n} := n^{-1} \sum_{i=1}$.

The counterfactual terms of the decomposition of the joint distribution can be estimated via a plug-in rule. For example, an estimator of the second integral on the right-hand side of (ref) can be constructed as follows: $$ \frac{1}{n_0} \sum_{i=1}^n (1-D_i) \Phi_2(X_i'\widehat \mu^1_y,X_i'\widehat \nu^1_w;g(X_i'\widehat \delta^1_{yw})), $$ where $\widehat \mu^1_y, \widehat \nu^1_w, \widehat \delta^1_{yw}$ are estimators of $\mu^1_y, \nu^1_w, \delta^1_{yw}$ in the subsample with $D=1$ and $n_0 := \sum_{i=1}^n (1-D_i)/n$. \footnote{Alternatively, a doubly-robust estimator can be formed as $$ \frac{1}{n_0} \sum_{i=1}^n \frac{(1-D_i)}{1-\widehat p(X_i)} 1\{Y_i \leqslant y, W_i \leqslant w\} - \frac{1}{n_0} \sum_{i=1}^n \frac{(D_i-\widehat p(X_i))}{1-\widehat p(X_i)} \Phi_2(X_i'\widehat \mu^1_y,X_i'\widehat \nu^1_w;g(X_i'\widehat \delta^1_{yw})), $$ where $\widehat p(x)$ is an estimator of the propensity score $p(x) = {\mathrm{P}}(D=1 \mid X=x)$.}

Inference

We use the exchangeable bootstrap for inference. The following algorithm describes how to obtain counterfactual draws of the BDR estimator.

Exchangeable bootstrap draws of the estimators of the functionals of interest can be obtained similarly.

Asymptotic Theory

To simplify the expressions of some of the assumptions and limit processes, it is convenient to introduce the following notation for the tails. Let $\bar r_r := \arg \min_{r' \in \bar{\mathcal{R}}} |r-r'|$ for $r \in \{y,w\}$ and $\bar \mathcal{R} \in \{\bar \mathcal{Y},\bar \mathcal{W}\}$. Note that $\bar r_r = r$ if $r \in \bar \mathcal{R}$, $\bar r_r = \bar r$ if $r \ge \bar r$ and $\bar r_r = \underline r$ if $r \leqslant \underline r$, where $\bar r := \sup (\bar \mathcal{R})$ and $\underline r := \inf(\bar \mathcal{R})$. We use a superscript $d$ to denote variables, parameters and samples from the group $d \in \{0,1\}$. We also use the following notation for partial derivatives: $\partial_x f(x) := \partial f(x)/\partial x$ and $\partial_{x x} f(x) := \partial^2 f(x)/(\partial x \partial x')$.

Assumption (ref) implicitly imposes common support for the variables $Y^d$, $W^d$ and $X^d$ across $d$. This condition guarantees that all the counterfactual distributions that we consider are well-defined.

Let $\widehat \theta^d_{yw} := ( \widehat \mu_y^{d'}, \widehat \nu_w^{d'}, \widehat \delta_{yw}^{d'})'$ and $\widehat \theta^{d*}_{yw} := ( \widehat \mu^{d*'}_y, \widehat \nu^{d*'}_w, \widehat \delta^{d*'}_{yw})'$ be the estimator and corresponding bootstrap draw of the parameter $\theta^d_{yw} := ( \mu_y^{d'},\nu_w^{d'},\delta_{yw}^{d'})'$, for $y \in \bar \mathcal{Y}$ and $w \in \bar \mathcal{W}$. Similarly, let $\widehat \xi_y^d := ( \widehat \alpha_y^d,\widehat \mu_y^{d'})'$ and $\widehat \xi_w^d := ( \widehat \alpha_w^d, \widehat \nu_w^{d'})'$ be the estimators of the tail parameters $\xi_y^d := ( \alpha_y^d, \mu_y^{d'})'$ and $\xi_w^d := ( \alpha_w^d, \nu_w^{d'})'$, for $y \in\{\underline{y} , \bar y\}$ and $w \in\{\underline{w} , \bar w\}$, and $\widehat \xi^{d*}_y := ( \widehat \alpha^{d*}_y,\widehat \mu^{d*'}_y)'$ and $\widehat \xi^{d*}_w := ( \widehat \alpha^{d*}_w, \widehat \nu^{d*'}_w)'$ be the corresponding bootstrap draws of the estimators. Let $Z_n \leadsto Z $ in $\mathbb{D}$ denote weak convergence of a stochastic process $Z_n$ to a random element $Z$ in a normed space $\mathbb{D }$, as defined in van1996weak. For example, $\mathbb{D } = \ell^{\infty}(\bar \mathcal{Y} \bar \mathcal{W})$, the set of measurable and bounded functions on $\bar \mathcal{Y} \bar \mathcal{W}$.

In order to state the results about bootstrap validity formally, we follow the notation and definitions in van1996weak. Let $D_{n}$ denote the data vector and $M_{n}$ be the vector of random variables used to generate bootstrap draws given $D_{n}$. Consider the random element $ \mathbb{Z}^{*}_{n} = \mathbb{Z}_{n}(D_{n}, M_{n})$ in a normed space $ \mathbb{E}$. We say that the bootstrap law of $\mathbb{Z}^{*}_{n}$ consistently estimates the law of some tight random element $\mathbb{Z}$ and write $\mathbb{Z}^{*}_{n} \leadsto_{{\mathrm{P}}} \mathbb{Z} $ in $\mathbb{D}$ if:

where $\text{BL}_{1}(\mathbb{D})$ denotes the space of functions with Lipschitz norm at most 1 and ${\mathrm{E}}_{M_{n}}$ denotes the conditional expectation with respect to $M_{n}$ given the data $D_{n}$; and $ \rightarrow_{{\mathrm{P}}}$ denotes convergence in (outer) probability.

The quantities of interest are functionals of the joint distribution of $(Y,W)$ conditional on $X$ and the marginal distribution of $X$. To state the results about these functionals it is convenient to introduce some notation. The following definitions apply to all $j,k,l,m \in \{0,1\}$. Let $F^m_X$, $\widehat F^m_X$ and $\widehat F^{m*}_X$ be the joint distribution of $X^m$, empirical counterpart and bootstrap draw; and $F^{(jkl)}_{ywx} := \Phi_2(x'\mu^j_y,x'\nu^k_w;g(x'\delta^l_{yw}))$, $\widehat F^{(j,k,l)}_{ywx} := \Phi_2(x'\widehat \mu^j_y,x'\widehat \nu^k_w;g(x'\widehat \delta^l_{yw}))$ and $\widehat F^{(jkl)*}_{ywx} := \Phi_2(x'\widehat \mu^{j*}_y,x'\widehat \nu^{k*}_w;g(x'\widehat \delta^{l*}_{yw}))$ be the counterfactual conditional distribution, its estimator and bootstrap draw. Consider the empirical processes and their bootstrap draws, $(y,w,x) \mapsto \widehat Z^{(jkl)}_{ywx} := \sqrt{n}(\widehat F^{(jkl)}_{ywx} - F^{(jkl)}_{ywx})$, $(y,w,x) \mapsto \widehat Z^{(jkl)*}_{ywx} := \sqrt{n}(\widehat F^{(jkl)*}_{ywx} - \widehat F^{(jkl)}_{ywx})$, $f \mapsto \widehat G^m_X(f) := \sqrt{n} \int f \mathrm{d} (\widehat F^m_X - F^m_X)$ and $f \mapsto \widehat G^{m*}_X(f) := \sqrt{n} \int f \mathrm{d} (\widehat F^{m*}_X - \widehat F^m_X)$, for $f \in \mathcal{F}$, where $\mathcal{F}$ is a class of measurable functions that (i) includes $F^{(jkl)}_{ywx}$, the indicators of all the rectangles in $\bar{\mathbb{R}}^{d_x+2}$, for $\bar{\mathbb{R}} = \mathbb{R} \cup \{-\infty,+\infty\},$ the extended real line, and (ii) is totally bounded under the metric: $$ \lambda(f; \tilde f) = \left[\int (f - \tilde f)\mathrm{d} F_X \right]^{1/2}, \quad f,\tilde f \in \mathcal{F}. $$ For $d \in \{0,1\}$, the selection matrix $\boldsymbol{S}_d^{(jkl)}$ picks up the elements of $Z_{ywx}^{d\theta}$ corresponding to the components used in the construction of the counterfactual conditional distribution. For example, $\boldsymbol{S}_0^{(010)} = {\rm diag}(\boldsymbol{I}_{d_{\mu}},\boldsymbol{0}_{d_{\nu}},\boldsymbol{I}_{d_{\delta}})$ and $\boldsymbol{S}_1^{(010)} = {\rm diag}(\boldsymbol{0}_{d_{\mu}},\boldsymbol{I}_{d_{\nu}},\boldsymbol{0}_{d_{\delta}})$, where $\boldsymbol{0}_{p}$ is a $p\times p$ matrix of zeros and $\boldsymbol{I}_{p}$ is the identity matrix of size $p$.

Let $$ \widehat F_{Y,W}^{(j,k,l,m)}(y,w) = \dfrac{1}{n_m} \sum_{i=1}^{n_m} \Phi_2(X_i^{m'}\widehat \mu_y^j, X_i^{m'}\widehat \nu_w^k;g(X_i^{m'}\widehat \delta^l_{yw})), $$ be an estimator of the counterfactual joint distribution in (ref) and $$ \widehat F_{Y,W}^{*(j,k,l,m)}(y,w) = \dfrac{1}{n_m} \sum_{i=1}^{n_m} \omega_{in_m}^m \Phi_2(X_i^{m'}\widehat \mu_y^{j*}, X_i^{m'}\widehat \nu_w^{k*};g(X_i^{m'}\widehat \delta^{l*}_{yw})), $$ be the corresponding bootstrap draw. The following result follows from Theorem (ref), Lemma D.1 of Chernozhukov2013inference, which establish the Hadamard differentiability of the counterfactual map, and the functional delta method.

Empirical Illustration

We now employ BDR to study intergenerational income mobility. BDR is a particularly useful and powerful tool for this area of research since it can model the joint distribution of incomes for different generations of family members. We focus on the joint distribution of a child’s and their father's labor incomes. We begin by documenting the raw data. This provides some insight into the presence of intergenerational mobility in these data. We then employ BDR to estimate the joint conditional distribution of these two outcomes conditional on covariates for both the child and the father. We examine the role of the estimated local dependence parameter in generating these observed outcomes. Finally, we decompose the observed differences in dependence between fathers/sons and fathers/daughters into composition, sorting and marginals effects.\footnote{We acknowledge the use of the term "sorting" in this potential empirical context is somewhat inappropriate given that one would not interpret the process by which fathers and children are paired is best characterized by a sorting process. However, we retain this description to avoid the introduction of additional terminology.}

Data

Our analysis examines data from the Panel Study of Income Dynamics (PSID), a rich longitudinal dataset including family and individual-level income data for individuals in the United States since 1968. The survey was conducted yearly until 1997 and biannually thereafter. The PSID has been employed in previous examiniations of intergenerational mobility (see, for example, Callaway2019 or landerso2017scandinavian).

Although the PSID data has 904,796 observations, the number for which we can construct father/children pairs is 6,6767 comprising 3,424 daughters and 3,243 sons. Our empirical analysis is based on this smaller sample. A pair is included if both the father and the child are observed at least once between the ages of 25 and 50. The outcome variables are their respective average annual labor earnings. This encompasses wages, business income, bonuses, and overtime payments. Labor earnings have been standardized to 1982 dollars. Note that we exclude observations with zero labor income earnings for that year. The average frequency of observations is the same for both father-daughter and father-son pairs. That is, 9 years for children and 15 years for fathers. The equal frequency for males and females is puzzling given the higher labor force participation rates of males. This might reflect the construction of the data set but probably also be due to our relatively weak requirement that the child only needs to be observed with earnings in one year. The higher average annual work hours for sons (2011) than daughters (1540) is consistent with males having a relatively higher level of labor market engagement than females.

The inclusion of only pairs for which both the father and child report labor income may introduce some form of selection bias. Moreover, not accounting for the frequency that the pair is observed may also introduce an additional form of bias. While each of these issues is interesting, we delay their treatment to future work as they are beyond the scope of this paper.

Studies in the literature rely on various income measures to study intergenerational mobility. These include, for example, labor earnings, tax records, or lifetime income imputed over the life cycle. Mogstad2021 note that the estimated degree of mobility can vary considerably depending on the measure used. We focus on labor income, noting that intergenerational persistence tends to be stronger for broader income concepts landerso2017scandinavian.

The variables displayed in Table (ref) are all available since 1968 with the exception of “college”, which indicates that the individual has a college degree, which is available from 1975. Observations with missing entries for any of the variables are excluded. Table (ref) indicates that the father's information is generally collected during the 1980s, while the children's information is typically gathered in the early 2000s. Fathers are observed later in life, mostly between the ages of 34 and 47, while children are observed earlier between the ages of 27 and 36. This also explains the large difference in labor income in favor of fathers. Unsurprisingly, we find higher levels of labor income for fathers with a college degree. The same is true for their children. This can be partially explained by the high level of persistence in college education. For example, 30 percent of sons have a college degree, but that number rises to 61 percent among sons whose fathers have a college degree. Similar figures are found for daughters. We also find that white fathers, and their children, generally have higher labor income. Household size is smaller for children, which can be attributed to their younger age, although it may also capture broader demographic changes and trends.

Unconditional Mobility

Figure (ref) plots the transition matrices for the outcome variables. The cells of the transition matrices are based on the quintiles of the fathers in USD on the y-axis and the quintiles of the children in USD on the x-axis. The values within these cells represent the percentage of observations observed with a father-daughter or father-son combination within these cells. Due to the overall wage gap between women and men, sons earn more than daughters and the percentages of sons are relatively higher on the right side of the figure, while they are relatively higher on the left side of the figure for daughters. As relative darkness of the squares capture larger probabilities, it is immediately apparent that father-daughter pairs are relatively more frequently located in the upper left corner and the father-son pairs are more relatively frequent in the lower right corner. A high degree of mobility would imply equal probabilities across all cells within each row. Therefore, the observed deviation from this pattern illustrates the limited mobility in the United States.

BDR Estimation

We now focus on the estimates of the BDR model. We first estimate the model using the covariates described below for each of the marginals while also allowing $\rho$ to be a function of the covariates. We also estimate the model in which $\rho$ is not a function of the covariates although these latter results are relegated to the appendix.

The empirical results reported in the main text uses the following variables. For the marginal distribution of father's labor income we employ age, age$^2$ and decade fixed effects for both fathers and children, father's and child's race, father's and child's household size and indicators for a college degree of the father and the child's origin (measuring the degree of urbanization where the child was raised). We use the same set of variables for the marginal distribution of the children. When $\rho$ is a function of the covariates we employ all of the covariates listed above but exclude the decade fixed effects. We acknowledge that some covariates are potentially endogenous. Accounting for this endogeneity is beyond the scope of this paper and we defer this important issue to further work. However, in the appendix we provide the corresponding results when we exclude some variables from the conditioning set. We note that their exclusion does not have substantial implication for the results reported.

Model predictions and Counter Factual

We first report the capacity of the BDR estimates to reproduce the observed joint distribution of fathers' and children's earnings and examine the impact on this predicted joint distribution when we eliminate the dependence operating through $\rho$. Figure (ref) presents the estimated transition matrices based on the BDR model. To construct these matrices, we estimate the joint distribution of father-child earnings using the set of covariates listed above. We then integrate over the empirical distribution of the observed covariates to obtain the average joint distribution.

Panels (A) and (B) display the resulting transition matrices for daughters and sons, respectively. These matrices correspond to the estimated joint distribution. The values of these matrices are calculated based on equation (ref). Panels (C) and (D) show the same average distribution, but under the hypothetical assumption that the local correlation is set to zero. This counterfactual allows us to isolate the contribution of local dependence to intergenerational mobility. Panels (E) and (F) present the difference between the two distributions and quantifies the effect of local dependence.

It is useful to consider two issues before focusing on the results from this exercise. The first is the interpretation of $\rho$ in this context and the second is the resulting estimates of $\rho$. Given the nature of the outcomes, the parameter $\rho$ captures how the outcomes move together after conditioning out the impact of the covariates. If the estimated value of $\rho$ was zero, this implies no dependence between the outcomes after conditioning out the covariates. If the estimated $\rho$ is positive (negative) this implies that, conditional on the covariates, a higher value for one outcome is associated with a higher (lower) outcome for the other. In this particular context this dependence may capture some local correlation between unobservables affecting the respective outcomes.

As noted above an attractive feature of BDR is its capacity to flexibly model the local dependence. Although we do not report our estimates of $\rho$, a notable feature of the results is the substantial heterogeneity in these estimates across different points on the grid. This reflects that the level of dependence varies greatly at different points of the joint distribution. This highlights the importance of adopting a flexible approach to its estimation.

Several insights emerge from these results. First, the estimated transition matrices in Panels (A) and (B) closely resemble those obtained from the raw data (Figure (ref)), indicating that the model successfully captures the observed patterns of mobility. While this is reassuring it is not surprising given the nature of the estimator and the large number of parameters estimated. Second, setting the local correlation parameter ($\rho$) to zero leads to a considerable flattening of the dependence structure. This effect is evident for both sons and daughters and is particularly pronounced in the upper tail of the income distribution, where the level of dependence is relatively strong in the raw data. The influence of local dependence is quantitatively meaningful. For example, the probability that a son reaches the top decile is 6.9% under local dependence, compared to 5.8% under the counterfactual with $\rho = 0$. This difference is economically significant and explains a large fraction of the departure from the benchmark probability of 4%, which would be expected under full independence and in the absence of gender wage gaps and compositional effects. This highlights the role of local correlation in capturing some aspect of the observed persistent intergenerational advantages. Comparable patterns are observed for daughters, for whom the impact of $\rho$ is most pronounced in the upper-left region of the figure.

Note that we also estimated the model predictions and the counterfactuals with a specification that does not take account of the covariates in $\rho$. These results are presented in Figure (ref). The results are very similar to the results presented in Figure (ref). This suggests that the results in this particulalr context are not sensitive to the treatment of the local dependence.

Decomposition of Joint Conditional CDF

Figure (ref) presents a decomposition of the differences in transition matrices between sons and daughters. The decomposition follows equation (ref) and separates the overall difference in the transition matrices into three distinct components: (i) Composition effects, which capture differences in observable characteristics (e.g. education, region, cohort); (ii) sorting effects, reflecting differences in the estimated local dependence parameter $\rho$; and (iii) marginal effects, which account for gender-specific differences in how covariates relate to the marginal earnings distributions.

Panels (A) and (B) display the estimated transition matrices for sons and daughters, respectively, using the BDR model with the set of covariates as described in Section (ref). Panel (C) shows the element-wise difference between the two matrices. Panels (D), (E), and (F) then report how much of this difference can be explained by composition, sorting, and marginal effects, respectively, expressed as a percentage of the absolute difference in each cell.\footnote{Note that the difference in Panel (C) is, in some instances, very close to zero, which means that the composition, sorting, and marginal effects are divided by a number near zero. This partially explains the occasional occurrence of very large values in Panels (D)–(F).}

Several key insights emerge from this analysis. Marginal effects are by far the most important driver of the observed gender differences in intergenerational earnings mobility. In many cells, they account for the majority of the variation across gender, highlighting the importance of differential returns to characteristics in shaping the marginal earnings distributions of sons and daughters. These marginal effects capture the differences in the labor market experiences, reflecting both supply and demand factors, for females and males. Sorting effects, driven by differences in the estimated local correlation, explain roughly 10% of the overall difference. This suggests that sons and daughters differ not only in their marginal earnings prospects, but also in how tightly their outcomes are linked to their fathers’ incomes, even after controlling for observed characteristics. Composition effects contribute very little to the observed differences. This indicates that the sons and daughters in our sample are, on average, very similar in terms of the included covariates. This suggests that the large difference in intergenerational earnings mobility for daughters and sons cannot be explained by differences in observable characteristics. These results emphasize that gender differences in intergenerational earnings mobility arise primarily from structural differences in earnings determination. These arise through the impact of observables on earnings and the nature of dependence between father's and child's outcomes.

We also provide a robustness check that excludes the children's race and household size. The results are presented in Table (ref). We find no substantial deviations from the results shown in Figure (ref). While we do not consider this an exhaustive examination of the sensitivity of our results, it does not provide some indication of their sensitivity to the exclusion of variables which are potentially important.

Discussion

Our empirical findings highlight the importance of modeling both the marginal distributions and the local dependence structure when analyzing intergenerational mobility. The BDR framework enables a flexible and interpretable decomposition of mobility, and highlights that both observable characteristics and residual local association patterns contribute to the overall structure. Importantly, one can isolate the channels through which gender differences in income dependence emerge. The stronger association observed in father-son pairs is largely due to the differences in the marginal distributions. However, accounting for differences in sorting, capturing local dependence, is also important. Overall, the evidence suggests that policies aiming to promote intergenerational equity should consider not just differences in endowments, but also structural differences in how income is transmitted across generations.

Conclusion

We employ bivariate distribution regression to estimate the joint distribution of two outcome variables conditional on set(s) of covariates. Our approach is particularly valuable when the dependence between the outcomes persists even after controlling for observed covariates. We describe how BDR can be implemented and present some associated functionals of interest. The functionals of primary interest in this setting are those involving the unexplained dependence. We provide a decomposition of the joint distributions for different groups into composition, marginal and sorting effects. We provide algorithms for estimation and inference, along with the associated theory. We also provide a similar decomposition for the associated transition matrices Using data from the Panel Survey of Income Dynamics, we model the joint distribution of parents' and children's earnings and find that the local dependence structure is an important contributor in explaining intergenerational mobility.