Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
115,553 characters · 34 sections · 78 citation commands
Dynamic Heterogeneous Distribution Regression Panel Models, with an Application to Labor Income Processes$^*$
\global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long\global\long \global\long \global\long \global\long\global\long \global\long \global\long \global\long \global\long \global\long \global\long\global\long\global\long \global\long \global\long\def\mc#1{\mathscr{#1}}
\global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long\def\abs#1{\left|#1\right|}
\global\long\def\norm#1{\left\Vert #1\right\Vert }
\global\long\def\rest#1{\left.#1\right|}
\global\long\def\bracket#1#2{\left\langle #1\middle\vert#2\right\rangle }
\global\long\def\sandvich#1#2#3{\left\langle #1\middle\vert#2\middle\vert#3\right\rangle }
\global\long\def\turd#1{\frac{#1}{3}}
\global\long \global\long\def\sand#1{\left\lceil #1\right\vert }
\global\long\def\wich#1{\left\vert #1\right\rfloor }
\global\long\def\sandwich#1#2#3{\left\lceil #1\middle\vert#2\middle\vert#3\right\rfloor }
\global\long\def\abs#1{\left|#1\right|}
\global\long\def\norm#1{\left\Vert #1\right\Vert }
\global\long\def\rest#1{\left.#1\right|}
\global\long\def\inprod#1{\left\langle #1\right\rangle }
\global\long\def\ol#1{\overline{#1}}
\global\long\def\ul#1{#1}
\global\long\def\td#1{ \widetilde{#1}}
\global\long \global\long \global\long \global\long \global\long
Keywords: distribution regression, individual heterogeneity, panel data, uniform inference, labor income dynamics, incidental parameter problem, poverty traps
\onehalfspacing
Empirical studies increasingly feature analyses of data comprising repeated observations on the same or similar units. While the most typical example is panel data, many of its attractive features are shared by other data structures, such as clustered, network and spatial data. A common feature of the panel data literature is its somewhat limited treatment of parameter heterogeneity browning2010heterogeneity. Although random coefficient panel models allow heterogeneous coefficients between units, and some recent developments discussed below incorporate heterogeneous coefficients within units, there are relatively few studies that incorporate heterogeneous coefficients both between and within units.\footnote{Exceptions include chetverikov2016iv, okui2019panel, zhang2019quantile and chen2021quantile, which are discussed in the literature review.} Allowing both forms of heterogeneity is potentially important in many economic settings.
This paper employs fixed effects distribution regression to estimate a dynamic panel model with coefficient heterogeneity between and within units. The model captures within-unit heterogeneous relationships between outcome and covariates through function-valued coefficients, and between-unit heterogeneity via coefficients which can vary across units in an unrestricted fashion. The model is very flexible and facilitates the analysis of various new functionals of the coefficients, including linear projections on unit covariates and predicted distributions. Via the estimated parameters we can create interesting counterfactual scenarios by manipulating the values of the initial conditions of the outcome variable, or the covariates, to examine the impact on these functionals. We can also consider both one-period-ahead and stationary counterfactual distributions to measure the short and long term effects of these changes.
The need for a flexible treatment of heterogeneity is partially motivated by previous empirical evidence highlighting the absence of consensus regarding the degree of heterogeneity required to appropriately model labor income dynamic processes. For example, abowd1989covariance and macurdy1982use considered models with a limited allowance for heterogeneity between units, whereas browning2007heterogeneity and browning2010modelling found some features of income processes, such as variances of the shocks, differ considerably between units. Moreover, browning2007heterogeneity and browning2010modelling noted that allowing for between unit heterogeneity has drastic implications for specific empirical questions. Others, including arellano2017earnings, allow for flexible heterogeneity within units to capture nonlinear persistence, but restrict the heterogeneity between units.
Modeling the correct degree of heterogeneity is also important for inference as the existing inference procedures might be only valid under specific circumstances. Here, rather than taking a firm position on the degree of heterogeneity, we provide a procedure that is uniformly valid over the degree of unknown heterogeneity of the coefficients. This covers “homogeneous models", “local heterogeneous models" (i.e. heterogeneity concentrated on subpopulations or parts of the distribution), and “strongly heterogeneous models" as special cases. The main theoretical challenge is that the rate of convergence of the estimators depends on the degree of heterogeneity, which is often unknown. This creates a uniform validity problem for the existing inference methods that require exact knowledge of the rate of convergence of the estimators. More formally, we establish that standard analytical plug-in methods are not valid for inference uniformly with respect to the degree of heterogeneity. We address this problem via a cross-sectional bootstrap scheme that resamples from the empirical distribution of the estimated coefficients. We show that this bootstrap is valid uniformly over various degrees of heterogeneity. A similar uniformity problem arises with the inference on average partial effects in nonlinear panel data models. fernandez2016individual, for example, bypassed this problem by assuming strong heterogeneity on the partial effects.
We also establish a relationship between dynamic distribution regression models with discrete outcomes and finite-state Markov chains.\footnote{browning2007heterogeneity,browning2010heterogeneity,browning2014dynamic previously established a related connection between dynamic binary response models with random coefficients and two-state Markov chains.} This allows us to express objects such as stationary distributions, mobility probabilities and recurrence times as functionals of the model coefficients.
Our methodology is applicable to a wide range of settings and we employ it here to examine labor income dynamics. This is an important research area with a large literature, starting with champernowne53, hart1976, shorrocks1976 and lillard1978dynamic, but now also including a long list of papers featuring econometric innovations. We apply our model to data from the Panel Study of Income Dynamics (PSID) to perform economically interesting experiments.
First, we consider how a ceteris paribus reduction in annual labor income in a given year, implemented via a negative shock, affects future annual labor income. We find that the predicted effect on the cross-sectional distribution of labor income after one period varies substantially depending on whether we account for heterogeneity in the level and persistence of income. Our model predicts substantially smaller effects than existing autoregressive models that restrict between and/or within heterogeneity by imposing several forms of homogeneous coefficients.
Second, we address the existence of poverty traps. We explicitly model the conditional probabilities of individuals being in poverty in a specific year given they were in poverty in the previous year. We establish that the substantial cross-sectional heterogeneity in the level and persistence of annual labor income has important implications for an individual's tendency to remain in a certain location of the income distribution.
From a theoretical perspective, our paper is related to chernozhukov2013inference (CFM) and chernozhukov2018network (CFW). The former studies distribution regression for cross-sectional data and the latter for panel data with fixed effects. Both flexibly model and estimate counterfactual distributions. We introduce two substantial and important departures from this earlier work. First, whereas all coefficients in CFM and CFW except for the intercept are fixed, we treat all coefficients as random. This facilitates the analysis of many economically interesting functionals which cannot be analyzed in the CFM and CFW frameworks. Moreover, our evidence below indicates this coefficient heterogeneity is empirically important to study labor income dynamics. It also introduces the theoretical challenge of how to perform inference that remains uniformly valid with respect to the degree of coefficient heterogeneity. These issues were not considered in CFM and CFW. Second, our model is dynamic, whereas those in CFM and CFW are static. This feature allows us to estimate economically interesting objects related to persistence.
Our model differs from the traditional random coefficients models of swamy1970, hsiao2008random, arellano2012identifying, fernandez2013panel and su2016, among others, as we allow for heterogeneous coefficients both between and within units. Moreover, existing distribution and quantile regression models with fixed effects often allow the intercepts to vary across units but restrict the slopes to be homogeneous; e.g., Koenker2004, Galvao2011, GalvaoKato2016smooth, KatoGalvaoMontes-Rojas2012, ArellanoWeidner2017, and chernozhukov2018network. chetverikov2016iv and chen2021quantile develop models similar to ours, but focus on projections of coefficients as the objects of interest in static quantile regression models. Other related recent works are okui2019panel and zhang2019quantile noting that their models and objects of interest differ from ours.
Bias correction methods based on large-$T$ asymptotic approximations for fixed effects estimators of dynamic and nonlinear panel models have been previously studied in Nickell81, Phillips:1999p733, Hahn:2004p882, fernandez2009fixed, HahnKuersteiner2011, dhaene2015split, and fernandez2016individual, among others (see ArellanoHahn2007 and fernandez2018fixed for recent reviews). We extend these debiasing methods to new functionals of the coefficients.
Inference that is robust to unknown heterogeneity has also been studied by liao2018unifddorm and lu2022uniform in the context of linear random coefficient panel models estimated by least squares. There are two main differences with respect to our context. First, our distribution regression model is nonlinear and its coefficients are estimated by conditional maximum likelihood methods. This introduces the additional difficulty of dealing with the incidental parameter problem. Second and more importantly, the coefficients of our model are infinity dimensional, whereas the coefficients in liao2018unifddorm and lu2022uniform are finite dimensional. This difference makes the inference problem more challenging as the coefficient functions in our model might exhibit different degrees of heterogeneity at different values of their domain. Similar to liao2018unifddorm, we propose cross-sectional panel bootstrap to make inference that is robust to the degree of coefficient heterogeneity. This method was previously used for panel data as a resampling scheme that preserves the dependence in the time series dimension, e.g., kapetanios2008bootstrap, kaffo2014bootstrap, and gonccalves2015bootstrap. We demonstrate that it also has robustness properties in models with heterogeneous coefficients.
Although our empirical work is related to the existing literature on labor income or earning processes, many aspects of our results are novel. This literature has typically focused on allocating the total error variances into transitory and permanent components. A summary is provided in moffitt2018income and three important recent innovations are arellano2017earnings,arellano2018nonlinear and hu2019semiparametric. The first two examined nonlinear persistence in the permanent component and how it varies over the earnings distribution. The third allowed for a flexible representation of the distributions of both components. Our approach is not intended to supersede these methodologies. Rather, we examine earnings dynamics to illustrate how our approach can complement the existing literature. The approach most similar to ours is arellano2017earnings,arellano2018nonlinear, which provided evidence of nonlinearity in income dynamics. While they considered nonlinear persistence that can vary by location in the earnings distribution, they do not allow heterogeneity between units. In contrast, we study income persistence that not only varies by location in the earnings distribution but also across units. Moreover, we allow persistence to be a function of both observed and unobserved individual characteristics and target different objects including counterfactual distributions and mobility probabilities.
The representation of the income mobility model as a finite-state Markov chain is motivated by champernowne53 and shorrocks1976, which previously used homogeneous Markov chains to analyze the same issue. We instead estimate a separate Markov chain for each unit to allow for unrestricted unit heterogeneity. We can recover the associated transition probabilities and apply standard tools for Markov chains to study stationary distributions and recurrence times. Other related work on income dynamics includes hirano2002 and gu2017, which estimated autoregressive labor income processes using flexible semiparametric Bayesian methods, and hoffmann2019hip, which studied the robustness of model parameters across specifications of the earnings dynamics. Finally, lee2025identificationestimationdynamicrandom studies identification and estimation of dynamic random coefficient models in short panels. In that setting, coefficients and their functionals are only partially identified when the time dimension is held fixed.
Let $\mathcal{F}_{it}$ be the filtration defined in Section (ref). We make use of several expectations. For a sequence of random variables $\{y_{it}: 1 \leq i \leq N, 1 \leq t \leq T\}$, where $i$ indexes cross-sectional units and $t$ time periods, we denote the expectation with respect to the distribution of $y_{it}$ conditional on $\mathcal{F}_{it}$ by $\mathbb E_{it} y_{it} := \mathbb E( y_{it} \mid \mathcal{F}_{it})$, the cross-sectional expectation at $t$ as $\mathbb E_{t} y_{it} := \text{plim}_{N \to \infty} N^{-1} \sum_{i=1}^N y_{it}$, and the cross-sectional and temporal expectation as $\mathbb E y_{it} := \text{plim}_{N,T \to \infty} (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T y_{it}$, provided that they exist. In what follows we shall assume without qualification that an expectation exists whenever it is used. For two deterministic sequences $a_T, b_T$, we use the notation $a_T\ll b_T$ and $b_T\gg a_T$ if $a_T=o(b_T)$.
The rest of the paper is organized as follows. Section (ref) presents the model and objects of interest. Section (ref) discusses estimation and inference. We present the empirical application in Section (ref) and Section (ref) establishes the associated asymptotic theory for our estimation and inference methods. Section (ref) reports simulation evidence. Proofs and additional results are gathered in the Appendix.
We observe a panel data set $\{(y_{it},\boldsymbol{x}_{it}) : 1 \leq i \leq N, 1 \leq t \leq T\}$, where $i$ typically indexes observational units and $t$ indexes time periods. The scalar variable $y_{it}$ represents the outcome or response of interest; and $\boldsymbol{x}_{it}$ is a $d_x$-vector of covariates, which includes a constant, lagged outcome values, and other predetermined covariates denoted by $\boldsymbol{v}_{it}$. That is, $$ \boldsymbol{x}_{it} =(1, y_{i(t-1)},...,y_{i(t-L)}, \boldsymbol{v}_{it}')'. $$
Let $\mathcal{F}_{it}$ be a filtration to which $\boldsymbol{x}_{it}$ and any time invariant variable for unit $i$ are adapted. We model the distribution of $y_{it}$ conditional on $\mathcal{F}_{it}$ as, for any $y \in \mathbb{R}$,
where $\Lambda: \mathbb R\mapsto[0,1]$ is a known, strictly increasing, and four times continuously differentiable link function (e.g., the standard normal or logistic CDF), and $y \mapsto - \boldsymbol{x}_{it}'\boldsymbol{\beta}_i(y)$ is increasing almost surely (a.s).\footnote{We can replace $\boldsymbol{x}_{it}$ by $P(\boldsymbol{x}_{it})$ where $P$ is a vector of transformations with good approximation properties such as polynomials, splines or interactions. Our asymptotic theory covers the case where the dimension of $P$ is fixed with the sample size. One could also allow the dimension of $P$ to grow with the sample size and $\Lambda$ to be unknown using semiparametric methods, but we do not pursue those extensions here.} This is a distribution regression model for panel data with heterogeneous coefficients, which we call a heterogeneous distribution regression model (HDR).
When $y_{it}$ is continuous, the model (ref) has the following representation as an implicit nonseparable model by the probability integral transform $$ \Lambda(-\boldsymbol{x}_{it}'\beta_i(y_{it})) = u_{it}, \ \ u_{it} \mid \mathcal{F}_{it} \sim U(0,1), $$ where the error $u_{it}$ represents the unobserved ranking of the observation $y_{it}$ in the distribution of unit $i$ at time $t$. The parameters of the model are related to derivatives of the conditional quantiles. Let $Q_{y_{it}}(u \mid \boldsymbol{x}_{it})$ be the $u$-quantile of $y_{it}$ conditional on $\boldsymbol{x}_{it}$ defined as the left-inverse of $y \mapsto F_{y_{it} }(y \mid \boldsymbol{x}_{it})$ at $u$, namely $$ Q_{y_{it} }(u \mid \boldsymbol{x}_{it}) := \inf\{y \in \mathbb R : F_{y_{it} }(y \mid \boldsymbol{x}_{it}) \geq u\}. $$ Then, it can be shown that if $y \mapsto F_{y_{it} }(y \mid \boldsymbol{x}_{it})$ is strictly increasing in the support of $y_{it}$ and $y \mapsto \boldsymbol{\beta}_i(y)$ is differentiable with $\dot \boldsymbol{\beta}(y) := \mathrm{d} \boldsymbol{\beta}(y)/\mathrm{d} y$,
The HDR coefficients are thus proportional to the derivatives of the conditional quantile function, and ratios of the DR coefficients correspond to ratios of the derivatives. The conditional quantile function $\boldsymbol{x}_{it} \mapsto Q_{y_{it} }(u \mid \boldsymbol{x}_{it})$ is therefore nonlinear and heterogeneous across $i$ and $u$.
By iterating expectations, the cross-sectional distribution of the observed outcome at time $t$ can be written in terms of the model coefficients as
This representation serves several purposes. First, it is the basis for a specification test of the model where an estimator of $F_{t}(y)$ based on the right hand side of (ref) is compared with the cross-sectional empirical distribution of $y_{it}$. Second, when $\boldsymbol{x}_{it}$ only includes lagged values of $y_{it}$, we can construct one-period-ahead predicted distributions by setting $t=T+1$. These distributions are useful for forecasting. Third, we can analyze dynamics of the distribution of $y_{it}$ over time. For example, below we analyze labor income mobility and the persistence of poverty traps. Fourth, we can consider the impact of interventions by comparing the counterfactual distribution after some intervention in $\boldsymbol{x}_{it}$ or $\boldsymbol{\beta}_i(y)$ with the actual distribution.
The main innovation of (ref) is that all model coefficients are random functions, $$ y \mapsto \boldsymbol{\beta}_i(y), $$ where the variation across $i$ and $y$ captures between-unit heterogeneity and within-unit heterogeneity, respectively. Moreover, we can explore if the two sources of heterogeneity are associated with observed unit characteristics using linear projections. Let $\boldsymbol{w}_{i}, \boldsymbol{z}_{i} \in \mathcal{F}_{i1}$ denote time invariant covariates such that $\dim(\boldsymbol{w}_{i})\geq \dim(\boldsymbol{z}_{i})$ and $\mathbb{E}(\boldsymbol{w}_{i} \boldsymbol{z}_{i}')$ has full column rank. Consider a linear regression
which covers the standard linear projection by setting $\boldsymbol{w}_{i} = \boldsymbol{z}_{i}$, and instrumental variables regression by setting $\boldsymbol{w}_i$ as the instrument. The coefficient $\boldsymbol{\theta}(y)$ is informative about which covariates are associated with the heterogeneity in $\boldsymbol{\beta}_i(y)$ across $i$, where we allow these relationships to vary within the distribution as indexed by $y$. For example, our empirical application explores whether initial income, education, race and year of birth are associated with differences in the level and persistence of labor income at different locations of the income distribution.
Our model can be used to construct flexible counterfactual distributions resulting from changing the values of the covariates and coefficients
where $h_{it}$ is a possibly data dependent transformation, and $ \boldsymbol{\beta}_i^g(y)$ is a transformation of the random coefficients. Specifically, we consider $$ \boldsymbol{\beta}_i^g(y) = \boldsymbol{\theta}(y)g(\boldsymbol{z}_{i}) + \boldsymbol{\gamma}_i(y) = \boldsymbol{\beta}_i(y) + \boldsymbol{\theta}(y)[g(\boldsymbol{z}_{i})-\boldsymbol{z}_{i}], $$ for a known transformation $g$ of the time invariant covariates $\boldsymbol{z}_i$. This transformation allows us to study the effect of changing the values of the covariates on the cross-sectional distribution through their impact on the random coefficients.
For instance, consider a hypothetical scenario where at time $t$ we increase the number of years of schooling to 12 for any worker who has less. If $\boldsymbol{z}_{i} = (z_{1i},\boldsymbol{z}_{-1,i}')'$, where $z_{1i}$ is the observed years of schooling of worker $i$ and $\boldsymbol{z}_{-1,i}$ includes the remaining components of $\boldsymbol{z}_i$, this counterfactual scenario is implemented via the transformation
$G_t(y)$ would then represent the counterfactual distribution of labor income at $t$ after the change. Another example is
which corresponds to giving an additional year of schooling to all workers.
We can also study the impact of shocks in dynamic models. For example, suppose at time $t-1$ a shock reduces income by $100*\kappa \%$ for individuals with income higher than certain (known) threshold $\tau_0$. This corresponds to the transformation
where $y_{i(t-1)}$ is measured in logarithmic scale. $G_t(y)$ now represents the counterfactual income distribution resulting from this income shock at time $t$.
When the actual and counterfactual cross-sectional distributions of $y_{it}$ are continuous, we consider quantiles of both distributions, and define quantile effects as their difference. Given a univariate distribution $F$, the quantile (left-inverse) operator is $$ \phi(F,\tau) := \inf\{y \in \mathbb{R}: F(y)\geq\tau\}, \quad \tau \in [0,1], $$ with the convenction $\inf \{ \emptyset \} = + \infty$. We apply this operator to the cross-sectional distributions defined above to obtain the quantile effects as $$ \mathsf{QE}_{t}(\tau) := \phi(G_{t}, \tau) -\phi(F_t, \tau), \quad \tau \in [0,1]. $$ The quantile effect measures the contemporaneous impact of the hypothetical policies at different parts of the outcome distribution, and are based on comparisons between counterfactual and actual marginal distributions.
Assume that $y_{it}$ is discrete with a finite support and the process $\{y_{i1}, \ldots, y_{iT}\}$ is ergodic for each $i$. Then the distribution of $y_{it}$ conditional on $\boldsymbol{x}_{it}$ can be represented by a time-homogeneous finite-state Markov chain for each unit. The cross-sectional stationary distribution can then be characterized by aggregating the transition matrices of all the units.
Let the discrete support of $y_{it}$ be $\mathcal{Y}_i = \{y_i^1 < \cdots < y_i^K\}$, which might be different for each unit, and $\boldsymbol{x}_{it}$ include only the first lag of $y_{it}$, i.e. $\boldsymbol{x}_{it} =(1, y_{i(t-1)},\boldsymbol{v}_{it}')'$. For each $i$, let $\boldsymbol{P}_i$ be the $K \times K$ transition matrix with $\boldsymbol{v}_{it}$ fixed at a value $\boldsymbol{v}_i$, which might be different for each unit. The typical element of this matrix can be expressed as
where $\boldsymbol{x}_i^k = (1,y_i^k,\boldsymbol{v}_i)'$. By standard theory for Markov Chains, see, e.g., hamilton2020time, the ergodic probabilities $\boldsymbol{\pi}_i = (\pi_{i1}, \ldots, \pi_{iK})$ are $$ \boldsymbol{\pi}_i = (\boldsymbol{A}_i'\boldsymbol{A}_i)^{-1} \boldsymbol{A}_i'\boldsymbol{e}_{K+1}, \quad \boldsymbol{A}_i = \left(
\right), $$ where $\boldsymbol{I}_K$ is the identity matrix of size $K$, $\boldsymbol{1}_K$ is a $K$-vector of ones, and $\boldsymbol{e}_{K+1}$ is the $(K+1)$th column of $\boldsymbol{I}_{K+1}$. The cross-sectional stationary distribution is $$ F_{\infty}(y) = \operatorname{plim}_{N\to \infty} \frac{1}{N} \sum_{i=1}^N F_{i,\infty}(y), \quad F_{i,\infty}(y) = \sum_{k: y_i^k \leq y} \pi_{ik}, $$ where $F_{i,\infty}$ is a step function with steps at the elements of $\mathcal{Y}_i$.\footnote{Note that $F_{i,\infty}$ is measurable because we assume that the cross-sectional population indexed by $i$ is countable. }
Stationary counterfactual distributions can be formed by replacing $\boldsymbol{\beta}_i(y_i^j)$ by $\boldsymbol{\beta}_i^g(y_i^j)$ in (ref). That is
We denote the resulting cross-sectional stationary distribution as $G_{\infty}$. Note that changes in $y_{i(t-1)}$ do not affect the stationary distribution by the ergodicity assumption. The stationary distribution is useful for analyzing dynamics of the distribution of $y_{it}$ in the long run. We employ it below to examine labor income mobility and the persistence of poverty traps.
In practical applications where the actual distribution of $y_{it}$ may be continuous, we discretize the support $\mathcal Y_i$ separately for each unit. For instance, in the income dynamic study, we divide the income into finitely many categories and cluster income data into the nearest category.
Our estimation procedure can be conducted in two stages. In stage 1, we estimate the coefficients by distribution regression applied separately to the time series dimension of each unit (HDR estimator), and debias the resulting estimates to address the incidental parameter problem.\footnote{Unit-by-unit estimation of random coefficient models without bias correction has been previously considered in the literature; see, e.g., hsiao2012 and pesaran2015time.} In stage 2, we estimate functionals via the plug-in method and further debias if the functionals are nonlinear.
We start with the HDR estimator of $\boldsymbol{\beta}_i(y)$. That is $$ \widetilde\boldsymbol{\beta}_i(y)=\arg\max_{\beta \in \mathbb{R}^{d_x}}Q_{y,i}(\beta),\quad y \in \mathcal{Y}_i, \quad i=1,...,N, $$ where $$ Q_{y,i}(\beta)=\sum_{t=1}^T1\{y_{it}\leq y\}\log \Lambda(-\boldsymbol{x}_{it}'\beta)+\sum_{t=1}^T1\{y_{it} > y\} \log [1-\Lambda(-\boldsymbol{x}_{it}'\beta)], $$ and $\mathcal{Y}_i$ is the set of observed values of the outcome for unit $i$, i.e. $\mathcal{Y}_i := \{y_{i1}, \ldots, y_{iT}\}$. If $\Lambda$ is the standard normal or logistic link, these are standard logit or probit estimators. We then obtain $\widetilde\boldsymbol{\beta}_i(y)$ for other values of $y$ noting that $y \mapsto \widetilde\boldsymbol{\beta}_i(y)$ is a vector of step functions with steps at the elements of $\mathcal{Y}_i$.
Two complications arise in this first stage. First, $\widetilde\boldsymbol{\beta}_i(y)$ is well-defined only if $ y \in [\underline{y}_i, \overline{y}_i)$, where $\underline{y}_i = \min_{1 \leq t \leq T} y_{it}$ and $\overline{y}_i = \max_{1 \leq t \leq T} y_{it}$. Let $N_0(y)$ be the number of indexes $i$ for which $y < \underline{y}_i$, $N_1(y)$ be the number of indexes $i$ for which $y \geq \overline{y}_i$, and $N_{01}(y) = N - N_0(y) - N_1(y)$ denote the number of indexes $i$ for which $\widetilde\boldsymbol{\beta}_i(y)$ exists. Without loss of generality we rearrange the index $i$ such that $\widetilde\boldsymbol{\beta}_i(y)$ exists for all $i = 1,\ldots,N_{01}(y)$. We show below how to adjust the plug-in estimators of the functionals to incorporate the units $i > N_{01}(y)$.
Second, the first stage estimator $\widetilde\boldsymbol{\beta}_i(y)$ is $\sqrt{T}$-consistent but possesses a nonlinear bias of order $T^{-1}$. It is necessary to remove the bias from the coefficients in order that it does not affect the subsequent functionals that employ them. We debias $\widetilde\boldsymbol{\beta}_i(y)$ using analytical methods. That is
where $\widehat B_{i,T}(y)$ is a consistent estimator of the bias of $\widetilde\boldsymbol{\beta}_i(y)$. The specific expressions of the bias and its estimator are presented in the Appendix, where we also consider alternative debiasing methods based on the Jackknife dhaene2015split. While our theory applies to both analytical and Jackknife methods, we focus on analytical methods as they have less demanding data requirements and performed better in our numerical simulations. By removing the $T^{-1}$-bias we can deal with the setting in which $T/N\to\rho\in[0,\infty]$. More precisely, we consider asymptotic sequences where $N=o(T^2)$.
We provide estimators for all the functionals of interest.
\paragraph{Projections of coefficients} A plug-in estimator of $\boldsymbol{\theta}(y)$ corresponds to applying two-stage least squares to (ref) replacing $\boldsymbol{\beta}_i(y)$ by $\widehat\boldsymbol{\beta}_i(y)$. This yields,
where $$ \widehat{\boldsymbol{z}}_i(y) := \sum_{j=1}^{N_{01}(y)} z_j \boldsymbol{w}_{j}' \left(\sum_{j=1}^{N_{01}(y)} \boldsymbol{w}_{j} \boldsymbol{w}_{j}' \right)^{-1} \boldsymbol{w}_{i}. $$ When $\boldsymbol{w}_{i} = \boldsymbol{z}_{i}$, the estimator simplifies to the OLS estimator with $\widehat{\boldsymbol{z}}_i(y) = \boldsymbol{z}_{i}$.
\paragraph{Actual and counterfactual distributions} The plug-in estimators of the actual and counterfactual distributions are
where
Here $\widehat B(y)$ and $\widehat B_G(y)$ are estimators of the first-order bias coming from the nonlinearity of $F_{t}$ and $G_t$ as a functional of $\boldsymbol{\beta}(y)$, $\widehat \Sigma_i(y)^{-1}$ is an estimator of the asymptotic variance matrix of $\sqrt{T}(\widetilde\boldsymbol{\beta}_i(y) - \boldsymbol{\beta}_i(y))$, $\mathsf{tr}$ is the trace operator, and $\ddot\Lambda$ is the second derivative of $\Lambda$. For units for which $\widehat\boldsymbol{\beta}_i(y)$ is not well-defined we set $\Lambda(-\boldsymbol{x}_{it}'\boldsymbol{\beta}_i(y)) = \Lambda(-h(\boldsymbol{x}_{it})'\boldsymbol{\beta}^g_i(y)) = 1$ if $y < \underline{y}_i$ and $\Lambda(-\boldsymbol{x}_{it}'\boldsymbol{\beta}_i(y)) = \Lambda(-h(\boldsymbol{x}_{it})'\boldsymbol{\beta}^g_i(y)) = 0$ if $y \geq \overline{y}_i$.
\paragraph{Quantile effects}
The estimators of the quantile effects are:
where $\widetilde\phi$ is the generalized inverse or rearrangement operator $$ \widetilde\phi(F,\tau) = \int_0^{\infty} 1\{F(y) \leq \tau\} \mathrm{d} y - \int_{-\infty}^0 1\{F(y) \geq \tau\} \mathrm{d} y $$ which monotonizes $y \mapsto F(y)$ before applying the inverse operator.
\paragraph{Stationary distributions}
We start with the empirical transition matrix as a preliminary plug-in estimator of $\boldsymbol{P}_i$, which we modify to enforce that all entries are non-negative and the rows add to one. More precisely, we define the $K \times K$ matrix $\widehat{\boldsymbol{Q}}_i$ with typical element
For each row of $\widehat{\boldsymbol{Q}}_i$, we sort (rearrange) the elements in increasing order to form the matrix $\check{\boldsymbol{Q}}_i$ with typical element $\check Q_{i,jk}$. We then construct the empirical transition matrix $\widehat{\boldsymbol{P}}_i$ with typical element $$ \widehat P_{i,jk} = \check Q_{i,jk} - 1(j>1) \check Q_{i,(j-1)k}. $$ The empirical ergodic probabilities $\widehat{\boldsymbol{\pi}}_i = (\widehat \pi_{i1}, \ldots, \widehat \pi_{iK})$ are now $$ \widehat{\boldsymbol{\pi}}_i = (\widehat{\boldsymbol{A}}_i'\widehat{\boldsymbol{A}}_i)^{-1} \widehat{\boldsymbol{A}}_i'\boldsymbol{e}_{K+1}- \frac{1}{T} \widehat B_{\boldsymbol{\pi}_i}, \quad \widehat{\boldsymbol{A}}_i = \left(
\right), $$ where $\widehat B_{\boldsymbol{\pi}_i}$ is an estimator of the bias coming from the nonlinearity of $\boldsymbol{\pi}_i$ as a functional of $(\boldsymbol{\beta}_i(y_i^1), \ldots, \boldsymbol{\beta}_i(y_i^K))$. We give the expression of $\widehat B_{\boldsymbol{\pi}_i}$ in the Appendix.
The estimator of the stationary distribution is $$ \widehat F_{\infty}(y) = \frac{1}{N} \sum_{i=1}^{N} \widehat F_{i,\infty}(y), \quad \widehat F_{i,\infty}(y) = \sum_{k: y_i^k \leq y} \widehat \pi_{ik}. $$ Estimators of stationary counterfactual distributions can be formed by replacing $\widehat \boldsymbol{\beta}_i(y_i^j)$ by $\widehat \boldsymbol{\beta}_i^g(y_i^j)$ and modifying the bias estimator, $\widehat B_{\boldsymbol{\pi}_i}$, in (ref). The modified expression of the bias estimator is given in the Appendix. The resulting estimator of $G_{\infty}$ is denoted by $\widehat G_{\infty}$.
Inference under flexible heterogeneity, as measured by the variance of the coefficients, comes with a novel uniformity challenge associated with the incidental parameter problem. We now show that the rate of convergence becomes unknown and variable, affecting the asymptotic distribution.
We illustrate the problem through a simple linear model with scalar coefficient to highlight that the standard analytical plug-in methods for inference are not uniformly valid with respect to the degree of heterogeneity.
Consider the model $$ y_{it} = \beta_i+ e_{it}, \quad \mathbb E(e_{it}\mid \beta_i)=0, \quad \mathbb E(\beta_{i}) = \theta, $$ where we allow $\mathsf{Var}(\beta_i)\in[0,C]$, $C>0$, with zero as an admissible value. This class of data generating processes captures different degrees of heterogeneity that might arise in empirical applications. For simplicity, we assume $e_{it}$ and $\beta_i$ are both i.i.d. sequences in both $i$ and $t$ and mutually independent. The estimator of $\theta$ is $$ \widehat\theta= \frac{1}{N}\sum_{i=1}^N\widehat\beta_i,\quad \widehat\beta_i=\frac{1}{T}\sum_{t=1}^Ty_{it} = \beta_i + \frac{1}{T}\sum_{t=1}^T e_{it}. $$ The goal is inference about $\theta$ based on $\widehat\theta$ that remains uniformly valid over $\mathsf{Var}(\beta_i)\in[0,C]$.
Let $\overline{\beta} = \sum_{i=1}^N \beta_i/N$. The asymptotic distribution of $\widehat\theta$ is determined by two components: $$ \widehat\theta-\theta= (\widehat \theta - \overline{\beta}) + (\overline{\beta} - \theta), $$ where
While both terms admit central limit theorems, they may have different rates of convergence. The rate of convergence of $\overline{\beta} - \theta$ depends on the degree of heterogeneity, $\mathsf{Var}(\beta_i)$, which is unknown and this gives rise to a new issue related to estimating the incidental parameters $\beta_1,...,\beta_N$. To illustrate this, consider three special cases:
Any degree of heterogeneity between the above three cases would lead to an unknown rate of convergence $\widehat\theta-\theta=O_P(\xi_{NT})$ where $\xi_{NT}\in[(NT)^{-1/2}, N^{-1/2}]$. Moreover, this unknown rate of convergence has consequences for the properties of standard inferential methods. Note that
A common method to estimate this variance is to plug in sample analogs of $\mathsf{Var}(e_{it})$ and $\mathsf{Var}(\beta_{i})$, $$ \widehat{\mathsf{Var}}(\widehat\theta)=\frac{1}{TN} \widehat{\mathsf{Var}}(e_{it}) + \frac{1}{N}\widehat{\mathsf{Var}}(\beta_i), $$ where $$ \widehat{\mathsf{Var}}(e_{it}) = \frac{1}{TN} \sum_{i=1}^N \sum_{t=1}^T (y_{it} - \widehat\beta_i)^2 = \mathsf{Var}(e_{it}) + o_P(1), $$ and $$ \widehat{\mathsf{Var}}(\beta_i) = \frac{1}{N} \sum_{i=1}^N (\widehat\beta_i - \widehat \theta)^2 = \mathsf{Var}(\beta_i) + O_P\left( \frac{1}{T} \vee \sqrt{\frac{\mathsf{Var}(\beta_i)}{T}}\right). $$
In the moderate and local heterogeneity cases $$ \mathsf{Var}(\widehat \theta) = O\left( \frac{1}{NT}\right), \quad \widehat\mathsf{Var}(\widehat \theta)-\mathsf{Var}(\widehat \theta) = O_P\left( \frac{1}{NT}\right). $$ This leads to incorrect inference of the standard confidence intervals $$ \mathsf{CI}_{1-p}(\theta) = \widehat \theta \pm \Phi^{-1}(1-p/2) \sqrt{\widehat\mathsf{Var}(\widehat \theta)} = \widehat \theta \pm \Phi^{-1}(1-p/2) \sqrt{\mathsf{Var}(\widehat \theta) + O_P((NT)^{-1})}, $$ because $\mathsf{CI}_{1-p}(\theta)$ is scaled by a quantity of the same order as the length of the interval leading to asymptotic distortion $$ \Pr(\theta \in \mathsf{CI}_{1-p}(\theta)) = 1 - p + O(1). $$ The source of the problem is that the estimation error of $\widehat{\mathsf{Var}}(\beta_i)$ does not adapt to degree of heterogeneity because $\widehat{\mathsf{Var}}(\beta_i) - \mathsf{Var}(\beta_i) = O_P(T^{-1})$ when $\mathsf{Var}(\beta_i) = o(T^{-1})$. Note that the solution of setting $\widehat \mathsf{Var}(\beta_i) = 0$ leads to asymptotic under-coverage in the strong heterogeneity case. In the next section, we propose a bootstrap method that is robust to the degree of heterogeneity and is convenient for simultaneous inference on function-valued parameters.
We now develop a simple cross-sectional bootstrap scheme that is uniformly valid over a large class of data generating processes that include both local and strong heterogeneity. We introduce the method in the context of the example from the previous section and provide implementation algorithms for the functionals of interest in our model in Appendix (ref). The formal theoretical results on the validity of cross-sectional bootstrap are given in Theorem (ref).
The cross-sectional bootstrap is based on resampling with replacement of the estimated coefficients $\widehat \beta_i$ instead of the observations $y_{it}$. We call this a cross-sectional bootstrap because it is equivalent to resampling the entire time series $\{y_{i1}, \ldots, y_{iT}\}$ of each cross-sectional unit. Let $\{\widehat\beta_i^*: i=1,..., N\}$ be random sample with replacement from $\{\widehat\beta_i: i=1,..., N\}$. The bootstrap draw of $\widehat \theta$ is $$ \widehat\theta^* =\frac{1}{N}\sum_{i=1}^N\widehat\beta_i^*. $$ We approximate the asymptotic distribution of $\widehat \theta - \theta$ by the bootstrap distribution of $\widehat\theta^* - \widehat\theta$.
Figure (ref) provides a numerical comparison of analytical and cross-sectional bootstrap estimators of the standard deviation of $\widehat \theta$ using a design where $e_{it} \sim \mathcal{N}(0,1)$, $\beta_i \sim \mathcal{N}(\theta, \mathsf{Var}(\beta_i))$, $\mathsf{Var}(\beta_i) \in \{0, 0.1, \ldots, 1\}$, $\theta=1$, $N=100$, and $T=10$. It reports the (true) standard deviation of $\widehat \theta$, based on $\mathsf{Var}(\widehat \theta) = \mathsf{Var}(e_{it})/(NT) + \mathsf{Var}(\beta_i)/N$, as a function of $\mathsf{Var}(\beta_i)$; together with averages over $5,000$ simulations of the following estimators: (1) Standard plug-in: based on $$ \widehat{\mathsf{Var}}(\widehat \theta) = \frac{1}{N^2T^2} \sum_{i=1}^N \sum_{t=1}^T (y_{it} - \widehat\beta_i)^2 + \frac{1}{N^2} \sum_{i=1}^N (\widehat\beta_i - \widehat \theta)^2. $$ This estimator is labeled as “over". (2) Plug-in that omits the heterogeneity in $\beta_i$, based on the first term of the previous expression. This estimator is labeled as “under". (3) Cross-sectional bootstrap interquartile range rescaled by the iterquartile range of the standard normal based on $1,000$ draws.
We find that the standard analytical plug-in estimator overestimates the standard error for any degree of heterogeneity, whereas the analytical plug-in estimator that omits the heterogeneity in $\beta_i$ underestimates the standard error in the presence of any heterogeneity. As predicted by the asymptotic theory, the mean of cross-sectional bootstrap estimator is very close to the standard error uniformly for all the degrees of heterogeneity considered.
The bootstrap algorithms for the model functionals presented in Appendix (ref) are designed to construct confidence bands that cover the functionals simultaneously over the region of points of interest. For example, if we are interested in the scalar function $y \mapsto \xi(y)$ over $y \in \mathcal{Y}$, the asymptotic $p$-confidence band $\mathsf{CI}_p(\xi(y)) := [\widehat \xi_l(y), \widehat \xi_u(y)]$ is defined by the data dependent end-point functions $y \mapsto \widehat \xi_l(y)$ and $y \mapsto \widehat \xi_u(y)$ that satisfy $$ \Pr\left(\widehat \xi_l(y) \leq \xi(y) \leq \widehat \xi_u(y), y \in \mathcal{Y}\right) \to p \text{ as } N,T \to \infty. $$ We illustrate in Section (ref) how these confidence bands can be used to test multiple hypotheses about the sign and shape of the functionals. Pointwise confidence intervals are special cases obtained by setting the region $\mathcal{Y}$ to include only one point.
We employ data from the Panel Study of Income Dynamics for the years 1967 to 1996 survey2020panel. The sample selection follows hu2019semiparametric which restricts the sample to male heads of household working a minimum of 40 weeks.\footnote{This sample is commonly employed in this literature as it represents full time full year workers.} We drop the worker-year observations where labor income is above the 99th sample percentile or below the 1st sample percentile, and keep workers observed for a minimum of 15 years. This selection results in an unbalanced panel with 1,629 workers and 33,338 worker-year observations.
The variables used in the analysis include measures of labor income, years of schooling, number of children, marital status, year of birth, survey year and an indicator denoting the individual is white. The years of schooling variable is constructed from the categorical variable highest grade completed with the following equivalence: 0-5 grades = 5 years, 6-8 grades = 7 years, 9-11 grades = 10 years, 12 grades = 12 years, some college = 14 years, and college degree = 16 years. Following the literature on labor income processes, we construct the outcome, $y_{it}$, as the residuals of the pooled regression of the logarithm of annual real labor income in 1996 US dollars, deflated by the CPI-U-RS price deflator, on indicators for marital status, number of children, year of birth and survey year. We refer to these residuals as labor income.
Unlike most of the previous work in the literature, our model does not explicitly decompose labor income into permanent and transitory components. If we are interested in modeling permanent income and the transitory component accounts for measurement error, working with raw labor income might render the our estimators inconsistent as our model and estimators are nonlinear. A possible solution to this problem is to first separate the permanent component from the transitory component using deconvolution methods and then work with the permanent component. lee2025identificationestimationdynamicrandom adopts this approach and finds similar results for persistence with raw labor income and the extracted permanent component using a dynamic random coefficient model applied to data from the PSID. Based on this evidence, we do not pursue this approach and work with raw labor income.
We estimate the HDR model (ref) with $\boldsymbol{x}_{it} = (1, y_{i,t-1},\boldsymbol{v}_{it})'$, where $\boldsymbol{v}_{it}$ is the age of the worker $i$ at time $t$. We denote the model coefficients by $\boldsymbol{\beta}_i(y) = (\alpha_i(y), \rho_i(y), \pi_i(y))'$, where we refer to $y \mapsto \alpha_i(y)$ as the intercept or level function and $y \mapsto \rho_i(y)$ as the slope or persistence function, and their bias corrected estimates by $\widehat \boldsymbol{\beta}_i(y) = (\widehat \alpha_i(y), \widehat \rho_i(y), \widehat \pi_i(y))'$. These estimates are obtained using (ref). The left panel of Figure (ref) (Between Median) plots the kernel density of the estimated slope function $\widehat\rho_i(y)$ at a fixed value of $y$ corresponding to the sample median of $y_{it}$ pooled across workers and years. We find substantial heterogeneity between workers in this parameter. The density of the persistence coefficient includes both positive and negative values corresponding to positive and negative state dependencies in the labor income process at the median. The right panel of Figure (ref) (Within Median) plots the pointwise sample median of the function $y \mapsto \widehat \rho_i(y)$ over a region $\mathcal{Y}$ that includes all the sample percentiles of the sample values of $y_{it}$ pooled across workers and years. The function is plotted with respect to the probability level of the sample percentile to facilitate interpretation. We find substantial heterogeneity in the slope within the distribution of the median worker. The slope is increasing with the percentile level indicating higher persistence parameter at the upper tail of the distribution. The two figures combined illustrate the existence of substantial heterogeneity in income dynamics both between and within workers.
We compare the estimates from our model with other alternatives in terms of measures of income persistence. In particular, we compare estimates of the quantile derivative with respect to the lagged dependent variable in (ref) obtained from our model with estimates obtained from the quantile-based model of arellano2017earnings.
For our HDR model, we estimate the derivative function $\dot\boldsymbol{\beta}_i(y)$ by linear interpolation from the grid $\{\dot{\widehat{\boldsymbol{\beta}}}_i(y_i^{(m)}), m=1,...,T\}$, where $ y_i^{(1)}\leq ...\leq y_i^{(T)} $ is the sorted sample of $\mathcal{Y}_i$ and $$ \dot{\widehat{\boldsymbol{\beta}}}_i(y_i^{(m)})=
$$ The derivative in \eqref{eq:qder} needs to be evaluated at $y= Q_{y_{it} }(\tau \mid \boldsymbol{x}_{it})$. Since the HDR model is not designed to estimate conditional quantiles directly, we estimate $ Q_{y_{it} }(\tau \mid \boldsymbol{x}_{it})$ using the quantile regression model in ((ref)).
For the quantile regression (QR) model, we use the specification
where $P(\boldsymbol{x}_{it}) = (1,y_{i(t-1)},y_{i(t-1)}^2,\boldsymbol{v}_{it})'$, and, unlike arellano2017earnings, we model total labor income rather than the persistent component. The quantile derivative corresponding to this model is $$ \dfrac{\partial Q_{y_{it}}(\tau \mid \boldsymbol{x}_{it})}{\partial y_{i(t-1)}} = \dfrac{ \partial P(\boldsymbol{x}_{it})}{\partial y_{i(t-1)}}'\boldsymbol{\beta}_t(\tau), \quad \dfrac{ \partial P(\boldsymbol{x}_{it})}{\partial y_{i(t-1)}} = (0,1,2y_{i(t-1)},0)', $$ which is nonlinear in $y_{i(t-1)}$ and heterogeneous over $t$, but is homogeneous with respect to $i$ for $y_{i(t-1)}=y$.
Figure (ref) plots estimates of the quantile derivative surface $(y_{i(t-1)},\tau) \mapsto \partial Q_{y_{it} }(\tau\mid \boldsymbol{x}_{it})/\partial y_{i(t-1)}$ evaluated at $\tau=(0.05,0.1,...,0.95)$, and $y_{i(t-1)}$ over its sample percentiles in the cross-sectional distribution at $t-1=1986$. The value of age in $\boldsymbol{x}_{it}$ is set at its observed value in year $t=1987$.
The top panel reports the median of the estimated derivatives obtained from our HDR model. We report the median quantile surface because the quantile derivative is heterogeneous across $i$. The bottom panel is produced by the nonlinear quantile regression (QR) model in (ref) fixing $t=1987$. As a middle-ground comparison, the center panel plots the estimated quantile derivative using the following nonlinear distribution regression (NDR) model with time-varying coefficient:
Without accounting for sample variation, These models produce quite different results. Both DR and QR produce nonlinear derivative maps, but HDR delivers smaller median quantile derivatives. Comparing our model with results produced by the other two models, we conclude that most of the differences can be explained by between-unit heterogeneity, as QR and HDR produce derivatives of similar magnitudes.\footnote{In results not reported, we find that the estimates of NDR do not change much if we replace $P(\boldsymbol{x}_{it})$ by $\boldsymbol{x}_{it}$. In other words, the difference between the upper panel and the other panels is not driven by nonlinearity of the index function $\boldsymbol{x}_{it} \mapsto P(\boldsymbol{x}_{it})'\boldsymbol{\beta}_t(y)$.} This result is consistent with the stylized fact in the literature that heterogeneous income profile models find more persistence than restricted income profile models drewianka2020earnings.\footnote{In results not reported, we find that the HDR model produces similar patterns of conditional skewness to arellano2017earnings. In particular, we find that the conditional income distribution is skewed to the right for individuals at the bottom; and is skewed to the left for individuals at the top of the distribution.}
We explore if specific worker characteristics are associated with the heterogeneity in the level and persistence coefficients using projections. In particular, we apply (ref) with $\boldsymbol{z}_i$ including a constant, the initial labor income, years of schooling, a white indicator and year of birth, and $\boldsymbol{w}_i=\boldsymbol{z}_i$.\footnote{There might be common unobservable variables affecting the coefficients and covariates $\boldsymbol{z}_i$. For example, unobserved ability might affect both the level coefficient and education. Ideally, we would like to implement an IV approach, but it is not feasible due to lack of credible instruments in the dataset.}
Figure (ref) reports the estimates and 90% confidence bands of the projection coefficient function $y \mapsto \boldsymbol{\theta}(y)$ for education over a region $\mathcal{Y}$ that includes all the sample percentiles of the pooled sample of $y_{it}$ with probability levels $\{0.10, 0.11, \ldots, 0.90\}$, plotted with respect to these probability levels. We find the education level is associated with coefficient heterogeneity at some locations of the distribution. For example, the persistence parameter $\rho_i(y)$ is negatively associated with education at the bottom of the distribution, whereas the level parameter $\alpha_i(y)$ is positively associated with education in the middle of the distribution. The effect of education on $\rho_i(y)$ is increasing with $y$, although this pattern should be interpreted carefully as the function is not very precisely estimated, as reflected by the width of the confidence band.
An important implication of the HDR representation of labor income is that an individual's location in the income distribution in a specific time period partially depends on his location in previous periods. Moreover, the nature of this dependence varies by worker. This indicates that a shock to current labor income will determine the path of future income.
To further illustrate the presence and heterogeneity of this dependence we examine the impact on future income resulting from a negative shock to initial income. We implement the shock by reducing labor income in 1985 by 25 percent simultaneously for all individuals.\footnote{We choose 1985 as the base year because it is the year with the largest number of observations in the dataset.} We interpret this as an unanticipated shock in that we change the level of initial income but keep all other aspects of the model constant. Specifically, we estimate the counterfactual distribution (ref) for the transformation $$h_{it}(\boldsymbol{x}_{it}) = (1, y_{i(t-1)}+\log(1-\kappa), \boldsymbol{v}_{it})'$$ with $\kappa = 0.25$. This transformation yields a counterfactual distribution of labor income in $t=1986$. We also estimate the actual distribution and the corresponding quantile effects. We compare our estimates with those from the following alternative models:
(a) Hetero.DR: the proposed model: $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it}) = \Lambda (\boldsymbol{\beta}_i(y)'\boldsymbol{x}_{it} )$.
(b) Homo.DR: the homogeneous DR: $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it}) = \Lambda (\boldsymbol{\beta}(y)'\boldsymbol{x}_{it} )$.
(c) Hetero.AR: the heterogeneous AR model: $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it}) = \Lambda ((y-\boldsymbol{\beta}_i'\boldsymbol{x}_{it})/\sigma ).$
(d) Homo.AR: the homogeneous AR model: $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it}) = \Lambda ((y-\boldsymbol{\beta}'\boldsymbol{x}_{it})/\sigma ).$
(e) AR-fixed effect: the AR model with fixed effects: $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it}) = \Lambda ((y-\boldsymbol{\beta}'\boldsymbol{x}_{it} - \alpha_i)/\sigma ).$
(f) QR: the quantile regression model in (ref).
The parameters of the AR models are estimated by least squares, the parameters of the DR models are estimated with $\Lambda$ equal to the standard logistic distribution, and the parameters of the QR model are estimated by quantile regression. We estimate the cross-sectional CDF of $y_{i,t}$ for QR as $$ \frac{1}{N}\sum_{i=1}^N\int 1\{ \widehat \boldsymbol{\beta}_{t}(\tau)'P(\boldsymbol{x}_{it})\leq y\} \mathrm{d} \tau $$ where $\widehat \boldsymbol{\beta}_{t}(\tau)$ is the QR estimate of $\boldsymbol{\beta}_{t}(\tau)$.
Figure (ref) reports estimates and 90% confidence bands of the quantile effects from Hetero.DR, together with the estimates obtained from the alternative models. The confidence bands are computed using Algorithm (ref) with $p=.90$, $B=500$ and $\mathcal{T} = \{.05, .06, \ldots, .95\}$. The estimates show that the fully homogeneous location-shift and DR models predict that the income shock reduces next period income almost on a one-for-one basis throughout the distribution. The linear AR models with fixed effects lower the effect to about 15% and 10%, whereas the HDR model further ameliorates it to about 5%. Quantile regression (QR) provides estimates similar to the fully homogeneous models. The confidence bands of the HDR model indicate there is no evidence of heterogeneous effects across the distribution. Moreover, they do not fully cover the estimates of the other three models. In results not reported, we find that joint confidence bands from these models do not fully overlap with the confidence bands of the HDR model.\footnote{The confidence level of the joint bands is corrected by the union bound to $97.5\%$, in order to preserve the joint coverage to at least 90%.} We can formally reject the homogeneity restrictions imposed by the alternative models. From the comparison of the difference estimators, we conclude that while between-unit heterogeneity matters the most, nonlinearity also plays a significant role in assessing the impact of the negative shock. Consistent with the results for the quantile derivative in Figure (ref), QR estimates less persistence than HDR.
We now analyze labor income mobility and the existence of “relative poverty" traps. We evaluate the probability of remaining in lower locations of the residual distribution noting that we refer to this as relative poverty as we acknowledge that the total income level may not be below the poverty line. We do so via the model from Section (ref), where the conditional distribution is represented by a Markov chain. We treat income as discrete and set the states for each worker as the observed values of $y_{it}$, that is $\mathcal{Y}_i:=\{y_{it}: t=1,...,T\}$ and $K=T$, and set $\boldsymbol{v}_{it}$ to the median value of age in the sample ($\boldsymbol{v}_i = 36$ for all $i$).
Following hu2019semiparametric, consider the following probabilities to describe mobility $$ P_i(p,q,h):= \Pr(y_{i(t+h)}< y_p \mid y_{it}< y_q, \mathcal{F}_{it} ),\quad i=1,...,N, $$ where $y_p$ and $y_q$ are the $p$-quantile and $q$-quantile of the distribution of labor income. These probabilities correspond to the following experiment: If we exogenously set labor income below $y_q$ at time $t$, then $P_i(p,q,h)$ is the probability labor income is below $y_p$ after $h$ years.\footnote{The probability $P_i(p,q,h)$ is identified if $y_{it}$ is observed below $y_p$ for some $t$. We restrict the sample to workers that satisfy this condition in the sample period to estimate these probabilities.} For example, if we define the poverty line as the $10$-percentile, then $P_i(0.1,0.1,5)$ is the probability that worker $i$ would remain in poverty after 5 years if he falls below the poverty line due to, for example, a negative income shock.
Our model allows the probabilities $P_i(p,q,h)$ to be heterogeneous across workers. To summarize this heterogeneity, we can examine the average probability
For instance, $\bar P(0.3, 0.1, 1)$ is the probability that a randomly chosen worker is below the 30-percentile if in the previous year he was below the 10-percentile. We also examine quantiles of the probabilities such as $$ Q_{\tau}(p, q, h) $$ which denotes the $\tau$-quantile of $\{P_i(p, q,h):, i=1,...,N\}$ for fixed $(p, q, h)$. For example, $Q_{0.25}(0.3, 0.1, 1)$ is the first quartile of the probability that a worker is below the 30-percentile if in the previous year he was below the 10-percentile.
The upper panel of Figure (ref) plots $p \mapsto \bar P(p,q,h)$ for $p \in [0,0.5]$, $q \in \{0.1,0.25,0.5\}$ and $h \in \{1,2,5\}$. We find heterogeneity with respect to the initial condition that vanishes with time due to the ergodicity of the process. The probability that a randomly selected worker remains below the 10-percentile after one year is more than 50%, whereas this probability decreases by about half if the worker was initially below the median. This difference in probabilities reduces after two years and almost vanishes after five years. The lower panel of Figure (ref) plots $p \mapsto Q_{\tau}(p,q,h)$ for $p \in [0,0.5]$, $q =0.1$, $h \in \{1,2,5\}$ and $\tau \in \{0.1,0.5,0.9\}$. We uncover significant heterogeneity across workers that is hidden in the analysis of the mean worker. Even after 5 periods the deciles of the probability of remaining below the 10-percentile range from $0$ to over $0.9$. This illustrates the importance of accounting for heterogeneity in understanding the probability of escaping poverty.
Let $h_i(p)$ denote the recurrence time of $y_p$. That is, starting from $\{y_{it} < y_p\}$, the number of years $h$ until the first occurrence of $\{y_{i(t+h)} > y_p\}$. For example, if $y_{0.10}$ is the poverty line, $h_i(0.10)$ is a random variable that measures the number of years that worker $i$ takes to escape from poverty. Then, $$ \Pr(h_i(p) = h) = \Pr(y_{i(t+h)}>y_p, y_{i(t+h-1)}<y_p,...,y_{i(t+1)}<y_p \mid y_{it}<y_p, \mathcal{F}_{it}), $$ which can be expressed as a functional of the parameters of the HDR model. Another interesting quantity is $$H_i(p) =\sum_{h} h \Pr(h_i(p) = h),$$ which gives the expected recurrence time for each individual. In the previous example, $H_i(0.10)$ gives the expected number of years that worker $i$ would take to escape from poverty. Figure (ref) plots a histogram of the estimated $H_i(0.10)$. More than 60% of the workers would escape from the poverty in two or less years, but about 10% of the workers would stay for more than 20 years. Table (ref) reports several quantiles of the estimated $H_i(0.1)$ for groups stratified by education and race. We find substantial heterogeneity between workers associated with education and race. Whereas the deciles of the expected recurrence time range from 1 to 7 years for workers with at least high school, the corresponding value of 176 years indicates there are more than 10% of workers with less than high school that would never escape poverty. The distribution of the expected recurrence time also differs by race. The upper decile of the expected recurrence time is about 20 years higher for nonwhite than for white workers. This heterogeneity in the persistence of poverty has clear implications for the design of poverty alleviation policies. As they employ a different sample to ours and employ a different definition of “relative poverty" we do not directly compare these results to lillard1978dynamic. However, in addition to confirming the dependence in labor income documented in their study, we illustrate the remarkable difficulty facing some workers in escaping relative poverty.
Finally, we examine the goodness of fit of our HDR model to the PSID data Figure (ref) compares the empirical distributions of $y_{it}$ in 1981 and 1991 with the corresponding distributions predicted by the HDR model. The model provides a remarkably close fit to the empirical distribution for all the values of $y$, including the tails. In results not reported, we find that the HDR model also provides a good fit of income dynamics by comparing model-based and empirical estimates of the autocorrelation of income.
This section develops asymptotic theory for the estimators of the functionals of interest. We start by introducing some notation. Recall that the loss function for the estimation of the coefficients is: $Q_{y,i}(b)=T^{-1}\sum_{t=1}^Tq_{y,it}(b)$, where $$ q_{y,it}(b) := 1\{y_{it}\leq y\}\Lambda(-\boldsymbol{x}_{it}'b)+1\{y_{it} > y\}[1-\Lambda(-\boldsymbol{x}_{it}'b)]. $$ Let
where all terms are evaluated at the true value of $\boldsymbol{\beta}_i(y)$. For a generic function $q(\boldsymbol{\beta})$ where $\boldsymbol{\beta}$ is a $d_{\beta}$-dimensional vector and $d_{\beta}:=\text{dim}(\boldsymbol{\beta})$, $\nabla q(\boldsymbol{\beta})$ denotes the gradient $d_{\beta}$-dimensional vector whose components are the partial derivatives of $\boldsymbol{\beta} \mapsto q(\boldsymbol{\beta})$; and $\nabla^2 q(\boldsymbol{\beta})$ is the Hessian matrix. In addition, $\nabla^3 q(\boldsymbol{\beta})$ is a $d_{\beta}\times d_{\beta}^2$ matrix, defined as $(\nabla B_1(\boldsymbol{\beta}),...,\nabla B_{d_{\beta}}(\boldsymbol{\beta}))$, where $\nabla B_j(\boldsymbol{\beta})$ is the $d_{\beta}\times d_{\beta}$ Jacobian of the $j$ th row of $\nabla^2 q(\boldsymbol{\beta})$.
The following assumptions relate to the properties of the sampling process. Recall that $\mathcal F_{i1}\subset...\subset \mathcal F_{iT}$ is the sequence of filtrations over time, and $ \boldsymbol{x}_{it}$ is updated to $\mathcal F_{it}$. The underlying probability space is equipped with a probability measure $P_T$. Let $\mathcal P$ denote the collection of all DGPs $P_T$ where our model holds
We assume that the following assumptions hold on $P_T$. Throughout these assumptions, we let $c$ and $C$ be absolute constants, which means they do not depend on the specific DGP $P_T\in \mathcal P$. Under different $P_T$, the degree of heterogeneity will vary and lead to different rates of convergence. We aim to establish inference results which are uniformly valid in $\mathcal P$.
We will allow the distribution of $y_{it}$ to be either continuous or discrete. When it is continuous, let $\mathcal Y$ be a compact subset of the support of $y_{it}$ on which the density of $y_{it}$ conditional on $\boldsymbol{x}_{it}$ is bounded away from zero. When $y_{it}$ is discrete, let $\mathcal Y$ be the set of discrete values of the support.
For a given integer $M>0$, let $Y_M=(y_1,...,y_M)'$ be an arbitrary $M$-dimensional vector on $\otimes_{i=1}^M\mathcal Y.$ Let $S_{wz}:= C_{1}^{-1}C_{2}(C_{2}'C_{1}^{-1}C_{2})^{-1}$ where $C_{1}=\mathbb{E} \boldsymbol{w}_{i}\boldsymbol{w}_{i}'$, $C_{2}=\mathbb{E} \boldsymbol{w}_{i}\boldsymbol{z}_{i}'$, and
Consider a covariance kernel is given by the limit of the elements of the following $M\times M$ matrix $$ H_{\eta,NT}=(H_{\eta, NT}(y_k, y_l))_{M\times M} $$ which is an $M\times M$ matrix with the $(k,l)$ element as: $$ H_{\eta,NT}(y_k, y_l)=\frac{ \eta' \Sigma_{NT}(y_k, y_l)\eta}{ [ \eta' \Sigma_{NT}(y_k)\eta]^{1/2} [ \eta' \Sigma_{NT}(y_l)\eta]^{1/2} } $$ and $\eta\in\mathbb R^{\dim(\mathsf{vec}\boldsymbol{\theta})}$. We make the following assumption regarding this covariance kernel:
For a generic estimator $\widehat F(y)$ of $F(y)$, which is either $F_t$ or $G_t$, one can show that it has the following expansion (proved in ((ref))): $$ \widehat F(y)-F(y) =\frac{1}{N }\sum_{i=1}^N[\frac{1}{\sqrt{T}}d_{\psi,i}(y) + d_{\boldsymbol{\gamma},i}(y)]+ o_P(\zeta_{NT}(y)) $$ where $\zeta_{NT}(y)= (NT)^{-1/2}\mathsf{Var}_t(d_{\psi,i}(y))^{1/2}+ N^{-1/2}\mathsf{Var}_t(d_{\boldsymbol{\gamma},i})^{1/2}$, and the two leading terms $d_{\psi,i}(y) $ and $d_{\boldsymbol{\gamma},i}(y)$ are asymptotically independent, and respectively capture the sampling variation from the first-stage and second-stage. The quantile effect has similar expansions
where $\bar \zeta_{NT}(\tau)= (NT)^{-1/2}\mathsf{Var}_t(p_{\psi,i}(\tau))^{1/2}+ N^{-1/2}\mathsf{Var}_t(p_{\boldsymbol{\gamma},i}(\tau))^{1/2}$, and $p_{\psi,i}(\tau) $ and $p_{\boldsymbol{\gamma},i}(\tau)$ are zero-mean uncorrelated terms. The formal definitions of $(d_{\psi,i}, d_{\boldsymbol{\gamma}, i},p_{\psi,i}, p_{\boldsymbol{\gamma}, i})$ depend on the specific $F\in\{F_t, G_{t}\}$, which are given in the Appendix.
The following condition bounds the moments. For notational simplicity, we write $$ V_{\boldsymbol{\gamma}}(y ):= V_{\boldsymbol{\gamma}}(y, y),\quad V_{\psi}(y ):= V_{\psi}(y, y). $$ Recall $\Lambda(s)$ denotes the link function of the distribution regression. Let $\dot{\Lambda}(s)= \mathrm{d}\Lambda(s)/\mathrm{d} s$ and $\ddot{\Lambda}(s)=\mathrm{d}^2\Lambda(s)/\mathrm{d} s^2$.
Assumption (ref) requires that the data are cross-sectionally independent. Assumption (ref) imposes conditions regarding serial dependence. We impose two high level conditions regarding the empirical process for weakly dependent data. It requires some primitive conditions, e.g., mixing conditions, so that $\{(Y_{it}, \boldsymbol{x}_{it}): t\leq T\}$ is serially weakly dependent. Additionally, we do not assume stationarity when the analytical debias is used to address the incidental parameter problem.
Assumption (ref) is used to establish the finite dimensional distribution (f.i.d.i.) of $ \eta'\mathsf{vec}(\widehat\boldsymbol{\theta}(\cdot)- \boldsymbol{\theta}(\cdot)) $, which is required for a given $Y_M, M$ and $\eta$. Therefore, the constant $c_{Y_M,\eta}$ is allowed to depend on these parameters. To show that Assumption (ref) is reasonable even though the variance of $\boldsymbol{\gamma}_i(y)= \boldsymbol{\beta}_i(y) - \boldsymbol{\theta}(y)\boldsymbol{z}_{i}$ may vary across $y$ in the second-stage regression, we consider the following model
Here $\xi_{NT}(y)$ is a bounded non-stochastic sequence that may converge to zero, whose rate depends on $y$; $\bar \boldsymbol{\gamma}_i(y)$ is a random vector of “normalized" $\boldsymbol{\gamma}_i(y)$ , so $V_{\bar \boldsymbol{\gamma}}(y,y)$ can be understood as a normalized covariance matrix. Hence the strength of $\boldsymbol{\gamma}_i(y)$ is determined by the rate of convergence of $\xi_{NT}(y)$. Given this setting, consider the following special cases:
Thus each element has a limit given on the right hand side. With sufficient variation across $y_k$, the limit of the matrix $H_{\eta, NT}$ is non-degenerate and satisfies ((ref)).
Assumption (ref) (i) requires that the fourth moments of $\boldsymbol{\gamma}_i(y)$ and $d_{\boldsymbol{\gamma}, i}(y)$ are bounded by their second moment up to a constant, uniformly in $y$. To see the plausibility of this condition, again consider model ((ref)). Then the left hand side of condition (i) becomes $$ \mathbb{E}\left[ \sup_{y\in\mathcal Y}\left(\frac{\| \boldsymbol{\gamma}_i(y)\boldsymbol{w}_{i}' \| ^2 }{ \lambda _{\min}(V_{\boldsymbol{\gamma}}(y))}\right)^{2} \right] = \frac{ \frac{1}{N}\sum_{i=1}^N \mathbb{E}(\| \bar\boldsymbol{\gamma}_i(y)\boldsymbol{w}_{i}' \| ^{4}) }{ \inf_{y\in\mathcal Y}\lambda _{\min}^{2}(V_{\bar \boldsymbol{\gamma}}(y,y))} , $$ which is upper bounded by a constant provided $\frac{1}{N}\sum_{i=1}^N\mathbb{E}(\| \bar\boldsymbol{\gamma}_i(y)\boldsymbol{w}_{i}' \| ^{4})<C.$ Other conditions of this assumption are standard. Condition (ii) requires higher moments to be bounded. For instance, we need $\mathbb{E} \|\boldsymbol{w}_{i}\|^{8+c}<C$ where $\boldsymbol{w}_i$ is a vector of characteristics such as initial labor income, years of schooling and race. These are standardized so it is plausible to assume they have high moments.
Conditions (iii) and (iv) identify the parameters $\boldsymbol{\theta}(y)$ and $\boldsymbol{\beta}_i(y)$. To see this, note that the model implies $$ -\frac{1}{T}\sum_{t=1}^T\mathbb E_i\boldsymbol{x}_{it}\Lambda^{-1}\left( \Pr(y_{it}\leq y|\mathcal F_{it})\right)= \left(\frac{1}{T}\sum_{t=1}^T\mathbb E_i\boldsymbol{x}_{it}\boldsymbol{x}_{it}'\right)\boldsymbol{\beta}_i(y). $$ Inverting $\frac{1}{T}\sum_{t=1}^T\mathbb E_i\boldsymbol{x}_{it}\boldsymbol{x}_{it}'$ leads to the identification of $\boldsymbol{\beta}_i(y)$. In addition, $\mathsf{rank}(C_2)\geq \dim(\boldsymbol{z}_{i})$ implies the identification of $\boldsymbol{\theta}(y).$
When the support of $y_{it}$ is continuous, Assumption (ref) imposes continuity of moment, the link and the $\boldsymbol{\beta}_i(y)$ functions. In particular, as the condition is imposed on rescaled functions in Condition (ii), so that it is not affected by the unknown strength of $\boldsymbol{\gamma}_i(y)$ and $(d_{\boldsymbol{\gamma},i}(y),p_{\boldsymbol{\gamma},i}(y))$.
In the next theorem, $L$ denotes the number of lags used for the Newey-West truncation for long-run variance, which is needed for analytical bias corrections.
The theorems below additionally require Assumption (ref), which are based on some additional notation for the stationary distribution. We present them in the appendix.
Theorem (ref) refers to the estimated distributions. It allows the support of $y_{i,t}$ to be either continuous or discrete with finitely-many states such that the stationary distribution can be modeled using Markov chains with finite states. For the processes $F_n(\cdot)$ and $\mathbb G(\cdot)$, we define $F_n(\cdot)\Rightarrow \mathbb G(\cdot)$ in $\ell^{\infty}(\mathcal Y)$ as the weak convergence in the set of bounded functions on $\mathcal Y$.
Theorem (ref) refers to the estimated quantile effect when we further require the conditional distribution of $y_{it}$ given $\boldsymbol{x}_{it}$ is continuously differentiable. We do not consider the quantile effect of the stationary distribution. Let $\mathcal T\subset (0,1)$ be a set of quantile indices, such that $\{\phi(F_t,\tau): \tau\in\mathcal T\}\subset \mathcal Y$.
Theorem (ref) shows the uniform validity of cross-sectional bootstrap over a large class of data generating processes with varying degrees of coefficient heterogeneity.
We now provide some simulation evidence documenting the finite-sample performance of our method. The online appendix includes additional simulation results.
Consider the dynamic distribution regression model:
with $$\theta(y)= 3 \text{ sgn}(y-2)(y-2)^2, \ y\in \mathcal Y.$$ We set $\mathcal Y= \{1.7, 1.8,..., 2.3\}$, where the two endpoints of $\mathcal Y$ are chosen to avoid the estimation of extreme quantiles. The marginal probabilities $\Pr(y_{it}<1.7)$ and $\Pr(y_{it}>2.3)$ are both approximately 0.1. We generate the simulated data by independently drawing $(e_{it}, w_{i}, \bar\gamma_i)$ from: $$ e_{it}\sim\mathcal N(0,1),\quad w_{i} \sim \text{Uniform}(1.5,2.5), \quad \bar\gamma_i\sim \text{Uniform}(-0.5,0.5).$$ Finally, $y_{it}$ is initialized by $ y_{i0}\sim \text{Uniform}(0.52,1.52)$, and iteratively generated via $$ y_{it}=\theta^{-1}\left( \frac{e_{it}}{y_{i(t-1)}(w_{i}+\bar\gamma_i)}\right). $$ The parameters of this DGP are chosen so that $y_{i(t-1)}(w_{i}+\bar\gamma_i)>0$ for all $t$ a.s. Therefore, $ \Pr(y_{it}\leq y \mid \mathcal{F}_{it})=\Phi(y_{i(t-1)}\beta_i(y))$ is satisfied.
The object of interest is $\theta(y)$. Figure (ref) plots the variance of $\gamma_i(y) = \theta(y)\bar\gamma_i$, the noise level of $\beta_i(y)$, across $y\in\mathcal Y$. By construction, $\mathsf{Var}(\gamma_i(y))$ degenerates at $y = 2$, and increases as $y$ deviates from 2, which affects the rate of convergence for estimating $\theta(y).$ The right panel plots the true standard error of the estimator $\widehat\theta(y)$, along with three estimators: the proposed bootstrap interquartile range (IQR) $se^*(y)$ defined as $ se^*(y) = (q^*_{.75}(y)-q^*_{.25}(y))/(z_{.75}-z_{.25})$, where $q_p^*$ is the bootstrap $p$-quantile of $\boldsymbol{\eta}'\mathsf{vec}(\widehat\boldsymbol{\theta}^*_b(y)-\widehat\boldsymbol{\theta}(y))$ and $z_p$ is the $p$-quantile of the standard normal. The IQR is a consistent estimator for the asymptotic standard deviation, which is often used to replace the bootstrap variance, as it is challenging to show the latter is consistent in most cases.
The other two estimators are respectively “Plug-in-over" and “Plug-in-under", respectively defined below. The plug-in methods are clearly not robust to changes in $\mathsf{Var}(\gamma_i(y))$ across $y$.
We examine the coverage properties of $\theta(y)$ and compare four inferential methods:
(i) Proposed: the proposed uniform inference procedure using the interquartile range described in Remark (ref).
(ii) No-debias: this method does not debias, while all other steps are the same as the proposed method.
(iii) Plug-in-over: this method plugs in the estimated standard error, it uses the estimated $V_\psi(y)$ and $V_\gamma(y)$ by:
where computing the estimators $\widehat {\mathbb A}_{y,1i} $ and $\widehat\gamma_i(y)$ is straightforward. Meanwhile, we apply the Newey-West type estimator $\Xi(y)$ to estimate $\mathbb{E}( \frac{1}{T}\sum_{s,t\leq T} \psi_{it}(y_k) \psi_{it}(y_l)' \mid W) $.
(iv) Plug-in-under: this method also plugs in the estimated standard error, but replaces $\widehat\Sigma_{NT}(y)$ of the Plug-in-over method with $$ \widetilde\Sigma_{NT}(y) =\frac{1}{NT}\widehat V_{\psi}(y). $$
Table (ref) summarizes the coverage probabilities of $ \{\theta(y): y\in {\mathcal{Y}}\} $ out of 1,000 replications. The results are generally as expected. The no-debias method performs unsatisfactorily when $T\leq N $ due to the incidental parameter bias issue. The plugin-over method assumes that there is arbitrary heterogeneity in $\beta_i(y)$, so is quite conservative; the plugin-under method is the standard treatment in the varying-coefficient literature, which assumes that the heterogeneity in $\beta_i(y)$ can be fully captured by covariates $\boldsymbol{w}_i$. The confidence band resulting from $ \widetilde\Sigma_{NT}(y)$ undercovers $\theta(y)$.
We simulate a dynamic distribution regression model from a heterogeneous-coefficient autoregressive model with calibrated parameters using the PSID data. Specifically, we first estimate the following model
for each in-sample individual $i$ to calibrate $\beta_{0i},\beta_{1i}$ and $\sigma$, which we use to calibrate the second-stage model parameters $\left(\theta_{00},\theta_{01},\theta_{10},\theta_{11}\right)$ and $\sigma_{0},\sigma_{1}$ by estimating the regressions
where $y_{it}$ is the outcome variables (residual log income), and $\boldsymbol{w}_{1i}$ is the vector of individual characteristics consisting of the variables initial labor income, years of schooling, a white indicator and year of birth.
Let $\boldsymbol{x}_{it}:=\left(1,y_{i(t-1)}\right)'$, $\boldsymbol{w}_{i}:=\left(1,\boldsymbol{w}_{1i}'\right)'$ and assume $\epsilon_{it}\sim \mathcal N(0,1)$. We then rewrite (ref) and (ref) as following heterogeneous dynamic distribution regression model
where
We then: (1) simulate $ y_{it}$ according to models (ref) and (ref) based on the calibrated values of $\sigma$, $\sigma_{0}$, $\sigma_{1}$, $\theta_{00}$, $\boldsymbol{\theta}_{01}$, $\theta_{10}$ and $\boldsymbol{\theta}_{11}$; (2) calculate the implied distribution regression parameters $\widetilde{\beta}_{0,i}\left(y\right)$, $\widetilde{\beta}_{1,i}\left(y\right)$, $\widetilde{\boldsymbol{\theta}}_{0}\left(y\right)$ and $\widetilde{\boldsymbol{\theta}}_{1}\left(y\right)$ based on (ref); and (3) employ our proposed estimation methods to estimate the quantile treatment effects of a counterfactual increase of the years of schooling variable by 1 year for every unit in the sample. In this exercise, the size of the panel $N$ and each individual's length of observations $T_i$, are the same as the PSID data. That is, $N=1649$, and the average $T_i$ is 20.5.
Table (ref) compares the MSEs of our proposed estimators using analytical debiasing with the estimator which does not debias. We find that the analytical bias correction yields reductions in the mse between $21$ and $45\%$ depending on the quantile index. In results not reported, we find that the Jackknife debiasing does not reduce the mse of the uncorrected estimator in this case.
We develop estimation and inference methods for dynamic distribution regression panel models that incorporate heterogeneity both within and between units. Our model can be employed in a large number of empirical settings. An empirical investigation of labor income processes illustrates some economic insights our approach can provide. We find that accounting for individual heterogeneity is important in studying the potential impact of income shocks on future income and understanding income mobility and poverty persistence.
We could extend our model in several directions. For instance, we could explicitly include time fixed effects and covariates with homogeneous coefficients in the first stage. This could be useful in empirical applications which directly model an outcome variable with trends rather than the residuals. To reduce the number of estimated parameters, we could model the individual coefficients in HDR using factor structures as in chernozhukov2018inference. We could also reduce dimensionality by modeling the between and within heterogeneity though a pseudo-factor structure where the value $y$ plays the role of time. For example, in the empirical application we can model the persistence coefficient as $\rho_i(y) \approx \boldsymbol{\lambda}_i'\boldsymbol{f}_y$, where $\boldsymbol{\lambda}_i$ is a vector of loadings and $\boldsymbol{f}_y$ a vector of factors. Alternatively, we could use the grouped fixed effects approach of bonhomme2015grouped. Finally, while our focus here is a panel comprising repeated time series observations on the same unit our approach could be applied to a network setting in which there is contemporaneous dependence across units. We leave these extensions to future work.