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.
97,271 characters · 19 sections · 70 citation commands
Density Forecasts in Panel Data Models: A Semiparametric Bayesian Perspective
\thispagestyle{empty}
\setcounter{page}{1}
Panel data, such as a collection of firms or households observed repeatedly for a number of periods, are widely used in empirical studies. It can also be useful for forecasting individuals' future outcomes, which is interesting and important in many applications, for example, PSID for income dynamics Hirano2002,gu2014unobserved and bank balance sheet data for bank stress tests LiuMoonSchorfheide2015. This paper constructs individual-specific density forecasts using a dynamic linear panel data model with common and heterogeneous coefficients as well as cross-sectional heteroskedasticity.
In this paper, I consider young firm dynamics as the empirical application. For illustrative purposes, consider a simple dynamic panel data model as the baseline setup:
where $i=1,\cdots,N$, and $t=1,\cdots,T+1$. $y_{it}$ is the observed firm performance such as log employment, $\lambda_{i}$ is the unobserved skill of an individual firm, and $u_{it}$ is an i.i.d.\ shock. Skill is independent of the shock, and the shock is independent across firms and times. $\beta$ and $\sigma^{2}$ are common across firms, where $\beta$ represents the persistence of the dynamic pattern and $\sigma^{2}$ gives the size of the shocks. Based on the observed panel from period $0$ to period $T$, I am interested in forecasting the future performance of any specific firm in period $T+1$.
The panel considered in this paper features a large cross-sectional dimension $N$ but short time series $T$. For instance, the number of observations for each young firm is restricted by its age. Good estimates of\textcolor{red}{ }the unobserved skill $\lambda_{i}$ facilitate good forecasts of $y_{i,T+1}$. Because of the short $T$, traditional methods have difficulty in disentangling the unobserved skill $\lambda_{i}$ from the shock $u_{it}$, which contaminates the estimates of $\lambda_{i}$, even if $N$ goes to infinity.
To tackle this problem, I assume that $\lambda_{i}$ is drawn from an underlying skill distribution $f$ and estimate this distribution by combining information from the whole panel. In terms of modeling $f$, the parametric Gaussian density misses many features in real-world data, such as asymmetry, heavy tails, and multiple peaks. For example, as good ideas are scarce, the skill distribution of young firms may be highly skewed. This calls for a flexible modeling of $f$, and here I estimate $f$ via a nonparametric Bayesian approach where the prior is constructed from a mixture model and allows for correlation between $\lambda_{i}$ and the initial condition $y_{i0}$ (i.e.\ a correlated random effects model).
Conditional on $f$, we can treat it as a prior distribution and combine it with firm-specific data to obtain the firm-specific posterior via Bayes' theorem. In a special case where the common parameters are $\left(\beta,\sigma^{2}\right)=\left(0,1\right)$, the firm-specific posterior is
This firm-specific posterior helps better infer the firm-specific unobserved skill $\lambda_{i}$ and better forecast the firm-specific future performance, thanks to the estimated underlying distribution $f$ that integrates the information from the whole panel in an efficient and flexible way. This is only an intuitive explanation of why the skill distribution $f$ is crucial. In the actual implementation, the correlated random effect distribution $f$, common parameters $\left(\beta,\sigma^{2}\right)$, and firm-specific skill $\lambda_{i}$ are all inferred simultaneously.
It is natural to construct density forecasts based on the firm-specific posterior. In general, forecasting can be done in a point, interval, or density manner, with density forecasts giving the richest insight into future outcomes. By definition, a density forecast provides a predictive distribution of firm $i$'s future performance and summarizes all sources of uncertainties; hence, it is preferable in the context of young firm dynamics and other applications with large uncertainties and nonstandard distributions. In particular, for the baseline model in ((ref)), the density forecasts reflect uncertainties arising from the future shock $u_{i,T+1}$, unobserved individual heterogeneity $\lambda_{i}$, and estimation uncertainty of common parameters $\left(\beta,\sigma^{2}\right)$ and of skill distribution $f$. Moreover, once density forecasts are obtained, one can easily recover point and interval forecasts.
The contributions of this paper are threefold. First, I establish the theoretical properties of the proposed predictor when the cross-sectional dimension $N$ tends to infinity. To begin, I provide conditions for identifying the common parameters and the distribution of the individual heterogeneity in both cross-sectional homoskedastic and heteroskedastic models. Then, I prove that the proposed estimator achieves posterior consistency in cross-sectional homoskedastic cases. Compared with previous literature on posterior consistency in density estimation problems, there are several challenges in the panel data framework: (1) a deconvolution problem disentangling unobserved individual effects and shocks, (2) an unknown common shock size in cross-sectional homoskedastic cases, (3) strictly exogenous and predetermined variables (including lagged dependent variables) as covariates, and (4) correlated random coefficients addressed by flexible conditional density estimation. Based on the posterior consistency of the estimates, the discrepancy between the proposed density predictor and the oracle is arbitrarily small asymptotically. The oracle predictor is an (infeasible) benchmark defined as the individual-specific posterior predictive distribution, assuming known common parameters and a known distribution of the heterogeneous parameters.
Second, I develop a posterior sampling algorithm specifically addressing nonparametric density estimation of the unobserved individual effects. For a random coefficients model, which is a special case where the individual effects are independent of the conditioning variables, the $f$ part becomes an unconditional density estimation problem. I adopt a Dirichlet Process Mixture (DPM) prior for $f$ and construct a posterior sampler building on the blocked Gibbs sampler proposed by IshwaranJames2001,IshwaranJames2002. For a correlated random coefficients model, I further adapt the proposed algorithm to the much harder conditional density estimation problem using a probit stick-breaking process prior suggested by PatiDunsonTokdar2013.
Third, Monte Carlo simulations demonstrate improvement in density forecasts relative to alternative predictors with various parametric priors on $f$, evaluated by the log predictive score. An application to young firm dynamics also shows that the proposed predictor provides more accurate density predictions. The better forecasting performance is largely due to three key features (in order of importance): the nonparametric Bayesian prior, cross-sectional heteroskedasticity, and correlated random coefficients. The estimated model also helps shed light on the latent heterogeneity structure of firm-specific coefficients and cross-sectional heteroskedasticity, as well as whether and how the unobserved heterogeneity depends on the initial condition of the firms.
Moreover, the proposed method is applicable beyond forecasting. Here estimating heterogeneous parameters is important because we want to generate good individual-specific forecasts, but in other cases, the heterogeneous parameters themselves could be the objects of interest. For example, the technique developed here can be adapted to infer individual-specific treatment effects.
\paragraph{Related Literature }
First, this paper contributes to the literature on individual forecasts in a panel data setup, and is closely related to LiuMoonSchorfheide2015 and GuKoenker2016,gu2014unobserved. LiuMoonSchorfheide2015 focus on point forecasts. They utilize the idea of Tweedie's formula to steer away from the complicated deconvolution problem in estimating $\lambda_{i}$ and establish the ratio optimality of point forecasts. Unfortunately, the Tweedie shortcut is not applicable to the inference of the underlying $\lambda_{i}$ distribution and therefore not suitable for density forecasts. In addition, this paper addresses cross-sectional heteroskedasticity where $\sigma_{i}^{2}$ is an unobserved random quantity, while LiuMoonSchorfheide2015 incorporate cross-sectional and time-varying heteroskedasticity via a deterministic function of observed conditioning variables.
gu2014unobserved address the density estimation problem, but with a different method. This paper infers the underlying $\lambda_{i}$ distribution via a full Bayesian approach (i.e.\ adopting a prior on the $\lambda_{i}$ distribution and updating the prior belief by the observed data), whereas they employ an empirical Bayes approach (i.e.\ choosing the $\lambda_{i}$ distribution by maximizing the marginal likelihood of data). In principle, the full Bayesian approach is preferable for density forecasts, as it captures all sources of uncertainties, including estimation uncertainty of the underlying $\lambda_{i}$ distribution, which has been omitted by the empirical Bayes approach. In addition, this paper features correlated random coefficients allowing the cross-sectional heterogeneity to interact with the initial conditions, whereas gu2014unobserved focus on random effects models without this interaction.
In their recent paper, GuKoenker2016 also compare their method with an alternative semiparametric Bayesian estimator featuring a Dirichlet Process (DP) prior under a set of fixed scale parameters. There are two major differences between their DP setup and the DPM prior used in this paper. First, the DPM prior provides continuous individual effect distributions, which could be the case in many empirical setups. Second, unlike their set of fixed scale parameters, this paper incorporates a hyperprior for the scale parameter and updates it via the observed data, hence let the data choose the complexity of the mixture approximation, which can essentially be viewed as an “automatic” model selection.
Earlier works on full Bayesian analyses with parametric priors on $\lambda_{i}$ can be found in Lancaster01072002 (orthogonal reparametrization and a flat prior); Chamberlain1999, chib1999mcmc, and sims2000using (Gaussian prior); and chib2008panel (student-t and finite mixture priors). There have also been empirical works on the DPM model with panel data, but they mostly focus on empirical studies rather than theoretical analyses. For example, Hirano2002 and jensen2015mutual use linear panel models with setups different from this paper. Hirano2002 considers flexibility in the $u_{it}$ distribution instead of the $\lambda_{i}$ distribution. jensen2015mutual assume random effects instead of correlated random effects. burda2013panel and 10.2307/j.ctt5hhrfp use a panel probit model and a panel logit model, respectively.
In the frequentist literature, LI1998139, delaigle2008deconvolution, evdokimov2010identification, and hu2017econometrics, among others, have studied a similar deconvolution problem and estimated the $\lambda_{i}$ distribution. Also see ECTJ:ECTJ12068 for a review of frequentist applications of mixture models. However, the frequentist approach misses estimation uncertainty, which matters in density forecasts, as mentioned previously.
Second, this paper also relates to the literature on nonparametric Bayesian methods in density estimation problems ghosh2003bayesian,Hjort2010,ghosal2017fundamentals. In particular, for unconditional density estimation, a recent paper by Canale2017 relaxed the tail conditions to accommodate multivariate location-scale mixtures. For conditional density estimation, the mixing probabilities can be characterized by a multinomial choice model norets2010approximation,norets2012bayesian, a kernel stick-breaking process ECT:9258097,pelenis2014bayesian,Norets2017, or a probit stick-breaking process PatiDunsonTokdar2013. I adopt the PatiDunsonTokdar2013 approach and establish posterior consistency for a multivariate conditional density estimator featuring infinite location-scale mixtures with a probit stick-breaking process.
To account for deconvolution, I construct an inversion inequality that links the convergence of the distribution of observables to the convergence of the distribution of the unobserved individual heterogeneity. The latter is in the Wasserstein metric, which is useful in handling deconvolution problems as found in the recent literature. For example, nguyen2013convergence considers the unobserved distribution on a discrete support, and su2020nonparametric flexibly model a symmetric unimodal unobserved distribution using a mixture of symmetric uniforms where the bounds are drawn from a Dirichlet process location-mixture of Gammas. Their setups, however, differ from the current framework, which calls for a new inversion inequality developed in this paper. Then, I further take into account the dynamic panel data structure, as well as obtain the convergence of the proposed predictor to the oracle predictor.
Last but not least, the empirical application in this paper also relates to the young firm dynamics literature. AkcigitKerr2010 document that R&D intensive firms grow faster, especially for smaller firms. robb2014role examine the role of R&D in capital structure and performance of young firms. The empirical analysis of this paper builds on these findings. Besides more accurate density forecasts, I also obtain the latent heterogeneity structure of firm-specific coefficients and cross-sectional heteroskedasticity.
The rest of the paper is organized as follows. Section (ref) specifies the general panel data model, density forecasts, and nonparametric Bayesian priors; Section (ref) establishes the posterior consistency of the estimates and the convergence of the density forecasts to the oracle; Section (ref) conducts Monte Carlo simulations; Section (ref) presents the empirical application to young firm dynamics; and Section (ref) concludes. Notations, proofs, algorithms, and additional results are in the Appendix.
The general panel data model with (correlated) random coefficients and potential cross-sectional heteroskedasticity can be specified as
where $i=1,\cdots,N$, and $t=1,\cdots,T+h$. Similar to the baseline setup in ((ref)), $y_{it}$ is the observed individual outcome, such as young firm performance. The main goal of this paper is to estimate the model using the sample from period $0$ to period $T$ and forecast the future distribution of $y_{i,T+h}$ for any individual $i$. In the remainder of the paper, I focus on the case where $h=1$ (i.e.\ one-period-ahead forecasts) for notational simplicity, and the discussion can be extended to multi-period-ahead forecasts via either a direct or an iterated approach marcellino2006comparison.
$w_{i,t-1}$ is a vector of observed covariates that have heterogeneous effects on the outcomes, with $\lambda_{i}$ being the unobserved heterogeneous coefficients. $w_{i,t-1}$ is strictly exogenous and captures key sources of individual heterogeneity. If $w_{i,t-1}=1$, $\lambda_{i}$ is reduced to an individual-specific intercept, e.g.\ firm $i$'s skill level in the baseline model ((ref)). More generally, $w_{i,t-1}$ can contain individual-specific variables (e.g.\ firm-specific R&D) and aggregate variables (e.g.\ a recession dummy). I focus on the former case below for notational simplicity. In the latter case, all theoretical analyses would be further conditioned on the aggregate observations.
$x_{i,t-1}$ is a vector of observed covariates that have homogeneous effects on the outcomes, and $\beta$ is the corresponding vector of common parameters. I decompose $x_{i,t-1}=\left[x_{i,t-1}^{O\prime},x_{i,t-1}^{P\prime}\right]^{\prime}$, where $x_{i,t-1}^{O}$ is strictly exogenous and $x_{i,t-1}^{P}$ is predetermined. One example of $x_{i,t-1}^{P}$ is the lagged outcome $y_{i,t-1}$ capturing the persistence. Both $x_{i,t-1}^{O}$ and $x_{i,t-1}^{P}$ can include other control variables, such as firm characteristics and general economic conditions. Let $x_{i,t-1}^{P*}$ denote the subgroup of $x_{i,t-1}^{P}$ excluding lagged outcomes, then $x_{i,t-1}=\left[x_{i,t-1}^{O\prime},x_{i,t-1}^{P*\prime},y_{i,t-1}\right]^{\prime}$ with $\beta=\left[\beta^{O\prime},\beta^{P*\prime},\rho\right]^{\prime}$. Here, the distinction between homogeneous effects $\beta^{\prime}x_{i,t-1}$ and heterogeneous effects $\lambda_{i}^{\prime}w_{i,t-1}$ helps model the key latent heterogeneities while avoiding the curse of dimensionality. Combining information from the covariates, the conditioning set at period $t$ is defined as $c_{i,t-1}=\left(x_{i,0:t-1}^{P},x_{i,0:T}^{O},w_{i,0:T}\right).$ We further define $D=\left(\left\{ D_{i}\right\} _{i=1}^{N}\right)$, where $D_{i}=c_{iT}$, as the data used for estimation and the conditioning set for posterior inference.
$u_{it}$ is an individual-time-specific shock characterized by zero mean and potential cross-sectional heteroskedasticity $\sigma_{i}^{2}$, with cross-sectional homoskedasticity being a special case where $\sigma_{i}^{2}=\sigma^{2}$. In a unified framework, I denote the common parameters by $\vartheta$, the individual heterogeneity by $h_{i}$, and the underlying distribution of $h_{i}$ by $f$. For instance, $\vartheta=\beta,\;h_{i}=\left(\lambda_{i},\sigma_{i}^{2}\right)$ under cross-sectional heteroskedasticity. In many empirical applications, such as the young firm example, the size of risk may vary over the cross-section, so cross-sectional heteroskedasticity could contribute to better density forecasts.
As stressed in the motivation, the underlying distribution of individual effects is the key to better density forecasts. In the literature, there are usually two types of assumptions on this distribution. One is the random coefficients model, where the individual effects $h_{i}$ are independent of the conditioning variables $c_{i0}=\left(x_{i0}^{P},x_{i,0:T}^{O},w_{i,0:T}\right)$. The other is the correlated random coefficients model, where $h_{i}$ and $c_{i0}$ could be correlated. This paper considers both models while focusing on the latter---although the former is more parsimonious and easier to implement, the latter is more realistic for young firm dynamics as well as many other empirical setups. In practice, it is more feasible to only take into account a subset of $c_{i0}$ or a function of $c_{i0}$ that is relevant for the specific study.
This subsection formally defines the infeasible optimal oracle predictor and the feasible semiparametric Bayesian predictor proposed in this paper. Both definitions rely on the conditional predictor,
which provides the density forecasts of $y_{i,T+1}$ conditional on the common parameters $\vartheta$, underlying distribution $f$, and individual $i$'s data $D_{i}$. The first term $p\left(\left.y\right|h_{i},\vartheta,w_{iT},x_{iT}\right)$ captures individual $i$'s uncertainty due to the future shock $u_{i,T+1}$. The second term \[ p\left(h_{i}\left|\vartheta,f,D_{i}\right.\right)=\frac{\prod_{t=1}^{T}p\left(\left.y_{it}\right|h_{i},\vartheta,w_{i,t-1},x_{i,t-1}\right)f\left(h_{i}\left|c_{i0}\right.\right)}{\int\prod_{t=1}^{T}p\left(\left.y_{it}\right|h_{i},\vartheta,w_{i,t-1},x_{i,t-1}\right)f\left(h_{i}\left|c_{i0}\right.\right)dh_{i}} \] is the individual-specific posterior. It characterizes individual $i$'s uncertainty due to unobserved individual heterogeneity that arises from insufficient time-series information to infer individual $h_{i}$. The common distribution $f$ helps regulate this source of uncertainty and hence contributes to individual $i$'s density forecasts.
The infeasible oracle predictor is defined as if we knew all the elements that can be consistently estimated. Specifically, the oracle knows the common parameters $\vartheta_{0}$ and the underlying distribution $f_{0}$, but not the individual effects $h_{i}$. Then, the oracle predictor is formulated by plugging the true values $\left(\vartheta_{0},f_{0}\right)$ into the conditional predictor in ((ref)),
In practice, $\left(\vartheta,f\right)$ are unknown and need to be estimated, thus introducing another source of uncertainty. For the common parameters $\vartheta$, I adopt a conjugate prior (e.g.\ mulitvariate normal for cross-sectional heteroskedastic cases) in order to stay close to the linear regression framework. For the distribution of individual heterogeneity $f$, I resort to the nonparametric Bayesian prior (specified in the next subsection) to flexibly model this underlying distribution, which could better approximate the true distribution $f_{0}$, and the resulting feasible predictor would be close to the oracle. Then, I update the prior belief using the observations from the whole panel and obtain the posterior. The semiparametric Bayesian predictor is constructed by integrating the conditional predictor over the posterior distribution of $\left(\vartheta,f\right)$,
The conditional predictor reflects uncertainties due to future shocks and unobserved individual heterogeneity, whereas the posterior of $\left(\vartheta,f\right)$ captures estimation uncertainty. Note that the inference of $\left(\vartheta,f\right)$ combines information from the whole panel. Once conditioned on $\left(\vartheta,f\right)$, we have that individuals' outcomes are independent across $i$ and that only individual $i$'s data are further needed for its density forecasts.
A prior on the distribution $f$ can be viewed as a distribution over a set of distributions. Among other options, I formulate the nonparametric Bayesian prior using mixture models, because mixture models can effectively approximate a general class of distributions while being relatively easy to implement. The specific functional form depends on whether $f$ is characterized by a random coefficients model or a correlated random coefficients model.
In cross-sectional heteroskedastic cases, I incorporate another flexible prior on the distribution of $\sigma_{i}^{2}$. Define $l_{i}=\log\frac{\bar{\sigma}^{2}\left(\sigma_{i}^{2}-\underline{\sigma}^{2}\right)}{\bar{\sigma}^{2}-\sigma_{i}^{2}},$ where $\underline{\sigma}^{2}$ ($\bar{\sigma}^{2}$) is some small (large) positive number. This transformation ensures that the support of $f_{\sigma^{2}}$ is bounded by $\left[\underline{\sigma}^{2},\bar{\sigma}^{2}\right]$ for numerical stability, whereas the support of $l_{i}$ is unbounded so a similar prior structures can be applied to both $\lambda_{i}$ and $l_{i}$. We assume $\lambda_{i}$ and $\sigma_{i}^{2}$ are conditionally independent conditioning on $c_{i0}$, so their mixture structures can be modeled separately. For a concise exposition, I define a generic variable $z$ that can represent either $\lambda$ or $l$, and include $z$ in the subscript as an indicator. When there is no confusion, $z$ and $i$ in the subscript are suppressed.
In the random coefficients model, the individual heterogeneity $z_{i}\left(=\lambda_{i}\text{ or }l_{i}\right)$ is assumed to be independent of the conditioning variables $c_{i0}$, so the inference of the $f$ part can be considered as an unconditional density estimation problem, and then the DPM prior is a typical choice in the nonparametric Bayesian literature. With component label $k$, component probability $p_{k}$, and component parameters $\left(\mu_{k},\Omega_{k}\right)$, one draw from the DPM prior can be written as an infinite location-scale mixture of normals,
Different draws from the DPM prior are characterized by different combinations of $\left\{ p_{k},\mu_{k},\Omega_{k}\right\} $, which lead to different shapes of $f$. This is why the DPM prior is flexible enough to approximate a wide range of continous distributions. The component parameters $\left(\mu_{k},\Omega_{k}\right)$ are drawn from the base distribution $G_{0}$, which is chosen to be a conjugate multivariate-normal-inverse-Wishart distribution, or a normal-inverse-gamma distribution for scalar $z_{i}$. The component probability $p_{k}$ is constructed via a stick-breaking process governed by the scale parameter $\alpha$.
The scale parameter $\alpha$ controls the number of unique components in the mixture density and thus determines the flexibility of the mixture density. One advantage of the nonparametric Bayesian framework is its ability to flexibly elicit the tuning parameter, such as $\alpha$, from the data. Namely, we can set up a relatively flexible hyperprior for $\alpha\sim\text{Ga}\left(a_{\alpha,0},b_{\alpha,0}\right),$ and update it based on the observations, which “automatically” chooses the complexity of the mixture structure.
To accommodate the correlated random coefficients model where the individual heterogeneity $z_{i}\left(=\lambda_{i}\text{ or }l_{i}\right)$ can be correlated with the conditioning variables $c_{i0}$, it is necessary to consider a nonparametric Bayesian prior that is compatible with the much harder conditional density estimation problem. One issue is associated with the uncountable collection of conditional densities, and PatiDunsonTokdar2013 circumvent it by linking the properties of the conditional density to the corresponding ones of the joint density without explicitly modeling the marginal density of $c_{i0}$. As suggested in PatiDunsonTokdar2013, I utilize the Mixtures of Gaussian Linear Regressions (MGLR\textsubscript{x}) prior, a generalization of the Gaussian-mixture prior for conditional density estimation, and extend it to the multivariate setup. Conditioning on $c_{i0}$,
Similar to the DPM prior, the component parameters can be directly drawn from the base distribution, $\left(\mu_{k},\Omega_{k}\right)\sim G_{0}.$ $G_{0}$ is again specified as a conjugate matricvariate-normal-inverse-Wishart form (or a multivariate-normal-inverse-gamma distribution for scalar $z_{i}$). Now the mixture probabilities are characterized by a probit stick-breaking process
where stochastic function $\zeta_{k}$ is drawn from Gaussian process $\zeta_{k}\sim GP\left(0,V_{k}\right)$ for $k=1,2,\cdots$. rodriguez2011nonparametric demonstrate the flexibility and computational simplicity of the probit stick-breaking prior.
This setup has three key features: component means are linear in $c_{i0}$; component covariances are independent of $c_{i0}$; and mixture probabilities are flexible functions of $c_{i0}$. This framework is relatively parsimonious for finite sample implementation and, at the same time, general enough to accommodate a broad class of conditional distributions. Intuitively, it is similar to approximating the conditional density via Bayes' theorem but does not explicitly model the distribution of the conditioning variables $c_{i0}$. The infinite mixture structure and flexible mixture probabilities could absorb dependency on $c_{i0}$, so we would not need further dependency of component means and covariances on $c_{i0}$ beyond the MGLR\textsubscript{x} specification (see details in the Appendix).
In general, it is desirable to ensure that the prior belief does not dominate the posterior inference asymptotically. For Bayesians with different prior beliefs, the asymptotic properties ensure that they will eventually agree on similar predictive distributions BlackwellDubins1962,DiaconisFreedman1986. For frequentists, the asymptotic properties can be viewed as a frequentist justification for the Bayesian method---as the sample size increases, the updated posterior recovers the unknown true data generating process (DGP). Also, the conditions for posterior consistency provide guidance in choosing better-behaved priors.
In the context of infinite dimensional analysis such as density estimation, posterior consistency cannot be taken as given---the null set for the prior can be topologically large, and hence the true model can fall beyond the scope of the prior Freedman1963,Freedman1965. Therefore, it is crucial to find reasonable conditions on the joint behavior of the prior and the true density to establish the posterior consistency result.
Although identification may not be necessary to ensure the convergence of the density forecasts to the oracle predictor, identification is essential to ensure the posterior consistency of the estimates so that the proposed method could be general to problems beyond forecasting, e.g.\ heterogeneous treatment effect. Here, I present the identification result in terms of the correlated random coefficients model with cross-sectional heteroskedasticity, where random coefficients and cross-sectional homoskedasticity can be viewed as special cases and will be discussed in Remark (ref).
Despite the conditional independence in condition 1-d, $\lambda_{i}$ and $\sigma_{i}^{2}$ can potentially relate to each other through $c_{i0}$. The setup could be further extended, such as relaxing the conditional independence between $\lambda_{i}$ and $\sigma_{i}^{2}$ and allowing for more general $v_{it}$ distributions (discussed in the Appendix).
The argument is similar to ArellanoBover1995 and Arellano2012, except for the treatment of cross-sectional heteroskedasticity---here $\sigma_{i}^{2}$ is an unobserved random quantity. First, the identification of common parameters $\beta$ in panel data models is standard in the literature Baltagi1995,ArellanoHonore2001,Arellano2003,Hsiao2014. For example, the rank condition helps identify $\beta$ via orthogonal forward differencing. Second, as $\lambda_{i}$ is additively separable from the shocks, I follow the standard proof based on characteristic functions to identify $f_{\lambda}$. Finally, note that unlike $\lambda_{i}$, $\sigma_{i}^{2}$ interacts with the shocks in a multiplicative way. The Fourier transform is not suitable for disentangling products of random variables, so I resort to the Mellin transform GalambosSimonelli2004 to obtain the identification of $f_{\sigma^{2}}$.
Most of the previous nonparametric Bayesian literature focuses on density estimation problems (see Related Literature) without deconvolution and dynamic panel data structures. In this subsection, I first provide general sufficient conditions that ensure posterior consistency of the estimated common parameters $\vartheta$ and the estimated (conditional) distribution of individual effects $f$ in a general semiparametric setup, and then I specify and verify these conditions in cases of (correlated) random coefficients models.
\paragraph{General Semiparametric Model.}
Let $\varTheta$ be the space of the common parameters $\vartheta$, $\mathcal{F}$ be a set of the underlying distributions $f$ with finite second moments, $\Pi\left(\cdot,\cdot\right)$ be a joint prior on $\mathit{\Theta}\times\mathcal{F}$ with marginal priors being $\Pi_{\vartheta}\left(\cdot\right)$ and $\Pi_{f}\left(\cdot\right)$, and $\Pi\left(\cdot,\cdot|D\right)$ be the corresponding joint posterior. The individual specific likelihood takes a general “convolution” form
where $\left.D_{i}\right\backslash c_{i0}$ denotes the set difference, and $q_{0}\left(c_{i0}\right)$ is the true marginal density of $c_{i0}$.
The posterior consistency results are established with respect to the Wasserstein metric on $f$. Let $\Gamma\left(f_{1},f_{2}\right)$ be the collection of all joint measures with marginals $f_{1}$ and $f_{2}$. We define the second Wasserstein distance, $W_{2}\left(f_{1},f_{2}\right)=\left(\inf_{\gamma\in\Gamma\left(f_{1},f_{2}\right)}\int\left\Vert h_{1}-h_{2}\right\Vert _{2}^{2}d\gamma\left(h_{1},h_{2}\right)\right)^{1/2}$. Note that convergence in the $W_{2}$ metric is equivalent to weak convergence plus convergence of the second moment santambrogio2015optimal.
When $f$ is a conditional distribution, it is helpful to link the properties of the conditional density to the corresponding joint density $f\left(h,c_{0}\right)=f\left(h|c_{0}\right)q_{0}\left(c_{0}\right)$ without explicitly modeling $q_{0}$, which circumvents the difficulty associated with an uncountable set of conditional densities PatiDunsonTokdar2013. Note that $q_{0}$ is only for theoretical derivation, and there is no need to estimate it in practice.
Then, the posterior achieves consistency at $\left(\vartheta_{0},f_{0}\right)$, i.e. for all $\epsilon,\delta>0$, as $N\rightarrow\infty$, \[ \Pi\left(\left.\left(\vartheta,f\right):\;\left\Vert \vartheta-\vartheta_{0}\right\Vert _{2}<\delta,\;W_{2}\left(f,f_{0}\right)<\epsilon\right|D\right)\rightarrow1, \] in probability with respect to the true DGP.
Intuitively, let $\Theta_{\delta}^{c}=\left\{ \left\Vert \vartheta-\vartheta_{0}\right\Vert _{2}\ge\delta\right\} $, $\mathcal{F}_{\epsilon}^{c}=\left\{ W_{2}\left(f,f_{0}\right)\ge\epsilon\right\} $, and the likelihood ratio $R_{N}\left(D,\vartheta,f\right)=\prod_{i=1}^{N}\frac{g\left(\left.D_{i}\right|\vartheta,f\right)}{g\left(\left.D_{i}\right|\vartheta_{0},f_{0}\right)}$, the posterior probability of the alternative region can be decomposed as
and we want to show that the whole expression tends to zero as $N$ goes to infinity. First, for the denominator, the KL property (condition 1-a) implies that the prior puts positive weight around neighborhoods of the true DGP, so the likelihood ratio integrated over the whole space is large enough. Second, the exponentially consistent sequence of tests (condition 2) takes an infimum over the alternative region $\Theta_{\delta}^{c}\times\mathcal{F}$, so it ensures that the first term in the numerator is arbitrarily small. Third, the sieve property on $f$ (condition 3) ensures that the sieve expands to the alternative region and puts an asymptotic upper bound on the number of balls that cover the sieve. As the likelihood ratio is small in each covering ball, the integration over the alternative region is still sufficiently small Canale2017.
When $g$ is observed instead of $f$, we need to further address convolution and common parameters. In terms of convolution, it preserves the $L_{1}$-norm as well as the number of balls that cover the sieve. Moreover, the inversion inequality in condition 1-c helps identify the underlying $f$ based on the observed $g$. I use the Wasserstein metric on $f$ because there are technical difficulties in establishing a similar inversion inequality in the $L_{1}$-norm, whereas recent literature found that the Wasserstein metric circumvents the issue nguyen2013convergence,su2020nonparametric. We can extend both condition 1-c and the posterior consistency result to the $W_{p}$ metric with $p\ge1$. In terms of the common parameters, when $\vartheta$ is close to $\vartheta_{0}$ but $f$ is far from $f_{0}$, condition 1-b makes sure that the deviation generated from $\vartheta$ is small enough so that it cannot offset the difference in $f$. Therefore, conditions 1-b,c and 3 together guarantee that the data are informative enough to differentiate the true distribution from the alternatives, so the second term of the numerator can be arbitrarily small as well.
Note that the estimated individual effects $h_{i}$ are not consistent because information is accumulated only along the cross-sectional dimension but not along the time dimension. Also, the result only guarantees pointwise convergence in the space of the distributions. For uniform forecasting performance in dynamic panel data models, see LiuMoonSchorfheide2015, which considers an empirical Bayes setup with a nonparametric kernel estimate of the marginal distribution of data.
\paragraph{Random Coefficients Model.}
In this case, $f$ is an unconditional distribution. Here I focus on the cross-sectional homoskedastic case due to the difficulty in constructing a suitable mollifier in the cross-sectional heteroskedastic setup, which is left for future research. Then, the space for common parameters $\vartheta=\left(\beta,\sigma^{2}\right)$ is $\Theta=\mathbb{R}^{d_{x}}\times\left[\underline{\sigma}^{2},\;\bar{\sigma}^{2}\right].$ Let $\mathbb{E}_{f}\left[\mathfrak{g}\left(\lambda\right)\right]=\int\mathfrak{g}\left(\lambda\right)f\left(\lambda\right)d\lambda$ for a generic function $\mathfrak{g}\left(\lambda\right)$. To ensure condition 1 in Theorem (ref), we consider space $\mathcal{F}=\left\{ f:\;\mathbb{E}_{f}\left\Vert \lambda\right\Vert _{2}^{2\left(1+\eta\right)}\le M\right\} $, for some large $M>0$, and $\eta$ is defined in Assumption (ref)(1-e) below.
The conditions on $w_{i,0:T-1}$ help obtain an upper bound on the $W_{2}$-distance between $f$ and its convolution with a mollifier and hence ensure Theorem (ref)(1-c). Both conditions can be relaxed to “almost everywhere” with slight adjustments in the proofs. The moment conditions on $x_{i,t-1}$ ensure that the GMM estimates of the common parameters are asymptotically normal, so the exponentially consistent sequence of tests in Theorem (ref)(2) can be constructed accordingly. All three conditions also prevent a slight difference in $\beta$ from obscuring the difference in $f$, and are essential to Theorem (ref)(1-a,b).
First, condition 1 ensures that the true distribution $f_{0}$ is well-behaved, and a multivariate-normal-inverse-Wishart $G_{0}$ in condition 2 guarantees that the DPM prior is general enough to contain the true distribution, so the KL property on $f$ is established. Second, according to Corollary 1 in Canale2017, condition 2 further ensures the sieve property (Theorem (ref)(3)), where $2d_{w}$ controls the tail behavior of component mean $\mu$ and $\left(2d_{w}+1\right)\left(d_{w}-1\right)$ regulates the eigenvalue structure of component variance $\Omega$.
Then, the posterior achieves consistency at $\left(\vartheta_{0},f_{0}\right)$.
\paragraph{Correlated Random Coefficients Model.}
$f$ is now a conditional distribution, so the following discussion is based on the $q_{0}$-induced measure. Let $\mathcal{C}$ be the support of the conditioning variables, and $\mathcal{F}^{*}$ be a subset of conditional distributions such that mapping $c_{0}\mapsto f\left(\cdot|c_{0}\right)$ is a continous function from $\mathcal{C}$ to the space of Lebesgue integrable functions on $\mathbb{R}^{d_{w}}$. Similar to the above discussion on random coefficients models, I focus on the cross-sectional homoskedastic case and consider space $\mathcal{F}=\left\{ f:\;\left\{ \mathbb{E}_{f,q_{0}}\left\Vert \lambda\right\Vert _{2}^{2\left(1+\eta\right)}\le M\right\} \cap\mathcal{F}^{*}\right\} $, where $\mathbb{E}_{f,q_{0}}\left[\mathfrak{g}\left(\lambda,c_{0}\right)\right]=\int\mathfrak{g}\left(\lambda,c_{0}\right)f\left(\lambda|c_{0}\right)q_{0}\left(c_{0}\right)d\lambda dc_{0}$ for a generic function $\mathfrak{g}\left(\lambda,c_{0}\right)$. $M$ is some large positive constant, and $\eta$ is defined in Assumption (ref)(1-e) below.
The compactness ensures uniform convergence on $\mathcal{C}$ in the proof of the KL property. It is stronger than the $\mathcal{C}$ part in Assumption (ref)(1,3) for random coefficients models.
These conditions build on PatiDunsonTokdar2013 for posterior consistency under the conditional density topology and further extend it to multivariate conditional density estimation with infinite location-scale mixtures. The conditions on $f_{0}$ and $G_{0}$ can be viewed as conditional density analogs of the conditions in Assumption (ref). In terms of the stick-breaking process, the variability of $p_{k}\left(c_{0}\right)$ due to $c_{0}$ decreases with component index $k$ according to condition 3-b, so the first several “sticks” would be able to capture a large fraction of the dependence of $\lambda$ on $c_{0}$. Moreover, the tail of $A_{k}$ cannot be too fat according to condition 3-c.
Then, the posterior achieves consistency at $\left(\vartheta_{0},f_{0}\right)$.
Based on posterior consistency, we can bound the discrepancy between the proposed predictor and the oracle by estimation uncertainties in $\vartheta$ and $f$, and then show the asymptotic convergence of the density forecasts to the oracle forecast. Theorem (ref) in the Appendix established the convergence result in the general semiparametric setup, and the following theorem focuses on the (correlated) random coefficients models considered in the paper.
The asymptotic convergence of aggregate-level density forecasts can then be derived by summing individual-specific forecasts over different subcategories.
This section conducts two sets of Monte Carlo simulation experiments: the baseline setup with random effects, and the general setup with correlated random coefficients and cross-sectional heteroskedasticity. The main text focuses on density forecast results, whereas point forecast results are deferred to the Appendix.
The accuracy of the density forecasts is measured by the log predictive score (LPS) as suggested in Geweke2010, $LPS=\frac{1}{N}\sum_{i}\log\hat{p}\left(y_{i,T+1}|D\right),$ where $y_{i,T+1}$ is the realization at $T+1$, and $\hat{p}\left(y_{i,T+1}|D\right)$ represents the predictive likelihood with respect to the estimated model conditional on the observed data $D$. $\exp\left(LPS_{A}-LPS_{B}\right)$ gives the odds of future realizations based on predictor A versus predictor B. I performed a test combining AmisanoGiacomini2007 (for the LPS) and timmermann2019comparing (for panel data, see their Section 2.6 on general loss functions) to examine the significance in the LPS difference.
Different predictors can be interpreted as different priors on the distribution of $\lambda_{i}$. As these priors are distributions over distributions, Figure (ref) plots two draws from each prior. The homogeneous prior (Homog) implies an extreme kind of pooling, which assumes that all firms have the same skill level $\lambda^{*}$. It can be viewed as a Bayesian counterpart of the pooled OLS estimator. More rigorously, this prior is defined as $\lambda_{i}\sim\delta_{\lambda^{*}}$, where $\delta_{\lambda^{*}}$ is the Dirac delta function representing a degenerate distribution. The unknown $\lambda^{*}$ becomes another common parameter, similar to $\beta$, so I adopt a multivariate-normal-inverse-gamma prior on $\left(\left[\beta,\lambda^{*}\right]^{\prime},\sigma^{2}\right)$.
The flat prior (Flat) is specified as $f\left(\lambda_{i}\right)\propto1$, an uninformative prior with the posterior mode being the MLE estimate. Given the common parameters, there is no pooling from the cross-section, so we learn firm $i$'s skill $\lambda_{i}$ only from its own history.
The parametric prior (Param) combines cross-sectional information via a parametric distribution, such as a Gaussian distribution with unknown mean and variance, $\lambda_{i}\sim N\left(\mu,\omega^{2}\right)$. A normal-inverse-gamma hyperprior is further adopted for $\left(\mu,\omega^{2}\right)$. The parametric prior can be viewed as a limit case of the DPM prior when the scale parameter $\alpha\rightarrow0$, so there is only one component, and $\left(\mu,\omega^{2}\right)$ are directly drawn from the base distribution $G_{0}$. The choice of the hyperprior follows the suggestion by Basu2003 to match the parametric model with the DPM model such that “the predictive (or marginal) distribution of a single observation is identical under the two models.”
The nonparametric discrete prior (NP-disc) is modeled by a DP where $\lambda_{i}$ follows a flexible nonparametric distribution on a discrete support. This paper focuses on continuous $f$, which may be more sensible for the skills of young firms as well as other similar empirical studies. In this sense, comparing with NP-disc helps examine how much can be gained or lost from the continuity assumption and from the additional layer of mixture.
Finally, NP-R denotes the proposed nonparametric prior for random effects/coefficients models, and NP-C for correlated random effects/coefficients models. Both are flexible priors on continuous distributions, and NP-C allows $\lambda_{i}$ to depend on the initial condition of the firms.
The semiparametric predictors would reduce the estimation bias due to their flexibility while increasing the estimation variance due to their complexity. It is not transparent ex ante whether the parsimonious parametric predictors or the flexible semiparametric ones would perform better. Therefore, it is worthwhile to implement the Monte Carlo experiments and assess which predictor produces more accurate forecasts under which circumstances.
The specifications are summarized in Table (ref). $\beta_{0}$ is set to 0.8, as economic data usually exhibit some degree of persistence. The initial condition $y_{i0}$ is drawn from a standard normal distribution, which satisfies the moment condition in Assumption (ref)(3). Choices of $N=1000$ and $T=6$ are comparable with the young firm application. There are three experiments with different true distributions of $\lambda_{i}$. The first experiment features a degenerate $\lambda_{i}$ distribution, where all firms have the same skill level. Note that it does not satisfy Assumption (ref)(1-a) requiring the true $\lambda_{i}$ distribution to be continuous, and thus serves as a robustness check against the misspecification that the true $\lambda_{i}$ distribution is out of the prior support. The second experiment is based on a skewed distribution, a more realistic scenario in empirical studies. The third experiment incorporates a bimodal distribution with asymmetric weights on the two components. Various robustness checks are discussed in the Appendix.
I simulate 1,000 panel datasets in each setup. Forecasting performance, especially the relative rankings and magnitudes, is highly stable across repetitions. In each repetition, I generate 40,000 MCMC draws and discard the first 20,000 as burn-in. Based on graphical and statistical tests, the MCMC draws converge to a stationary distribution (see Appendix).
Table (ref) shows the forecasting comparison across predictors. When the $\lambda_{i}$ distribution is degenerate, Homog and NP-disc are the best, as expected. They are closely followed by NP-R and Param. Flat is considerably worse. When the $\lambda_{i}$ distribution is non-degenerate, there is a substantial gain from employing NP-R. In the bimodal case, NP-R far exceeds all alternatives. In the skewed case, Flat and Param are second best, yet still significantly inferior to NP-R. Homog and NP-disc yield the poorest forecasts, which suggests that their discrete supports may not be able to approximate the continuous $\lambda_{i}$ distribution in this case---even the nonparametric DP prior with countably infinite support may still be far from enough.
To investigate why we obtain better forecasts, Figure (ref) plots the posterior distribution of the $\lambda_{i}$ distribution for experiments Skewed and Bimodal. In the skewed case, NP-R better tracks the peak on the left and the tail on the right. In the bimodal case, NP-R nicely captures the M-shape. Therefore, the nonparametric prior flexibly approximates a vast set of distributions, which provides more precise estimates of the underlying $\lambda_{i}$ distributions and consequently more accurate density forecasts. This connection between distribution estimation and density forecasts reflects the theoretical results in Theorem (ref).
The general model accounts for three key features: multidimensional individual heterogeneity, cross-sectional heteroskedasticity, and correlated random coefficients. The exact specification is characterized and depicted in Table (ref).
In terms of multidimensional individual heterogeneity, $\lambda_{i}$ is now a 3-by-1 vector, and the corresponding covariates are composed of the intercept, time-specific $w_{t-1}^{(2)}$, and individual-time-specific $w_{i,t-1}^{(3)}$. In terms of correlated random coefficients, I adopt the conditional distribution following dunson2008kernel and ECT:9258097. They regard it as a challenging problem because this conditional distribution exhibits rapid changes in its shape, which considerably restricts the local sample size. Their original conditional distribution is one-dimensional, and I expand it to accommodate the three-dimensional $\lambda_{i}$ via a linear transformation. In terms of cross-sectional heteroskedasticity, I also let $\sigma_{i}^{2}$ interact with the initial conditions, and the functional form is modified from pelenis2014bayesian Case (ii). The modification guarantees that the $\sigma_{i}^{2}$ distribution is continuous with a large but bounded support above zero, and that the average signal-to-noise ratio is not far from 1. In addition, I consider the distribution of the innovations $v_{it}$ to be either normal or skewed. In the latter case, the normal likelihood function is misspecified. The $v_{it}$ distributions are standardized, i.e.\ $\mathbb{E}\left(v_{it}\right)=0$ and $\mathbb{V}\left(v_{it}\right)=1$, so we can identify $\sigma_{i}^{2}$.
The left two columns of Table (ref) describe the prior setups of $f_{\lambda}$ and $f_{\sigma^{2}}$. Due to cross-sectional heteroskedasticity and correlated random coefficients, the prior structures become more complicated. I further add Homosk-NP-C to examine whether it is practically relevant to model heteroskedasticity. The third column of Table (ref) assesses the forecasting performance under correct specification. Heterosk-NP-C is the most accurate density predictor. There are several messages if we compare density forecast performance across predictors. First, based on the comparison between Heterosk-NP-C and Homog/Homosk-NP-C, it is important to account for individual effects in both coefficients $\lambda_{i}$ and shock size $\sigma_{i}^{2}$. Second, comparing Heterosk-NP-C with Heterosk-Flat/Heterosk-Param, we see that the flexible nonparametric prior plays a significant role in enhancing density forecasts. Third, the difference between Heterosk-NP-C and Heterosk-NP-disc indicates that the discrete prior performs less satisfactorily when the underlying individual heterogeneity is continuous. Last, Heterosk-NP-R is less favorable than Heterosk-NP-C, which necessitates a careful modeling of the correlated random coefficient structure.
Under a misspecified $v_{it}$ distribution, the oracle knows the true distribution of $v_{it}$ and still serves as a legitimate benchmark for forecast evaluation. Although there is no theoretical guarantee, the proposed semiparametric method could still be helpful in density forecasts due to its flexibility---in the last column of Table (ref), the relative ranking is the same as the correctly specified case, and NP-C is still significantly better than the alternatives.
Studies have documented that young firm performance is affected by R&D and that different firms may react differently robb2014role,AkcigitKerr2010. In this empirical application, I examine this type of firm-specific latent heterogeneity from a density forecasting perspective. I use the confidential data from the Kauffman Firm Survey (KFS), which offers a large panel of startups (4,928 firms founded in 2004, nationally representative sample), a reasonable time span (2004-2011, one baseline survey and seven follow-up annual surveys), and detailed information on young firms. See robb2009overview for further description of the survey design.
I consider the general model with multidimensional individual heterogeneity in $\lambda_{i}$ and cross-sectional heteroskedasticity in $\sigma_{i}^{2}$. Following the firm dynamics literature, such as zarutskie2015did and AkcigitKerr2010, firm performance is measured by employment. From an economic point of view, young firms make a significant contribution to employment and job creation HaltiwangerJarminMiranda2012, and their struggle during the Great Recession may partly account for the jobless recovery afterward. Below, I focus on the following model specification, \[ \log\text{emp}_{it}=\beta\log\text{emp}_{i,t-1}+\lambda_{1i}+\lambda_{2i}\text{R\&D}_{i,t-1}+u_{it},\quad u_{it}\sim N\left(0,\sigma_{i}^{2}\right), \] where $\text{R\&D}_{it}$ is given by the ratio of a firm's R&D employment over its total employment. Other setups are discussed in the Appendix. An extension to a panel Tobit model as in tobit2018 could help accommodate firms' endogenous exit choice, which is left for future exploration.
The panel used for estimation spans from 2004 ($t=0$) to 2010 ($t=T$) with time dimension $T=6$. The data for 2011 ($t=T+1$) are reserved for pseudo-out-of-sample forecast evaluation. The sample is constructed as follows. First, for any $\left(i,t\right)$, if firm $i$'s R&D employment is greater than its total employment, there is an incompatibility issue, and the corresponding $\text{R\&D}_{it}$ is set to NA, which only affects 0.68% of the observations. Then, I only keep firms with long enough observations for identification in unbalanced panels. This results in a cross-sectional dimension $N=503$. The proportion of missing values is $\left(\#\mbox{missing obs}\right)/\left(NT\right)=9.32\%$. Here I consider unbalanced panels with randomly omitted observations (see Appendix), which helps incorporate more individuals into estimation and elicits more information for prediction. The descriptive statistics for $\log\text{emp}_{it}$ and $\text{R\&D}_{it}$ are summarized in Table (ref), and the corresponding densities are plotted in Figure (ref) in the Appendix. Both distributions are right skewed and may be multimodal, so we expect that the proposed predictors with nonparametric priors could perform well in this example.
The alternative priors are similar to those in the Monte Carlo simulation except for one additional prior, Heterosk-NP-C/R, where $\lambda_{i}$ can be correlated with $y_{i0}$ while $\sigma_{i}^{2}$ is independent with respect to $y_{i0}$. Then, I adopt an MGLR\textsubscript{x} prior on $f_{\lambda}$ and a DPM prior on $f_{l}$ for Heterosk-NP-C/R. The conditioning variable $y_{i0}$ is further standardized, which ensures numerical stability as the conditioning variables enter exponentially into the covariance function of the Gaussian process.
The first two columns in Table (ref) characterize the posterior estimates of the common parameter $\beta$. In most cases, the posterior means are mostly around $0.5\sim0.6$, which suggests that the young firm performance exhibits some degree of persistence, but the persistence is not strong. For Homog and NP-disc, their posterior means of $\beta$ are much larger. This may arise from the fact that homogeneous or discrete $\lambda_{i}$ structure may not be able to capture all individual effects, so these estimators may attribute the remaining individual effects to the persistence and thus overestimate $\beta$. NP-R also gives a large estimate of $\beta$. The reason is similar---if the true DGP features correlated random coefficients, the random coefficients model would miss the effect of the initial condition and misinterpret it as the persistence. In all scenarios, the posterior standard deviations are relatively small.
The last column in Table (ref) compares density forecasting performance. The overall best is Heterosk-NP-C/R. The main message is similar to the Monte Carlo of the general model---it is crucial to account for individual effects in both coefficients $\lambda_{i}$ and shock size $\sigma_{i}^{2}$ through a flexible nonparametric prior that acknowledges continuity and correlated random coefficients when the underlying individual heterogeneity has these features. Intuitively, the odds, given by the exponential of the difference in the LPS, indicate that Heterosk-NP-C/R produces density forecasts 32% (31%) more likely than Homog (Heterosk-Flat) does, on average.
Figures (ref) and (ref) (in the Appendix) provide the histograms of the probability integral transformation (PIT). While the LPS characterizes the relative ranks of predictors, the PIT complements the LPS and can be viewed as an absolute evaluation of how well the density forecasts coincide with the true (unobserved) conditional forecasting distributions given the current information set. Under the null hypothesis that the density forecasts coincide with the true DGP, the PITs are i.i.d.\ $U\left(0,1\right)$ and the histogram is close to a flat line diebold1998evaluating,amisano2013prediction. We can see that, in NP-C/R, NP-C, and Flat, the histogram bars are mostly within the confidence band, while other predictors yield apparent inverse-U shapes. The reason might be that the other predictors do not take correlated random coefficients into account but instead attribute their effects to the shock variance, which leads to more diffused predictive distributions.
Figure (ref) shows four types of firm-level predictive distributions: compared with Homog's Gaussian predictive distributions, NP-C/R is more concentrated in (a), more dispersed in (b), more skewed in (c), or exhibits extra kurtosis in (d). Figure (ref) in the Appendix regroups these predictive distributions by predictors. For Homog, all predictive distributions share the same Gaussian shape paralleling with each other. On the contrary, for NP-C/R, the predictive distributions exhibit fairly different shapes.
Figures (ref) and (ref) (in the Appendix) further aggregate the predictive distributions over sectors. It plots the predictive distributions of log average employment within each sector. Comparing Homog and NP-C/R across sectors, we can see several patterns. First, NP-C/R predictive distributions tend to be narrower. The reason is that NP-C/R tailors to each firm while Homog prescribes a general model to all the firms, so NP-C/R yields more precise predictive distributions. Second, NP-C/R predictive distributions have longer right tails, whereas Homog ones are in the standard bell shape. The long right tails in NP-C/R concur with the fact that good ideas are scarce. Finally, there is substantial heterogeneity in density forecasts across sectors. For sectors with relatively large average employment, e.g. construction, Homog pushes the forecasts down and hence systematically underpredicts their future employment, while NP-C/R respects this source of heterogeneity and significantly lessens the underprediction problem. On the other hand, for sectors with relatively small average employment, e.g. retail trade, Homog introduces an upward bias into the forecasts, while NP-C/R reduces this bias by flexibly estimating the underlying distribution of firm-specific heterogeneity.
The latent heterogeneity structure is presented in Figure (ref), which plots the joint distributions of the estimated individual effects and the conditional variable. For example, the pairwise relationship between $\lambda_{i1}$ and the standardized $y_{i0}$ is nonlinear and exhibits multiple components, which reassures our adoption of the nonparametric prior with correlated random coefficients. I also depict pairwise joint distributions involving $\hat{\sigma}_{i}^{2}$ in the Appendix. There does not seem to be much correlation between $\hat{\lambda}_{i}$ and $\hat{\sigma}_{i}^{2}$ and between $\hat{\sigma}_{i}^{2}$ and $y_{i0}$ (the latter is in line with the forecasting performance ranking where NP-C/R provides better density forecasts than NP-C does), which, together with sanity checks on (un)conditional correlation as well as a robustness check on density forecast performance (see Appendix), partially supports the assumption that conditioning on $y_{i0}$, $\lambda_{i}$ and $\sigma_{i}^{2}$ would be independent in this young firm sample.
This paper proposes a semiparametric Bayesian predictor, which performs well in density forecasts of individuals in a panel data setup. It considers the underlying distribution of individual effects and combines information from the whole panel in a flexible and efficient way. The full Bayesian procedure helps capture all sources of uncertainties and, together with the flexibility in the nonparametric Bayesian prior, cross-sectional heteroskedasticity, and correlated random coefficients, leads to more accurate density forecasts. The proposed method is theoretically appealing as the paper proves the posterior consistency of the estimates and the convergence of the density forecasts to the oracle in cross-sectional homoskedastic cases. The proposed method is also practically useful as demonstrated in the Monte Carlo simulations and an empirical application to young firm dynamics.