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.
176,848 characters · 14 sections · 66 citation commands
Bayesian Estimation of Panel Models under Potentially Sparse Heterogeneity
JEL CLASSIFICATION: C11, C23, C53, E20
KEY\ WORDS: Bayesian Analysis, Forecasting, Income Dynamics, Panel Data Models, Sparsity, Spike-and-Slab Priors
\thispagestyle{empty} \setcounter{page}{0}
Panel data models feature two types of coefficients: homogeneous coefficients that are common among all units $i=1,\ldots,N$, and heterogeneous coefficients that differ across units. In recent years, there has been a growing literature that focuses on the estimation of heterogeneous coefficients, which is challenging in environments in which the time series dimension is small. A natural way of adding information to the estimation of unit-specific coefficients is the use of prior distributions. The key insight in panel data applications is that one can extract information from the cross-section and equate the prior distribution with the (unobserved) cross-sectional distribution of unit-specific coefficients. The time series information for each unit can then be combined with the prior distribution, to generate more precise estimates of the heterogeneous coefficients. While in many models the cross-sectional coefficient distribution is in principle non-parametrically identified, in practice it is often desirable to impose parametric restrictions. In this paper, we consider a particular restriction: we assume that there is a core group of units that share a common parameter value, whereas the remaining units (deviators) have idiosyncratic values. If the core group contains a large number of units, then heterogeneity is sparse.
We assume that the heterogeneous coefficients are distributed according to a version of a spike-and-slab (S&S) prior: with probability $1-q$ the coefficient for unit $i$ equals the predominant value (spike), and with probability $q$ it deviates from the predominant value. This setup nests two important special cases: coefficient homogeneity ($q=0$) and full heterogeneity ($q=1$). We will say that heterogeneity is sparse if $q$ is close to zero and dense if $q$ is close to one. The probability $q$ and the parameters characterizing the distribution of the deviations are hyperparameters that are estimated from the panel data. Whether the concept of sparse heterogeneity is useful in practice, is an empirical question. We use the proposed prior distribution to estimate a dynamic model for labor earnings obtained from the Panel Study of Income Dynamics (PSID). We construct two types of samples from the PSID data: balanced panels that include individuals who have been working uninterruptedly for a long period of time, and an unbalanced panel that includes individuals with a variety of earnings and employment histories. The balanced panels do feature sparse heterogeneity, whereas the more diverse unbalanced panel is better described by full heterogeneity.
The paper makes the following contributions: first, we incorporate a version of an S&S prior, comprising a pointmass at zero (“spike”) and a Normal distribution around zero (“slab”) into a dynamic panel data framework to model coefficient heterogeneity. Such a prior has been used by Geweke1996 and, more recently, GiannoneLenzaPrimiceri2018, in the context of regressor selection problems. We extend the prior to handle heterogeneous variance parameters and show how to implement posterior inference using a Gibbs sampler. Second, we conduct a Monte Carlo study that compares the point estimation accuracy under the proposed (S&S) prior to that obtained by (incorrectly) imposing coefficient homogeneity ($q=0$) or full heterogeneity ($q=1$). If $q$ is close to zero or one, then estimates under the S&S prior are about as accurate as the estimates that impose the nearly-correct restriction. If $q$ is between 0.2 and 0.8 the estimates under the S&S prior clearly dominate the alternatives. Third, we show in our application that the concept of sparse heterogeneity is empirically useful.
In the empirical application we estimate income processes based on earnings data from the PSID. There is a long-standing debate in the literature whether earnings profiles are homogeneous or heterogeneous across individuals. Guvenen2007 labeled these specifications RIP (restricted income profiles) and HIP (heterogeneous income profiles), respectively. We offer a third option in our setup: sparsely heterogeneous income profiles, which could be abbreviated as SHIP. The literature on the estimation of idiosyncratic income processes has considered many different specifications. There are a few features that all of them share: the income process typically takes the form of an unobserved components model with deterministic and stochastic components. We consider three components: a linear function in experience, an autoregressive component with AR coefficient $\rho_i$, and an $iid$ component.\footnote{The unobserved components model is fitted to residuals from a regression of earnings on time fixed effects and some demographic characteristics.} There is typically time variation and cross-sectional variation in the income risk, i.e., the stochastic components. In our model we allow for time specific variance coefficients that also vary across individuals. While most of the literature assumes that the autocorrelation parameter of the persistent income component is common across units, we allow for heterogeneity.
The evidence for sparse heterogeneity depends on the construction of the panel data set that is used in the estimation. Using relatively small balanced samples, our posterior estimates imply that heterogeneity with respect to returns to experience and the autocorrelation of the persistent stochastic income component is sparse at best, with large core groups of individuals sharing identical coefficients and a small number of deviators. These estimates can be interpreted as evidence against a HIP version of the model. We also generate one-step-ahead forecasts for the entire cross-section of individuals using our baseline specification, a RIP specification, and a HIP specification. Using log predictive scores, the HIP version is clearly dominated, and the performance of our sparsely heterogeneous baseline specification is for most samples identical to that of the more restrictive RIP specification. Repeating the estimation on a large unbalanced panel, the parameter estimates lead to an almost fully heterogeneous specification that does not feature core groups of units with identical parameter values, which is evidence in favor of the HIP conjecture. In terms of one-step-ahead out-of-sample forecasting performance, however, even for the unbalanced panel the more parsimonious RIP specification leads to slightly better forecasts on average.
Third, we examine unit-level interval forecasts conditional on information from the full panel $Y_{1:N,0:T}$ accounting for shock and parameter uncertainty, the full panel ignoring parameter uncertainty, and the individual history $Y_{i,0:T}$ accounting for shock and parameter uncertainty, distinguishing between core group members and deviators. This comparison has an important message for the macroeconomic life-cycle model literature: subtle assumptions about the agents information set can have large effects on the income uncertainty that agents face and hence on their precautionary savings motive.
Our paper is related to several strands of literature. We are using a variant of an S&S prior, which was originally proposed for regressor selection problem. The prior proposed by MitchellBeauchamp1988 consisted of a point mass at zero (“spike”) and a uniform prior (“slab”). For computational reasons GeorgeMcCulloch1993,GeorgeMcCulloch1997 replaced the point mass by a Normal distribution with very small variance, and the uniform distribution by a Normal with large variance. Geweke1996 combined a pointmass at zero with a Gaussian distribution, which is the specification that was recently used by GiannoneLenzaPrimiceri2018 to examine whether regression models in economics tend to be sparse or dense. The regressor selection problem can be viewed as a special case of the fundamental problem of estimating the vector of means of a multivariate Normal distribution. JohnstoneSilverman2004 examined the large sample properties of Bayesian posteriors under a prior that is a mixture of a pointmass at zero and a continuous distribution that has fatter tails than a Normal distribution. The authors are particularly interested in sparse settings with a large number of zeros. The asymptotic concentration of posterior distributions in this environment is also studied by CastillovanderVaart2012, albeit under a somewhat different prior distribution.
Our paper is concerned with the estimation of potentially heterogeneous slopes and variance parameters in dynamic panel data models, which involves more complicated versions of the vector of means problem. While much of the panel data literature has traditionally focused on the estimation of the common or homogeneous parameter values, treating, for instance, heterogeneous intercepts as incidental nuisance parameters, the estimation of the heterogeneous coefficients can be important for treatment effect analysis or forecasting. In the context of our application, the prediction of idiosyncratic income profiles is essential for consumption and investment decisions and is at the core of every life-cycle model in macroeconomics. We are building on recent Bayesian approaches by ChamberlainHirano1999, GuKoenker2014,GuKoenkerJAE2016, LiuMoonSchorfheide2017, Liu2023, LiuMoonSchorfheide2018. The key point of departure from the earlier work is the use of an S&S type prior that lets us pick up patterns of either sparse or dense heterogeneity, bridging homogeneity and full heterogeneity.\footnote{The notion of group heterogeneity has received a lot of attention in the recent panel data literature, see BonhommeManresa2015 for the seminal paper and Zhang2023 for a Bayesian implementation that allows users to formulate priors about group memberships. Our paper considers a specification with one core group of units. In principle it could be extended to multiple core groups, but this extension is left for future research.}
There exists a large literature at the intersection of labor economics and macroeconomics on the estimation of idiosyncratic income processes. A longstanding debate in this literature is whether income profiles are homogeneous or heterogeneous. The heterogeneity typically refers the return to experience. Among others, MaCurdy1982, AboundCard1989, and Guvenen2007,Guvenen2009 provide evidence in favor of HIP, whereas LillardWeiss1979, Hause1980, Baker1997, Haider2001, Hryshko2012, and Hoffmann2019 favor the homogeneous RIP specification, which often involves a unit-root or highly persistent stochastic component. The importance of heterogeneity in income profiles, persistence of income processes, and the variance of income shocks is also emphasized by BrowningEjrnjesAlvarez2010 and BrowningEjrnjes2013. Coefficient heterogeneity in this literature is typically modeled through random effects. The latter two papers allow for these random effects to be correlated with each other, assuming they are determined by some underlying factors. While we do not allow for conditional heteroskedasticity in our model specification as, for instance, in MeghirPistaferri2004 or Hospido2012, we do allow for cross-sectional heterogeneity as in GuKoenker2014 and find it to be dense.
Traditionally, the earnings dynamics literature focuses on the estimation of (correlated) random effects distributions, characterizing the degree of coefficient heterogeneity in the population. Our Bayesian estimation approach generates such estimates along with estimates of the unit-level coefficients, under the assumption of potentially sparse heterogeneity. Whether or not the interest lies in the random effects distribution or the unit-level coefficients depends on the particular application. One novel contribution of our application below, is to also assess the fitted income-dynamics model by evaluating the pseudo-out-of-sample forecast accuracy for the individuals included in the panel data set, rather than a hypothetical population implied by the estimated model.
The remainder of the paper is organized as follows. Our approach of capturing potentially sparse heterogeneity through an S&S type prior distribution is described in the context of a vector of means estimation in Section (ref). This prior is then incorporated into two linear panel data models with heterogeneous intercept, slope coefficients, and shock standard deviations in Section (ref). Results from two Monte Carlo experiments are summarized in Section (ref) and the empirical analysis is presented in Section (ref). Finally, Section (ref) concludes. Details on the posterior samplers for the Monte Carlo study and the empirical analysis are provided in the Online Appendix. The Appendix also contains additional empirical results.
Throughout this paper, we pursue a Bayesian approach, place a prior distribution over the parameters, and generate draws from the posterior distribution. In the context of estimating a vector of means of a multivariate Normal distribution, our approach to sparse heterogeneity takes the following form. Consider the model
equipped with a prior distribution that is a mixture of a pointmass at zero and a Normal distribution:
This prior distribution has been used by Geweke1996 and, more recently, by GiannoneLenzaPrimiceri2018 in the context of regressor selection problems. The sparsity of the $\delta = [\delta_1,\ldots,\delta_N]'$ vector is determined by the hyperparameter $q$, which controls the probability of the coefficient for unit $i$ deviating from the common value. Thus, $1-q$ is the probability assigned to the pointmass at zero and hence the height of the “spike.” The second hyperparameter, $v_{\delta}$, controls the width of the “slab.” The larger $v_{\delta}$, the larger deviations from the common parameter on average. We will refer to units for which $\delta_i=0$ as the core group, whereas the other units are deviators.\footnote{ Alternatively, one could think of two groups: a $\delta_i=0$ group and a $\delta_i \not=0$ group.} The prior can be generalized in a straightforward manner by replacing the ${\cal N}(0,v_\delta)$ distribution with a ${\cal N}(\underline{\delta},v_\delta)$ distribution and treating the $\underline{\delta}$ as a third hyperparameter.
Let $Y_{1:N}$ denote $\{ y_1,\ldots, y_N\}$. Conditional on the hyperparameters $(q,v_\delta)$, the posterior distribution is also a mixture of a pointmass at zero and a Normal component:
where
The odds in favor of a Gaussian component are an increasing function of the squared deviations from zero, $y_i^2$. The posterior mean is given by $\mathbb{E}[\delta_i|Y_{1:N}] = q^* \delta_i^*$, which is never exactly equal to zero. However, it can be shown that the posterior median, denoted by $\mbox{med}(\delta_i|Y_{1:N})$, has the property that it can be exactly equal to zero for small values of $y_i^2$: $\mbox{med}(\delta_i|Y_{1:N}) = 0$ for $-c(\delta_i^*,q^*,v_\delta^*) \le y_i \le c(\delta_i^*,q^*,v_\delta^*)$; see, for instance JohnstoneSilverman2004.
The hyperparameters are determined by the marginal likelihood function
We will combine this marginal likelihood with a prior $p(q,v_\delta)$ for the hyperparameters and conduct posterior inference. For the posterior computations it is convenient to introduce the auxiliary variables $z_i$ such that $z_i=0$ if $\delta_i=0$ and $z_i=1$ if $\delta_i \not=0$. We will say that unit $i$ is a member of the core group if $z_i=0$ and it is a deviator if $z_i=1$. Then,
Given a prior for $(q,v_\delta)$, draws from the posterior distribution can be obtained via Gibbs iterating over the conditional posteriors $Z_{1:N}|(Y_{1:N},q,v_\delta)$, $q|(Y_{1:N},Z_{1:N},v_\delta)$ and $v_\delta|(Y_{1:N},Z_{1:N},q)$.
Notice that the argmax of $p(Y_{1:N},Z_{1:N}|q, v_\delta)$ with respect to $(Z_{1:N},q,v_\delta)$ satisfies:
where $\mathbb{I}\{x=a\}$ is the indicator function that is equal to one if $x=a$ and equal to zero otherwise. Thus, $\hat{z}_i=1$ if the likelihood ratio of $\delta_i \sim {\cal N}(0,\hat{v}_d)$ versus $\delta_i=0$ exceeds one, $\hat{q}$ corresponds to the fraction of $\hat{z}_i$s that are equal to one, and $\hat{v}_d$ is the variance of the observations attributed to the Gaussian component of the model (minus one).
The prior distribution in ((ref)) nests coefficient homogeneity ($q=1$) and full coefficient heterogeneity ($q=1$) as special cases. In practice, whether $q$ is close to zero or one is partly determined by the researcher when (s)he constructs the estimation sample. In our empirical application we use the raw panel data set to construct balanced panels in which units are fairly homogeneous and an unbalanced panel that tries to include as many units as possible. The former samples will feature sparse heterogeneity in some dimensions, whereas the latter sample will lead to $q$ estimates that are close to one.
We consider two dynamic linear panel data models in this paper. The first model, $M_1$, which will be used in the Monte Carlo analysis in Section (ref), takes the form of a dynamic linear regression. The second model, $M_2$, is an unobserved components model that can be written in state-space form and is used for the empirical analysis. In the remainder of this section we present the specifications of $M_1$ and $M_2$ and sketch the posterior sampler for the state-space model $M_2$, which nests the simulation algorithm for the regression model $M_1$ as a special case.
{\bf Linear Dynamic Panel Regression Model $M_1$.} Under $M_1$ the dependent variable $y_{it}$ evolves according to
This specification allows for cross-sectional heterogeneity in the intercept $\alpha$, the autoregressive coefficient $\rho$, and the innovation standard deviation $\sigma$ through the discrepancies $\delta^{\alpha}_i$, $\delta^{\rho}_i$, and $\delta^{\sigma}_i$. Units $i$ for which the discrepancies are equal to zero, share a common intercept, autocorrelation parameter, and shock standard deviation, respectively. The remaining units deviate from the parameters of the core group(s) and have heterogeneous coefficients. To induce sparsity, we use the prior in ((ref)) for $\delta_i^\alpha$ and $\delta_i^\rho$, equipping $q$ and $v_\delta$ with superscripts $\alpha$ and $\rho$. Because of the non-negativity of variances, we use a slightly different prior for $\delta_i^\sigma$:
where $IG$ is the Inverse Gamma distribution. We impose that $\mathbb{E}[\delta_i^\sigma]=1$ and reparameterize the prior distribution in terms of its variance $v_{\delta^\sigma}$ which leads to\footnote{The model could be restricted by requiring that $q^\alpha=q^\rho=q^\sigma=q$ and that the indicator functions $z_i^\alpha=z_i^\rho=z_i^\sigma = z_i$ are identical. This restriction implies that members of the core group have identical $\alpha$, $\rho$, and $\sigma^2$ parameter values, and deviators differ in all three dimensions.}
We use the following priors for the common parameters $\alpha$, $\rho$, and $\sigma^2$ of the core group(s):
For the hyperparameters of the S&S priors we use
where $B(a,b)$ is the Beta distribution. We define the parameter vector $\theta$ as
and stack the hyperparameters that index the prior distributions in ((ref)) and ((ref)) into the vector
We denote the prior density of $\theta$ conditional on the hyperparameters as $p(\theta|\lambda,M_1)$.
{\bf State-Space Model $M_2$.} This model contains a latent variable $s_{it}$ and takes the form of a state-space model, which is given by
The main difference between ((ref)) and ((ref)) is the presence of the latent state $s_{it}$. ((ref)) can be interpreted as the measurement equation of the state-space model, and ((ref)) is the state transition. Moreover, we let the shock standard deviations $\sigma_{u,t}$ and $\sigma_{\epsilon,t}$ vary over time. The parameter $\alpha$ is now a $k$-dimensional vector that interacts with the vector of regressors $x_{it}$.
For the $\delta_i$ discrepancies we use the priors ((ref)) and ((ref)) from above with the modification that $v_\delta^\alpha$ is a covariance matrix for $k>1$. The priors for $\alpha$ and $\rho$ are the same as in ((ref)), with the understanding that $\underline{v}_\alpha$ is also a covariance matrix. For the shock variances we now use independent IG priors:
The prior for the hyperparameters of the S&S prior in ((ref)) is modified as follows: because $\alpha$ is a vector for $k>1$ the IG distribution for $v_\delta^\alpha$ is replaced by an Inverse Wishard distribution $IW(\underline{\nu}_{\delta^\alpha},\underline{\Psi}_{\delta^\alpha})$. Moreover, the priors for $q^\sigma$ and $v_{\delta^\sigma}$ are replaced by separate priors for $q^{\sigma_u}$, $q^{\sigma_\epsilon}$, $v_{\delta^{\sigma}_u}$, and $v_{\delta^{\sigma}_\epsilon}$. Finally, we need a prior distribution for the initial state $s_0$. We assume that $s_0 \sim {\cal N}(\mu_{s_0}, v_{s_0})$, where
Notice that the state-transition equation ((ref)) induces the distribution of $s_1,\ldots,s_T$ conditional on $s_0$. This leads to the following definitions of $\theta$ and $\lambda$:
{\bf Restricted Model Specifications.} For $M_1$, we consider a homoskedastic version in Experiment 1 which is obtained by setting $q^\sigma=0$ and $\delta_i^\sigma=0$. Moreover, we also estimate a model with homogeneous coefficients which is obtained by setting $q^\alpha=q^\rho = 0$ and $\delta_i^\alpha=\delta_i^\rho=0$. Under this restriction every unit $i$ belongs to the core group. A fully heterogeneous model specification is obtained by letting $q^\alpha=q^\rho=1$, meaning there no longer are core groups of units that share identical coefficients. In Experiment 2 we consider the sparsely heteroskedastic model specification and the restricted versions that are obtained by either setting $q^\alpha=q^\rho=q^\sigma=0$ (homogeneous) or $q^\alpha=q^\rho=q^\sigma=1$ (fully heterogeneous). The empirical analysis is based on model $M_2$. Here we consider a homoskedastic version, a heteroskedastic version with homogeneous $\alpha$ and $\rho$ coefficients, and a heteroskedastic version with fully heterogeneous $\alpha$ and $\rho$ coefficients. As for $M_1$ the restricted specifications are obtained by setting the appropriate $q^l$, $l \in \{\alpha,\rho,\sigma_u,\sigma_\epsilon\}$, either equal to zero or one.
{\bf Posterior Sampling.} We outline a few important aspects of the posterior sampler for $M_2$. The sampler for $M_1$ is a special case in which there is no latent variable $s_{it}$. Detailed descriptions of the posterior samplers are provided in the Online Appendix.\footnote{Posterior computations for the restricted $q=0$ and $q=1$ versions of model $M_2$ have been conducted by NakataTonetti2015. Their paper shows through a simulation study that Bayes estimates are generally more accurate than the GMM/minimum distance estimates that are widely-used in the empirical literature.}
We use a Gibbs sampler to generate draws from the posterior distribution of the parameters $(\theta,\delta^\alpha,\delta^\rho,\delta^{\sigma_u},\delta_{\sigma_\epsilon})$ and the latent states $S_{i,0:T} = \{s_{i0},\ldots, s_{iT} \}$. To implement the Gibbs sampler we use indicator variables $z_i^l$ such that $z_i^l=0$ if $\delta_i^l=0$ and $\delta_i^l=1$ otherwise, where $l \in \{\alpha,\rho,\sigma_u,\sigma_\epsilon\}$. Because we found $\delta_i^\alpha$ to have a strong {\em a posteriori} correlation with $S_{i,0:T}$ we do not use a Kalman filtering/simulation smoothing step. Instead, we write $M_2$ as
where $Y_i$, $X_i$, and $U_i$ are matrices that stack $y_{it}$, $x_{it}'$, and $u_{it}$, respectively. $D_i$ is a diagonal matrix with elements $\sigma_{u,t}\sqrt{\delta_{i,u}}$ and the $T \times (k+T)$ design matrix $W_i$ comprises the regressor matrix $X_i$ and time dummies. For instance, if $k=1$ and $x_{it}=1$ such that $\alpha$ is an intercept,
and $\beta_i' = [\delta^{\alpha}_{i}, \; s_{i1}, \; ... \;, s_{iT}]$ stacks the $\alpha$ discrepancy and the latent states. The prior for $\beta_i$ combines the prior for $\delta_i^\alpha$ with the prior for the latent states generated by the state-transition equation ((ref)). In case of $K=1$ this leads to
where $D_{i,\epsilon}$ is a diagonal matrix with elements $\sigma_{\epsilon,t} \sqrt{\delta^{\sigma}_{i,\epsilon}}$. The $(t,\tau)$ elements of the $T \times T$ prior covariance matrix $\underline{V}_{s_i}$ can be calculated as follows:
with the understanding that $\underline{V}_{s_i}(0,0)=v_{s_0}$. The conditional {\em posterior} distribution of $\beta_i$ is also a mixture of Normals.
In the Monte Carlo experiment we generate panels of observations $Y_{1:N,0:T}$ from model $M_1$ and compute estimates $\hat{l}_i(Y_{1:N,0:T})$ of $l_i = l + \delta_i^l$, where $l \in \{\alpha,\rho\}$. We evaluate the estimates based on a quadratic loss function and report a Monte Carlo approximation of the compound risk
The Monte Carlo experiment is based on the following algorithm:
We compare the performance of four posterior mean estimators. The first one is the spike-and-slab (S&S) estimator defined as \[ \hat{l}_i^* = \mathbb{E}[l + \delta_i^l|Y_{1:N,0:T},\lambda]. \] This estimator treats the slab probability $q$ as unknown, but conditions on the hyperparameters $\lambda$. Second, the estimator $\hat{l}_i^0$ sets $q^l=0$, $l \in \{\alpha,\rho,\sigma\}$, and thereby imposes homogeneity on all of the coefficients. Third, we compute an estimator $\hat{l}_i^1$ for a fully heterogeneous estimator that is obtained by setting $q^l=1$, $l \in \{\alpha,\rho,\sigma\}$. Finally, we report the performance of the oracle estimator that is based on knowing the “true” value of $\theta$: \[ \hat{l}_i^o = \mathbb{E}[l + \delta_i^l|Y_{1:N,0:T},\theta]. \] We consider a homoskedastic design (Section (ref)) and a heteroskedastic design (Section (ref)).
For the first Monte Carlo experiment we consider a homoskedastic design, i.e., in model ((ref)) we set $q^\sigma=0$ and $\delta_i^\sigma = 1$ for all $i$. Baseline choices for the parameter vector $\theta$ and the hyperparameters are summarized in Table (ref). First, consider the elements of $\theta$. We fix $\alpha=1.0$, set $\rho=0.6$, and let $\sigma^2=0.8$. In the experiments below we consider various values for $q^\alpha=q^\rho=q$ and $v_\delta^\alpha$, which control the height of the spike and the width of the slab of the $\delta_i^\alpha$ distribution. The parameter $v_{\delta^\rho}$ is set to 0.09. The parameter $v_{\delta^{\sigma}}$ is only relevant for the heteroskedastic design.
Second, consider the elements of $\lambda$. We set the variance of $\alpha$, denoted by $\underline{v}_\alpha$, equal to one. Moreover, we let $\underline{v}_\rho =0.25$. The hyperparameters for the prior for $\sigma^2$ imply that $\mathbb{E}[\sigma^2]=1$ and $V[\sigma^2]=1/4$. Our choices of $(a,b)$ imply that $q \sim U[0,1]$. We choose $\nu_{\delta^\alpha}=6$ and $\tau_{\delta^\alpha}=4$, which implies $\mathbb{E}[v_{\delta^\alpha}]=1$ and $\mathbb{V}[v_{\delta^\alpha}]=1$. This distribution covers the range of $v_{\delta^\alpha}$ values in Table (ref). Suppose we choose $\nu_{\delta^\rho}=6$ and $\tau_{\delta^\alpha}=2$. Then, $\mathbb{E}[v_{\delta^\alpha}]=1/2$ and $\mathbb{V}[v_{\delta^\alpha}]=1/4$. The true value of $v_{\delta^\rho}=0.09$ can be viewed as a draw from the left tail of this distribution.
Using Algorithm (ref), we generate $N_{sim}=100$ panel data sets of size $N=500$ and $T=8$. For each panel data set we generate 5,000 parameter draws from the posterior distribution, discard the first 2,500 draws, and use the remaining draws to compute the posterior mean estimates $\hat{l}_i(Y_{1:N,0:T})$. Simulation results are summarized in Table (ref). In the top half of the table we report mean squared errors (MSEs) for the estimates $\hat{\alpha}_i$ of $\alpha_i = \alpha+\delta_i^\alpha$ for various combinations of $(v_{\delta^\alpha},q)$ and in the bottom half MSEs for $\hat\rho_i$. During the estimation we do not restrict $q^\alpha$ to be equal to $q^\rho$. Thus, the members of the core $\alpha$ and $\rho$ groups are not required to be the same. Formally, we do not impose the restriction that $z_i^\alpha = z_i^\rho$ for all $i$.
For $q=0$ the Bayes estimator with S&S achieves the performance of the estimator $\hat{l}_i^0$ that imposes coefficient homogeneity, and for $q=0.8$ and $q=1$ its compound risk is approximately equal to the risk associated with the estimate $\hat{l}_i^1$ that treats $l_i$ as fully heterogeneous. For values of $q \in \{0.2, 0.4, 0.6\}$ the S&S estimator achieves a lower MSE than the estimators that impose full homogeneity or full heterogeneity. Overall, $\hat{l}_i^*$ is equal or smaller than the risk of the best estimator among $\hat{l}_i^0$ and $\hat{l}_i^1$. The risk differentials between the fully homogeneous and fully heterogeneous estimators is increasing in $v_{\delta}^\alpha$ and the risk patterns are similar for the estimation of $l=\alpha_i$ and $l=\rho_i$. As a benchmark, we also result the risk of the oracle estimator that is obtained by conditioning on the “true” $\theta$. By construction the oracle estimator dominates the Bayes estimator $\hat{l}^*_i$, but overall the risk differentials are relatively small which means that $\theta$ is generally well estimated.
In the top and center panels of Figure (ref) we plot histograms of the cross-sectional distribution of the posterior mean estimates $\hat{\alpha}_i^*$ and $\hat{\rho}_i^*$ obtained conditional on data sets simulated with $q=0.2$ and $q=0.8$, respectively. Unlike the results reported in Table (ref), the histograms are generated from a single Monte Carlo iteration and not averaged across $N_{sim}$ runs. For $q=0.2$ the histogram of posterior means inherits the spike-and-slab shape of the prior distribution with a distinct spike near the “true” values of $\alpha$ and $\rho$, which are 1.0 and 0.6, respectively; see Table (ref). For $q=0.8$ the spike in the distribution of posterior means is much less pronounced and there is a lot more mass away from $\alpha=1.0$ and $\rho=0.6$. Finally, we plot the posterior densities for the slab probabilities $q^\alpha$ and $q^\rho$. While there is considerable uncertainty about $q$, the densities peak near the “true” values of $0.2$ and $0.8$, respectively.
In the second Monte Carlo design we allow for heteroskedasticity by setting $q^\sigma=q$ and sampling $\delta_i^\sigma$ from ((ref)). Recall that we imposed $\mathbb{E}[\delta_i^\sigma] =1$. To generate the data we let $v_{\delta^\sigma}=\mathbb{V}[\delta_i^\sigma]=1$. The hyperparameters for the $v_{\delta^\sigma}$ prior imply that $\mathbb{E}[v_{\delta^\sigma}] = 1$ and $\mathbb{V}[v_{\delta^\sigma}] = 1/4$. In all other aspects, the experiment is configured in the same as the homoskedastic experiment in Section (ref).
Results for the MSE of $\hat{\alpha}_i$ are presented in Table (ref). We consider two S&S estimators: (hetsk) is based on the heteroskedastic version of the benchmark model which is the correctly specified model in light of the data generating process (DGP). (homsk) is based on the misspecified homoskedastic version of the benchmark model, i.e., $\delta_i^\sigma=0$. The results are similar to the ones obtained in the first experiment. For $q \in \{0.2,\ldots,0.8\}$ the S&S(hetsk) estimator dominates the estimators obtained under full homogeneity, $q=0$, and full heterogeneity, $q=1$. The S&S(homsk) generally performs poorly, except for $q=0$, when there is no heteroskedasticity. Thus, ignoring the heteroskedasticity generates imprecise estimates also of the conditional mean parameters.
Figure (ref) extends the results shown in Figure (ref) to the heteroskedastic design. As before, we plot histograms of the cross-sectional distribution of the posterior mean estimates, which now include $\hat{\sigma}_i^*$. As for the homoskedastic design, for $q=0.2$ the cross-sectional distributions of the estimators inherit the spike-and-slab shape of the prior distribution, which is more pronounced for $\hat{\alpha}_i^*$ and $\hat{\rho}_i^*$ than for $\hat{\sigma}_i^{2*}$. The distributions of posterior means are generally centered near the “true” values of the homogeneous parameters.
We now use model $M_2$, comprising ((ref)) and ((ref)), to estimate income profiles from the PSID data. Model $M_2$ captures some of the key features of the empirical specifications considered in the literature over the past four decades; see for instance the model specification discussion in the survey paper BrowningEjrnjes2013. These features include a deterministic income component $x_{it}'\alpha$, and persistent and transitory stochastic components, $s_{it}$ and $u_{it}$, respectively. We allow for potential heterogeneity in the $\alpha$ vector, the autocorrelation $\rho$ of the persistent component, and the shock standard deviations.
More specifically, we let
where $h_{it}$ represents experience of individual $i$ in period $t$, defined as age minus years of education minus six. We partition the $2\times 1$ vector $\alpha = [\alpha_0,\alpha_1]$. We also divide experience by 10 to scale $\alpha_1$ and improve the numerical performance of the posterior sampler. If $\delta_i^\alpha=0$ for all $i$, then income profiles are homogeneous, which was previously labeled RIP and has been the preferred specification of LillardWeiss1979, Hause1980, Baker1997, Haider2001, Hryshko2012, and Hoffmann2019. The case of $\delta_i^\alpha\not=0$ for all individuals, corresponds to the HIP specification considered in the literature and is favored by, for instance, MaCurdy1982, AboundCard1989, and Guvenen2007,Guvenen2009. Unlike the previous literature, our model considers a third option, under which heterogeneity is potentially sparse: $\delta_i^\alpha\not=0$ for some, but not all, individuals.
Income has a persistent and a transitory stochastic component. The persistent component is given by an AR(1) process $s_{it}$ with autocorrelation coefficient $\rho$ and the transitory component is represented by the $iid$ process $u_{it}$. Through $\delta_{i,u}^\sigma, \, \delta_{i,\epsilon}^\sigma \, \not=0$ we allow income risk to be heterogeneous across agents. In addition, the standard deviations $\sigma_{u,t}$ and $\sigma_{\epsilon,t}$ are time dependent, which can capture business cycle fluctuations and trends in idiosyncratic risk. Finally, we allow for idiosyncratic autocorrelations $\rho_i$ which is typically not done in the literature, an exception being BrowningEjrnjesAlvarez2010 and BrowningEjrnjes2013.
There are a few model features that have been explored in the literature that we do not include in $M_2$ to keep the model parsimonious: a third stochastic component that features a unit root, idiosyncratic autoregressive conditional heteroskedasticity (see, for instance, MeghirPistaferri2004 for both features), time variation in the returns to experience (e.g., Haider2001), or age effects in the variances of the idiosyncratic shocks (e.g., Hoffmann2019). In addition to $M_2$ itself we consider three restricted versions: one without heteroskedasticity, one with homogeneous slope parameters (RIP), and one with fully heterogeneous slope parameters (HIP). A summary is provided in Table (ref).
Some authors have emphasized, for instance, cohort effects on earnings, education effects on earnings profiles, or age-dependent persistence in earnings. While we are not explicitly controlling for this type of heterogeneity through sample splitting or interaction terms, it can be captured {\em a posteriori} by estimates of $\delta_i$ that are different from zero.
Much of the previous literature has used (simulated) minimum-distance estimation of a subset of the parameters that we collect in the vector $\theta$, which involves minimizing the discrepancy between sample moments and model-implied moment. Typically, these moments are time-dependent autocovariances that are computed by cross-sectional averaging across all individuals in the panel. In a few instances, researcher have also used moments computed from unit-specific time series; see BrowningEjrnjes2013. Our likelihood-based estimation approach aims at tracking individual earnings histories, which implicitly encompasses both time series and cross-sectional moments. This approach has several advantages. In particular, we are easily able to handle non-Normalities introduced by the S&S prior and the cross-sectional heteroskedasticity. Morever, in addition to parameters that characterize the distribution of coefficients across agents, we can also make inference about the unit-specific coefficients themselves, and subsequently generate predictions about future income for individual units.
The empirical analysis is based on the PSID. We use raw data from 1968 to 1997, starting from 32,465 observations for 2,052 individuals. Up to 1997 the PSID data are annual, whereas subsequently the survey was conducted only biennially. The data set used for estimation is constructed following the conventions of the literature, e.g., MeghirPistaferri2004. The dependent variable is log real income of all sources, including: wages and salaries; bonuses, overtime, and/or commissions; income from professional practice and trade; the labor part of farm income and unincorporated business income; the labor portion of income from farming, market gardening, roomers and boarders. Nominal wages are converted into real wages, using the PCE deflator (base = 1993). We take account of the fact that the measure of income refers to the previous calender year. We eliminate observations with outlying earnings records, defined as a change in log real earnings greater than two or less than minus one.
As is common in the literature, we restrict the sample to male heads of household and keep individuals between ages 25 and 55.\footnote{If the household age is inconsistent across years, we compute the mode of the differences between reported age and year across periods. Age is then defined as the year plus the mode.} To avoid oversampling low-income households, members from the Survey of Economic Opportunity (SEO) are removed. We also drop the members of the Latino sample added in 1990 and self-employed individuals. We keep individuals who are in the labor force, have positive labor income, and have been in the sample for at least nine years. We remove individuals with missing race, marital status, region, state, and education information. In a first-stage regression we project log income on year dummies, and indicators for race, educational group,\footnote{High School dropouts: (those with less than 12 grades of schooling); High School graduates (those with at least a High School diploma, but no College degree); College graduates (those with a College degree or more)} region of residence, and marital status. The variable $y_{it}$ is the residual from this regression. We refer to the resulting panel data set as the full sample.
Based on the full sample, we construct estimation samples. The first set of estimation samples are balanced panels, denoted by ${\cal B}_{T,\tau}$, and used for the analysis in Sections (ref) and (ref). The balanced panels have a time dimension of $T+1$ and include only individuals for which $T+1$ observations are available. Typically, we use the first $T$ time periods in these samples for estimation and reserve the last time period, $t=T+1$, for one-step-ahead forecast evaluation. The index $\tau$ refers to the year that corresponds to period $T+1$. For instance, the sample ${\cal B}_{20,1991}$ starts in 1971 ($t=1$) and ends in 1991 ($t=20+1$).
In addition to the balanced sample, we also estimate model $M_2$ and its variants based on an unbalanced sample ${\cal U}$ in Section (ref). The unbalanced sample is similar to the full sample described above, except that we require that individuals have at least nine {\em consecutive} observations. Figure (ref) provides some information about the composition of ${\cal U}$. The left panel of the figure shows the cross-sectional dimension as a function of time. In 1968 it covers approximately 500 units. The number of units rises above 1,100 in 1988 and then declines to 800 in 1995. The right panel provides information about how long units stay in the panel. The modal duration is 10-12 years, but approximately 330 units stay for more than 20 years. {\em A priori}, we expect coefficient heterogeneity to be more prevalent in the unbalanced sample ${\cal U}$ than in the balanced samples ${\cal B}_{T,\tau}$.
Table (ref) summarizes the marginal prior distributions for the elements of $\theta$. We are assuming that the elements are {\em a priori} independent, such that the joint distribution is the product of the marginals. In addition to the $(a,b)$ parameterization of the prior, we also report the 5th and 95th percentiles, which delimit the 90% equal-tail probability credible intervals.
The prior for $\alpha$ is centered at zero with a variance of one. Recall that we divided experience by ten, which means that a value of $\alpha_1=0.2$ implies a return to a year of experience of $0.2/10 = 0.02$, or two percent. Thus, the 90% {\em a priori} credible interval covers all reasonable and some unreasonable values for the return to experience. The marginal distribution of $\rho$ is centered at 0.8, representing an {\em a priori} belief that the idiosyncratic income process $s_{it}$ is fairly persistent. The variance of one implies that the prior does not have a lot of curvature in the range of 0.6 to 1.0, which covers the estimates obtained in previous studies. We do allow for explosive dynamics {\em a priori}, but turn out to be extremely unlikely under the posterior.
The 90% credible interval for the innovation variances reaches from 0.02 to 0.12. Note that our prior implies that they are independent across time periods, because in our panel data set we want to identify the variances from the cross-sectional information. Our prior for $\mu_{s_0}$ and $v_{s_0}$ implies that $s_{i0}$ has mean zero, which is the long-run mean of $s_{it}$ provided that $\rho$ is less than one. Moreover, $s_{i0}$ falls in the interval -1.0 to 1.0 with high probability.
The prior for the probabilities of not belonging to the core group, $q^l$, are uniform on the unit interval. The prior for the deviation variance $v_{\delta_\alpha}$ implies that the covariance matrix of the $\alpha$ deviations is diagonal with values 0.24 and 0.048. At 0.048, the 95% credible interval for $\delta_{\alpha_1}/10$ ranges from roughly -0.044 to 0.044.\footnote{The credible interval for the off-diagonal element of $v_{\delta^\alpha}$ ranges from -0.1 to 0.1.} The prior for $v_{\delta^\sigma_l}$ implies that an individual at the 5th percentile at the prior has an idiosyncratic shock variance that is one fourth of that of someone at the 95th percentile.
{\bf Idiosyncratic Parameters.} Of particular interest are the estimates of the idiosyncratic parameters and the degree of heterogeneity they exhibit. The results are summarized in Figure (ref). The figure shows posterior estimates for $\alpha_{0i}$, $\alpha_{1i}$, $\rho_i$, $\delta_{i,u}^\sigma$, and $\delta_{i,\epsilon}^\sigma$ for two samples ${\cal B}_{T,\tau}$. The panels in the left column correspond to $\tau=1988$ and the panels in the right column to $\tau=1991$. In both cases $T=20$. We depict posterior medians (solid blue lines), means (dashed red lines), and 90% credible intervals (light blue bands). Units are sorted based on posterior median estimates. As discussed in Section (ref), posterior median estimates of $\delta_i^l$ can deliver exact zeros. To highlight the potential sparsity implied by the posterior medians of $\delta_i^l$, we plot $\mbox{med}(l|Y_{1:N,1:T})+\mbox{med}(\delta_i^l|Y_{1:N,1:T})$ instead of $\mbox{med}(l+\delta_i^l|Y_{1:N,1:T})$ in the first three rows for $l \in \{\alpha_0,\alpha_1\}$.
Sparse heterogeneity manifests itself through a long flat section of the posterior median line. For instance, the $\tau = 1988$ sample consists of 82 observations. The number of units with $\mbox{med}(\delta_i^l|Y_{1:N,0:T})=0$ is 74 for $\alpha$. Thus, the $\alpha_i$s feature sparse heterogeneity. From examining the posterior distribution of the $\alpha_{1i}$s, which are divided by 10 to undo the rescaling of experience in ((ref)), we find that some units have a return experience that is up to 0.5 percentage points less than that of the core group, whereas others have a return that is up to 1.5 percentage points above that of the core group. The $\tau=1991$ sample has more heterogeneity in the $\alpha$ dimension, with only 32 out of 88 units having $\mbox{med}(\delta_i^\alpha|Y_{1:N,0:T})=0$. Unlike the posterior median curves, the posterior mean curves have generally no flat segments, because the posterior means of $\delta_i^\alpha$ are never exactly equal to zero.
For the autocorrelation and the two shock variances a different picture emerges. For $\rho$ essentially all units have $\mbox{med}(\delta_i^l|Y_{1:N,0:T})=0$, meaning that the autocorrelations are homogeneous, whereas for the variances essentially all of the units have $\mbox{med}(\delta_i^l|Y_{1:N,0:T})\not=1$ (recall the different centering of the variance discrepancies). Thus, the $\sigma^2_i$s are fully heterogeneous, meaning there is strong cross-sectional heteroskedasticity, a point that is also emphasized, for instance, by GuKoenker2014.
Figure (ref) shows posterior densities for the sizes of the $\alpha$ and $\rho$ core groups, defined as $N-\mathbb{E} \big[ \sum_{i=1}^N z_i^l \, | \, Y_{1:N,0:T} \big]$. The posterior modes of the $\rho$ core group sizes are 80 (out of 82) and 82 (out of 88), respectively. For the 1988 sample the posterior histogram for the size of the $\alpha$ group peaks at around 70, but the probability mass is spread out between zero (homogeneity) and 82 (full heterogeneity). In the 1991 sample, the probability mass for the $\alpha$ core group size shifts to the left and much of it lies between 0 and 30. We do not plot the distribution of the $\sigma$ group sizes, as it is clear from Figure (ref) that there are no core groups and the coefficients are fully heterogeneous.
{\bf Posterior Distribution of $\theta$.} Table (ref) summarizes the posterior median and 90% highest posterior density (HPD) estimates for the elements of the $\theta$ vector. The point estimates and credible intervals for $q^l$ line up with the information in Figures (ref) and (ref). The posterior medians of $q^{\sigma_u}$ and $q^{\sigma_\epsilon}$ are close to one and those for $q^{\rho}$ are close to zero. The posteriors for $q^\alpha$ are more spread out. For the 1988 sample the 90% credible set ranges from 0 to 0.84 and for 1991 from 0.53 to 1.
As in Figure (ref), we rescale the $\alpha_1$ (and $v_{\delta^{\alpha_1}}$) entries in the table to undo the division of experience by 10 in the definition of $x_{it}$ in ((ref)) and to facilitate its interpretation. For the 1988 sample the posterior median estimate of the return to an additional year of experience is 1.8%. It is 1.4% for the 1991 sample. Both the time-series and the cross-sectional variation contributes to the identification of the return to experience. Examining the cross-sectional distribution of experience at the beginning of the 1988 sample, we find that the 25th percentile is seven years of experience, the median is ten years, and the 75th percentile is thirteen years. For the 1991 sample these numbers are very similar.
The size of the discrepancy $\delta_i^\alpha$ from the $\alpha$ value of the core group is controlled by the matrix $v_{\delta^\alpha}$. We report posterior summary statistics for the diagonal elements. For the 1988 sample, the posterior credible interval for $v_{\delta^\alpha_0}$ ranges from 0.03 to 0.21. The interval for $v_{\delta^\alpha_1}$ (divided by 100) extends from .0001 to .0004. To put these estimates into perspective we provide a comparison to Guvenen2009, who reports a $v_{\delta^{\alpha_0}}$ estimate of .022 with a standard error of .074, and a $v_{\delta^{\alpha_1}}$ estimate of .00038 with a standard error of 0.0008. Hryshko2012 reports $\hat{v}_{\delta^{\alpha_1}}=.0004$. These estimates are based on a different subset of the PSID data and obtained from a different estimation procedure. Our $v_{\delta^{\alpha_0}}$ estimate is slightly higher, and the $v_{\delta^{\alpha_1}}$ estimate is slightly lower than those obtained in the two earlier papers, but intervals overlap.
The HPD interval for the autocorrelation parameter $\rho$ ranges from 0.91 to 1.00 for the 1988 sample and from 0.84 to 0.95 for the 1991 sample. Guvenen2009 obtains $\hat{\rho}=0.988$ for RIP specification with homogeneous coefficients, and 0.82 for the HIP model with fully heterogeneous coefficients. Our posterior medians of 0.97 and 0.9 are in between those estimates. The $\tau=1988$ features less heterogeneity in the experience profiles and is therefore closer to Guvenen's RIP model, whereas the 1991 one sample exhibits more heterogeneity in the $\alpha_i$s. Thus, we are able to reproduce the finding that homogeneous coefficient specifications are associated with more persistence in the stochastic component that is persistent. Because the estimated $q^\rho$ is close to zero, the sample contains hardly any information about $v_{\delta^\rho}$ and the posterior is close to the prior.
The estimated mean of the initial state $s_{i0}$ is zero. Using the posterior median for $v_{s_0}$ from the 1988 sample and ignoring the estimation uncertainty, the standard deviation of $s_{i0}$ is 0.17. A second back-on-the-envelope calculation using posterior median estimates of $\rho$, $v_{\delta_\epsilon^\sigma}$ and the average level of $\sigma^2_{\epsilon,t}$ in Figure (ref) below, leads to an unconditional standard deviation of $s_{it}$ of 1.25. Thus, focusing on a particular cohort, over time the income inequality due to the persistent stochastic component increases. This also shows up in the cross-sectional dispersion of the posterior mean estimates $\hat{s}_{it}$, when comparing period $t=1$ to $t=T=20$: the estimates range from -0.2 to 0.2 at the beginning of the sample, and from -1 to 1 at the end of the sample.
Using our estimates, we can decompose income inequality in a cohort into the contribution of return to experience heterogeneity and heterogeneity in the {\em ex post} realization of the persistent stochastic process $s_{it}$. To construct the decomposition, we simulate a cohort of $N=10,000$ individuals for $T=20$ periods from the model $M_2$ using the posterior mean estimates presented in Figure (ref) and Table (ref). The cohort starts with an experience of $h_{i1}=1$. We subsequently conduct two more simulations in which we set either the return to experience heterogeneity $\delta_i^\alpha$ or the transitory shock variances $\sigma^2_{u,t}$ to zero, but use the same innovations as in the baseline simulation. We measure inequality as the cross-sectional variance of log income. For the baseline simulation, this variance is denoted by $V_t$, and for the alternative simulations they are denoted by $V_t(\delta^\alpha=0)$ and $V_t(u=0)$.
The ratio $[V_t-V_t(\delta^\alpha=0)]/V_t$ for the ${\cal B}_{20,1988}$ sample is plotted in Figure (ref) as solid blue line. Recall that $\hat{q}^\alpha = 0.34$. Initially, return to experience heterogeneity, due to the individuals deviating from the core group, explains about 40% of income inequality in the cohort. However, the share drops quickly and after seven years $\alpha$ heterogeneity is irrelevant and income inequality is solely due to the stochastic income components. The red line corresponds to the relative effect of the transitory shock on inequality. Initially, its is about 15% but it also quickly converges to zero. Thus, after seven years, inequality is essentially solely due to the persistent idiosyncratic income shocks.
{\bf Time-path of Average Idiosyncratic Income Risk.} Figure (ref) depicts posterior means and credible intervals for $\sigma^2_{u,t}$ and $\sigma^2_{\epsilon,t}$ as a function of time. We overlay the posteriors for the two samples: $\tau=1988$ and $\tau=1991$. Years with quarters that were classified by the NBER as recessions are indicated by the grey shaded areas. It is important to note that that the cross sections for the two panels in any particular year are not identical. This explains the variation in the estimates of the year-specific standard deviations.
While the volatility spikes do not line up perfectly with the NBER recessions, there is a positive correlation. For the 1988 sample $\sigma^2_{u,t}$ spikes in the 1970, 1974, and 1982 recessions. On the other hand, it does not spike in the 1980 recession, but does have a pronounced peak during the recovery in 1984. The innovation variance $\sigma^2_{\epsilon,t}$ peaks in 1970 and 1980, but not in the other recessions. Note that individuals in our sample are always employed. Thus, we are not capturing unemployment risk, which contributes to the cyclical variation of income risk.
{\bf How Are Deviators Different?} In Figure (ref) we compare characteristics of the core group members and the deviators. The top panels show two histograms and the bottom panels show a histogram and a density estimate. In terms of the age distribution and years of education the core group and deviators are quite similar. The core group contains slightly fewer individuals below the age of 50 and not as many individuals of age 54, but if the information is aggregated to a bin width of two or three years, many of these differences would vanish. The same is true for years of education. The distribution of time series means (computed prior to projecting onto demographics) for the deviators and the core group also look quite similar, except for the presence of a few high(er) income individuals among the deviators. A similar pattern emerges for the time series standard deviation of the residual income after projecting it unto the demographics. Overall, in terms of the distribution of age, years of education, mean earnings, and standard deviation of earnings residuals the core group members and deviators look quite similar.
We now examine the predictive properties of our estimated models. We consider two types of analyses. In the first part, we use the predictive performance as a way of assessing several restricted versions of model $M_2$. For a given sample ${\cal B}_{T,\tau}$, we estimate the model specifications based on data up to period $T$ and then we forecast income in period $T+1$ for the individuals $i=1,\ldots,N$ and evaluate the accuracy of the predictions, averaging across all individuals. In the second part of this section we examine the prediction problem from the perspective of an individual who would like to forecast her/his income to make economic decisions, e.g., a consumption/savings decision.
{\bf Assessing Model Specifications Using Predictive Performance.} In Table (ref) we listed three variants of the baseline model $M_2$: $M_2$(homosk), $M_2$(RIP), and $M_2$(HIP). Our pseudo-out-of-sample approach to compare competing model specifications is different from what is typically done in the income-dynamics literature. Since a large number of authors have estimated the income models using minimum-distance techniques, e.g., by minimizing the discrepancy between sample autocovariances $\widehat{\mbox{Cov}}[y_{it},y_{i,t-h}]$ and the model-implied population analogues, it is natural to consider tests of the hypothesis $v_{\delta^\alpha}=0$. In our Bayesian setting we could in principle use posterior odds tests of $q^l=0$ versus $q^l > 0$, but because the income processes are often embedded in life-cycle models in which agents have to predict their future income to make optimal decisions, we directly focus on the pseudo-out-of-sample performance.
We estimate the model specifications for a sequence of samples ${\cal B}_{T,\tau}$, starting with $\tau=1988$ and ending with $\tau=1996$. Throughout, we consider $T=20$. The forecast evaluation approach is essentially the same as in LiuMoonSchorfheide2017, LiuMoonSchorfheide2018, and Liu2023. As evaluation criteria we consider mean-squared-error (MSE) statistics and log predictive scores (LPS). The MSEs are computed as in Step 3 of (ref) and the LPS are defined as $\frac{1}{N} \sum_{i=1}^N \ln\hat{p}(y_{iT+1}|Y_{1:N,1:T})$, where $\hat{p}(\cdot)$ is the predictive density with a particular model specification.
The results are summarized in Figure (ref). In the left panel of the figure we show relative changes in MSEs [%] defined as $100 \cdot \big(\mbox{MSE}(M_2(xxx))-\mbox{MSE}(M_2)\big)/\mbox{MSE}(M_2)$, where $M_2(xxx)$ is one of the alternative specifications in Table (ref). Thus, a negative value indicates that model $M_2(xxx)$ attains a lower MSE than the benchmark model $M_2$. The LPS differentials are defined such that a negative value is an improvement, meaning that we report $\mbox{LPS}(M_2)-\mbox{LPS}(M_2(xxx))$. To convert them into predictive log odds, they would have to be multiplied by the cross-sectional dimension $N$, which various across samples but is between 82 and 127.
We draw the following conclusions from Figure (ref). First, the fully heterogeneous $M_2$(HIP) specification has a larger MSE than the baseline model $M_2$ for eight out of nine samples. Moreover, in terms of the LPS, $M_2$ clearly dominates $M_2$(HIP). Second, there are only three samples for which $M_2$(homosk) attains a lower MSE than $M_2$. In the other six samples $M_2$ delivers more accurate point forecasts. More importantly, the density forecasting performance of $M_2$(homosk) is dismal compared to $M_2$. This confirms that the cross-sectional heteroskedasticity is an important feature of the data. Finally, $M_2$(RIP) delivers a similar point forecast performance as $M_2$ in five and a similar LPS in six out of nine samples. In the remaining samples $M_2$(RIP) slightly dominates the baseline model. This is not surprising, because as we have seen in Section (ref) regression coefficient heterogeneity is generally quite sparse and imposing a restriction that is not quite correct often favorably trades off a slight increase in bias against a variance reduction.
{\bf Forecasting From an Individual's Perspective.} We previously generated simultaneous forecasts for all units to compare the predictive ability of competing model specifications. Now we consider the forecasting problem from the perspective on a particular individual $i$. We consider three different interval forecasts. The baseline forecast of individual $i$ conditions on the information $Y_{1:N,0:T}$ and accounts for parameter uncertainty. ChamberlainHirano1999 interpret this scenario as\footnote{This setup is also considered by GuKoenker2014.} “an individual seeking advice from a financial planner. The individual provides data on his earnings history and on various personal characteristics such as age and education. The planner has access to longitudinal data sets that provide data on earnings histories and personal characteristics for samples of individuals.” (Abstract, Page 211) The second set of interval forecasts is also based on the $Y_{1:N,0:T}$ information, but we shut down the parameter uncertainty by fixing $(\delta, \theta)$ at the posterior mean estimates.
Third, we generate forecasts based on the time series information of individual $i$: $Y_{i,0:T}$. To obtain forecasts based on $Y_{i,0:T}$, Equations ((ref)) and ((ref)) of model $M_2$ are re-estimated based on data from individual $i$ only. The prior distributions are summarized in Table (ref). The hyperparameters for the priors are selected based on the prior specification for $M_2$, as outlined in Table (ref). Specifically, the prior variance of $\delta^{\alpha}_i$ and $\delta^{\rho}_i$ are set to be equal to the prior mean of $v^{\delta^\alpha}$ and $v^{\delta^\rho}$, respectively. The prior mean of $\delta^{\sigma}_{i,u}$ and $\delta^{\sigma}_{i,\epsilon}$ are aligned with the prior means of $\sigma^2_{u,t}$ and $\sigma^2_{\epsilon,t}$. Both the prior mean and variance of $s_{i0}$ are matched with the prior mean of $\mu_{s_0}$ and $v_{s_0}$.
Figure (ref) depicts fan charts to summarize interval forecasts under the three scenarios for a member of the core group, the unit with ID 74, and one of the deviators, the unit with ID 59. We show their actual income for periods $t=1,\ldots,T=20$ and then interval forecasts for the next five years. A comparison of the panels in the first and second row of the figure suggests that shutting down parameter uncertainty shortens the predictive intervals somewhat. But, overall, the uncertainty appears to be dominated by future shocks. The forecasts in the third row of the figure are obtained by conditioning on the smaller information set $Y_{i,0:T}$, which widens the intervals relative to the baseline forecasts in the first row. The reason is that the cross-section provides information about the distribution of the idiosyncratic parameters that can be used as a prior to generate sharper time series estimates for unit $i$.
To conduct a more systematic analysis of how parameter uncertainty and the information set used for generating the forecasts affect interval and density forecasts, we plot 90% interval width ratios for all individuals $i$ in the ${\cal B}_{20,1991}$ sample in Figure (ref) as a function of years of experience. The denominators of the ratios are obtained from the forecasts based on $Y_{1:N,0:T}$ information with parameter uncertainty (baseline). The numerators for the left panels are the interval widths associated with the no-parameter uncertainty forecasts, whereas the numerators for the right panels are constructed from the $Y_{i,0:T}$. As suggested by Figure (ref), the ratios on the left are below one, whereas the ratios on the right are above one.
Ratios for the core group are indicated through blue dots and values for the deviators by red dots. For the no-parameter-uncertainty-to-baseline ratios in the left panels, we expect the average ratio for the deviators to be smaller than for the core group. Because the parameter estimates for the deviators do not benefit from pooling observations across units, parameter uncertainty should play a more important role for the predictive distribution. For $h=1$ the average for the deviators is 0.91 and for the core group is 0.92. For $h=5$ the discrepancy is more pronounced as the average of the deviators drops to 0.88. However, none of these differences are statistically significant.
For the $Y_{i,0:T}$-to-$Y_{1:N,0:T}$ information ratios plotted in the right panel, we expect the core group members to benefit more from including the information from other individuals. How much more depends on the size of the deviations, captured by the $v_\delta$ parameters. In principle, this translates into higher ratios for core group members than for deviators. It turns out, this effect is not reflected in the average ratios. For $h=1$ ($h=5$) the average core group ratio is 1.4 (1.3), whereas the average deviator ratio is 1.45 (1.31). As for the ratios in the left panels, the differences are not statistically significant. One reason for the insignificance could be that the classification of units into core members and deviators is not sharp. In general there is posterior uncertainty about the memberships, which might explain the insignificance.
Uncertainty about future income plays an important role in life-cycle models because it determines the magnitude of precautionary savings. Many models assume full information rational expectations, which in the context of our analysis means that there is no uncertainty about the parameters of the income process. The only source of uncertainty is the realization of future shocks. Some authors, e.g., Guvenen2007, explicitly model the Bayesian learning about heterogeneous parameters, and account for parameter uncertainty in addition to shock uncertainty. Our estimates highlight, that the degree of uncertainty varies considerably based on the assumptions how the income forecasts are generated. The no-parameter uncertainty case could be interpreted as rational expectations. Accounting for parameter uncertainty increases the width of the 90% predictive intervals by approximately 10%. Using a sparser information set that only includes the income trajectory of individual $i$ herself, raises the width of the intervals by an additional 40%.
{\bf Idiosyncratic Parameters.} We now re-estimate model $M_2$ based on the unbalanced panel ${\cal U}$. Posterior estimates of the idiosyncratic parameters $\alpha_{i0}$, $\alpha_{i1}/10$, $\rho_i$, and $\delta_{i,\epsilon}^\sigma$ are provided in Figure (ref). In each panel, we sort the units by the posterior median estimates. As before, heteroskedasticity, i.e., full heterogeneity of $\delta_i^\sigma$ is a key feature of the panel. Unlike in the case of the balanced panels in Figure (ref), there is no evidence of sparsity in the $\delta_{i}^\alpha$ estimates in this larger panel data set. The sorted posterior median estimates of $\alpha_{i0}$ and $\alpha_{i1}/10$ do not exhibit any flat segments. The only parameter that does have a flat segment in the posterior median function is $\rho_i$, but the flat segment is not as pronounced as under the ${\cal B}_{20,1988}$ and ${\cal B}_{20,1991}$ samples. Thus, as expected, expanding the cross-sectional dimension of the panel by including more individuals, pushes the estimates toward full heterogeneity.
{\bf Posterior Distribution of $\theta$.} Parameter estimates from the baseline model and the HIP and RIP versions are reported in Table (ref). Consistent with Figure (ref), the estimates of $q^\alpha$ and $q^\rho$ are very close to one, implying full heterogeneity. The posterior median estimates of $\alpha_0$ and the return to experience $\alpha_1/10$ are with -0.21 and and .017, respectively, very close to the estimates for the balanced samples ${\cal B}_{20,1988}$ and ${\cal B}_{20,1991}$ reported in Table (ref). Due to the larger sample size, the credible intervals for those parameters are much smaller for the unbalanced sample.
An important difference between the estimates from the ${\cal U}$ and ${\cal B}_{T,\tau}$ samples manifests itself in regard to the autocorrelation of the persistent income component $\rho$ and the return-to-experience heterogeneity captured by $v_{\delta^{\alpha_1}}$. For the unbalanced sample $\hat{\rho} =0.63$, which is considerably lower than the estimates obtained for the balanced samples and the numbers generally reported in the literature. In regard to $v_{\delta^{\alpha_1}}$, recall that the estimates for the two balanced panels reported in Table (ref) were .0002 and .0001. In the unbalanced panel the estimate is an order of magnitude larger: $\hat{v}_{\delta^{\alpha_1}} = .002$. Qualitatively, the low estimate of $\rho$ is consistent with the larger return-to-experience heterogeneity, as these two parameters have found to be negatively related in the existing literature.
The low estimate of $\rho$ and the high estimate of $\hat{v}_{\delta^{\alpha_1}}$ implies that return to experience heterogeneity and transitory shocks play a more important role for income inequality within cohorts. For the ${\cal B}_{20,1988}$ sample in Figure (ref) we found that with five years of experince, $\alpha$ heterogeneity generates about 8% of the income inequality and the transitory income shock about 4%. With eight years of experience the contribution is essentially zero. Using the estimates from the unbalanced panel instead, $\alpha$ heterogeneity and transitory shock explain 40% and 10% of the income inequality within the cohort after five years of work experience and 25% and 5% after ten years.
Comparing the $v_{\delta^\sigma}$ estimates from the unbalanced panel to those from the balanced panels in Table (ref) one also notices a substantial increase in the shock variances. For instance, for ${\cal B}_{20,1988}$ the posterior median estimate of $v_{\delta^\sigma_\epsilon}$ is 4.58, whereas it equals 35.85 for the unbalanced panel. To interpret the magnitude of $\delta_{i,\epsilon}^\sigma$ (and its variance), recall that it has to be multiplied by $\sigma^2_{\epsilon,t}$ which we plotted in Figure (ref) for the balanced panels. The values for the unbalanced panel are slightly smaller. Overall, the cross-sectional averages of $\sigma^2_{\epsilon,t} \delta_{i,\epsilon}^\sigma$ are quite similar across samples.
The remaining columns of Table (ref) contain estimates for $M_2$(HIP) and $M_2$(RIP). The HIP estimates, generated conditional on $q^\alpha=q^\rho=1$, are essentially equal to the baseline estimates, which feature $\hat{q}^\alpha =0.99$ and $\hat{q}^\rho = 0.89$. The RIP estimates feature a substantially larger $\hat{\rho} = 0.97$. Moreover, the estimated return to experience is about three percentage points larger than for $M_2$ and $M_2$(HIP). In the bottom rows of Table (ref) we report results for the one-step-ahead out-of-sample predictive performance. As previously in Figure (ref), we compute LPS differentials and relative RMSE differentials (in percent) with respect to the baseline model $M_2$. Negative values indicate an improvement. The differentials between $M_2$ and $M_2$(HIP) are essentially zero. The RIP specification with the more persistent income process leads to a small improvement in terms of point (RMSE) and density (LPS) forecasts. The RMSE gain is about 1.5%. Multiplying the LPS differential by the cross-sectional sample size of approximately 800 observations at the forecast origin, we obtain log odds of about 5.6 in favor of the RIP specification. Here the parsimony of the RIP specification seems to improve the out-of-sample forecast performance if we average across all the units.
Using balanced panels of individuals with a long employment history, we found that heterogeneity with respect to returns to experience and the autocorrelation of the persistent stochastic income component is sparse, with large core groups of individuals sharing identical coefficients and a small number of deviators. For a large unbalanced panel, we found evidence for full coefficient heterogeneity, with the caveat that the specification with homogeneous coefficients performs better in terms of one-step-ahead out-of-sample forecasting. We also examined forecasts for specific individuals conditional on information from the full panel $Y_{1:N,0:T}$ accounting for shock and parameter uncertainty, the full panel ignoring parameter uncertainty, and the individual history $Y_{i,0:T}$ accounting for shock and parameter uncertainty. This comparison has an important message for the macroeconomic life-cycle model literature: subtle assumptions about the agents information set can have large effects on the income uncertainty that agents face and hence on their precautionary savings motive.
We incorporated a version of a spike and slab prior, comprising a pointmass at zero (“spike”) and a Normal distribution around zero (“slab”) into a dynamic panel data framework to model coefficient heterogeneity. Our framework nests coefficient homogeneity and full heterogeneity as a special case, but can also capture sparse heterogeneity, meaning that there is a core group of units that share the same coefficients. The remaining units deviate from the core group. There is a straightforward extension of our framework to correlated random effects and one could also treat the slab part of the $\delta_i$ distribution nonparametrically as in Liu2023. However, in many empirical applications simple parametric specifications tend to work well and we hope that the proposed framework will prove to be useful beyond the application presented in this paper.
\setcounter{page}{1}