The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
97,240 characters
Density Forecasts in Panel Data Models: A Semiparametric Bayesian Perspective
\title{\textbf{Density Forecasts in Panel Data Models:}\\
\textbf{\Large{}A Semiparametric Bayesian Perspective}{\Large{}}\thanks{First version: November 15, 2016. Latest version: \protect\href{https://goo.gl/8zZZwn}{https://goo.gl/8zZZwn}.
I am indebted to my advisors, Francis X.\ Diebold and Frank Schorfheide,
as well as my committee members, Xu Cheng and Francis J.\ DiTraglia,
for much help and guidance throughout this project. I also thank Todd
Clark and Christian Hansen (Co-Editors), an anonymous Associate Editor,
and two anonymous referees for their constructive comments and suggestions.
I further benefited from many helpful discussions with St\'{e}phane
Bonhomme, Evan Chan, Benjamin Connault, Hyungsik R.\ Moon, Alexandre
Poirier, and seminar participants at University of Pennsylvania, FRB
Philadelphia, Federal Reserve Board, University of Virginia, Microsoft,
UC Berkeley, UC San Diego (Rady), Boston University, University of
Illinois at Urbana--Champaign, Princeton University, Libera Universit\`a
di Bolzano, University of Michigan, Universit\'{e} de Montr\'{e}al,
Emory University, Tilburg University, Erasmus University Rotterdam,
Tinbergen Institute, University of Toronto, as well as conference
participants at the Midwest Econometrics Group, NBER-NSF Seminar on
Bayesian Inference in Econometrics and Statistics, Bayesian Nonparametrics
Conference, Microeconometrics Class of 2017 Conference, Interactions
Workshop, First Italian Workshop of Econometrics and Empirical Economics,
International Association for Applied Econometrics Annual Conference,
and North American Winter Meeting of the Econometric Society. I would
also like to acknowledge the Kauffman Foundation and the NORC Data
Enclave for providing research support and access to the confidential
microdata. All remaining errors are my own.}}
\author{Laura Liu\thanks{Indiana University, \protect[email removed].}}
\date{September 26, 2021}
\maketitle
\vspace{-0.6cm}
\begin{abstract}
This paper constructs individual-specific density forecasts for a
panel of firms or households using a dynamic linear model with common
and heterogeneous coefficients as well as cross-sectional heteroskedasticity.
The panel considered in this paper features a large cross-sectional
dimension $N$ but short time series $T$. Due to the short $T$,
traditional methods have difficulty in disentangling the heterogeneous
parameters from the shocks, which contaminates the estimates of the
heterogeneous parameters. To tackle this problem, I assume that there
is an underlying distribution of heterogeneous parameters, model this
distribution nonparametrically allowing for correlation between heterogeneous
parameters and initial conditions as well as individual-specific regressors,
and then estimate this distribution by combining information from
the whole panel. Theoretically, I prove that in cross-sectional homoskedastic
cases, both the estimated common parameters and the estimated distribution
of the heterogeneous parameters achieve posterior consistency, and
that the density forecasts asymptotically converge to the oracle forecast.
Methodologically, I develop a simulation-based posterior sampling
algorithm specifically addressing the nonparametric density estimation
of unobserved heterogeneous parameters. Monte Carlo simulations and
an empirical application to young firm dynamics demonstrate improvements
in density forecasts relative to alternative approaches.
\end{abstract}
\begin{flushleft}
\textbf{JEL Codes: }C11, C14, C23, C53, L25
\par\end{flushleft}
\vspace{-0.5cm}
\begin{flushleft}
\textbf{Keywords: }Bayesian, Semiparametric Methods, Panel Data, Density
Forecasts, Posterior Consistency, Young Firm Dynamics
\par\end{flushleft}
\thispagestyle{empty}
\newpage{}
\setcounter{page}{1}
\section{Introduction\label{sec:Introduction}}
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 \citep{Hirano2002,gu2014unobserved} and bank
balance sheet data for bank stress tests \citep{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:
\begin{equation}
\underset{\text{performance}}{\underbrace{y_{it}}}=\beta y_{i,t-1}+\underset{\text{skill}}{\underbrace{\lambda_{i}}}+\underset{\text{shock}}{\underbrace{u_{it}}},\quad u_{it}\sim N\left(0,\sigma^{2}\right),\label{eq:motivation}
\end{equation}
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
\begin{align*}
p\left(\lambda_{i}\left|f,y_{i,0:T}\right.\right) & =\frac{p\left(\left.y_{i,1:T}\right|\lambda_{i}\right)f\left(\lambda_{i}\left|y_{i0}\right.\right)}{\int p\left(\left.y_{i,1:T}\right|\lambda_{i}\right)f\left(\lambda_{i}\left|y_{i0}\right.\right)d\lambda_{i}}.
\end{align*}
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{eq:motivation}), 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 \citet{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 \citet{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 \citet{LiuMoonSchorfheide2015}
and \citet{GuKoenker2016,gu2014unobserved}. \citet{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 \citet{LiuMoonSchorfheide2015} incorporate cross-sectional
and time-varying heteroskedasticity via a deterministic function of
observed conditioning variables.
\citet{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 \citet{gu2014unobserved} focus
on random effects models without this interaction.
In their recent paper, \citet{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 \citet{Lancaster01072002} (orthogonal
reparametrization and a flat prior); \citet{Chamberlain1999}, \citet{chib1999mcmc},
and \citet{sims2000using} (Gaussian prior); and \citet{chib2008panel}
(student-\emph{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,
\citet{Hirano2002} and \citet{jensen2015mutual} use linear panel
models with setups different from this paper. \citet{Hirano2002}
considers flexibility in the $u_{it}$ distribution instead of the
$\lambda_{i}$ distribution. \citet{jensen2015mutual} assume random
effects instead of correlated random effects. \citet{burda2013panel}
and \citet{10.2307/j.ctt5hhrfp} use a panel probit model and a panel
logit model, respectively.
In the frequentist literature, \citet{LI1998139}, \citet{delaigle2008deconvolution},
\citet{evdokimov2010identification}, and \citet{hu2017econometrics},
among others, have studied a similar deconvolution problem and estimated
the $\lambda_{i}$ distribution. Also see \citet{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 \citep{ghosh2003bayesian,Hjort2010,ghosal2017fundamentals}.
In particular, for unconditional density estimation, a recent paper
by \citet{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 \citep{norets2010approximation,norets2012bayesian}, a kernel
stick-breaking process \citep{ECT:9258097,pelenis2014bayesian,Norets2017},
or a probit stick-breaking process \citep{PatiDunsonTokdar2013}.
I adopt the \citet{PatiDunsonTokdar2013} approach and establish posterior
consistency for a multivariate conditional density estimator featuring
infinite location-scale mixtures with a probit stick-breaking process.
\label{paragraph: lit-inv-ineq}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, \citet{nguyen2013convergence} considers
the unobserved distribution on a discrete support, and \citet{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. \citet{AkcigitKerr2010} document
that R\&D intensive firms grow faster, especially for smaller firms.
\citet{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{sec:Model}
specifies the general panel data model, density forecasts, and nonparametric
Bayesian priors; Section \ref{sec:Theoretical-Properties} establishes
the posterior consistency of the estimates and the convergence of
the density forecasts to the oracle; Section \ref{sec:Simulation}
conducts Monte Carlo simulations; Section \ref{sec:Empirical-Application}
presents the empirical application to young firm dynamics; and Section
\ref{sec:Concluding-Remarks} concludes. Notations, proofs, algorithms,
and additional results are in the Appendix.
\section{Model\label{sec:Model}}
\subsection{General Panel Data Model\label{subsec:General-Panel-Data-1}}
The general panel data model with (correlated) random coefficients
and potential cross-sectional heteroskedasticity can be specified
as
\begin{gather}
y_{it}=\beta^{\prime}x_{i,t-1}+\lambda_{i}^{\prime}w{}_{i,t-1}+u_{it},\quad u_{it}\sim N\left(0,\sigma_{i}^{2}\right),\label{eq:general_panel}
\end{gather}
where $i=1,\cdots,N$, and $t=1,\cdots,T+h$. Similar to the baseline
setup in (\ref{eq:motivation}), $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 \citep{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{eq:motivation}).
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.
\subsection{Oracle and Feasible Predictors\label{subsec:oracle}}
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,
\begin{align}
f_{i,T+1}^{cond}\left(y\left|\vartheta,f,D_{i}\right.\right) & =\int\underset{\text{future shock}}{\underbrace{p\left(\left.y\right|h_{i},\vartheta,w_{iT},x_{iT}\right)}}\cdot\underset{\text{individual heterogeneity}}{\underbrace{p\left(h_{i}\left|\vartheta,f,D_{i}\right.\right)}}dh_{i},\label{eq:cond-pred}
\end{align}
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{eq:cond-pred}),
\begin{align*}
f_{i,T+1}^{oracle}\left(y\left|D_{i}\right.\right) & =f_{i,T+1}^{cond}\left(y\left|\vartheta_{0},f_{0},D_{i}\right.\right).
\end{align*}
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)$,
\begin{align*}
f_{i,T+1}^{sp}\left(y\left|D\right.\right) & =\int\underset{\text{shock \& heterogentity}}{\underbrace{f_{i,T+1}^{cond}\left(y\left|\vartheta,f,D_{i}\right.\right)}}\cdot\underset{\text{estimation uncertainty}}{\underbrace{d\Pi\left(\vartheta,f\left|D\right.\right)}}d\vartheta df.
\end{align*}
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.
\subsection{Nonparametric Bayesian Priors\label{subsec:Prior-Specification}}
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.
\subsubsection{Random Coefficients Model}
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,
\begin{align}
z_{i} & \sim\sum_{k=1}^{\infty}p_{k}N\left(\mu_{k},\Omega_{k}\right).\label{eq:dpm-def2}
\end{align}
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$.
\begin{equation}
\left(\mu_{k},\Omega_{k}\right)\sim G_{0},\text{ and }p_{k}\sim\zeta_{k}\prod_{j<k}\left(1-\zeta_{j}\right),\;\text{where }\zeta_{k}\sim\text{Beta}\left(1,\alpha\right).\label{eq:stick-breaking}
\end{equation}
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.
\subsubsection{Correlated Random Coefficients Model\label{subsec:mglr}}
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 \citet{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 \citet{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}$,
\begin{align*}
z_{i}|c_{i0} & \sim\sum_{k=1}^{\infty}p_{k}\left(c_{i0}\right)N\left(\mu_{k}\left[1,c_{i0}^{\prime}\right]^{\prime},\Omega_{k}\right).
\end{align*}
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
\begin{equation}
p_{k}\left(c_{i0}\right)=\Phi\left(\zeta_{k}\left(c_{i0}\right)\right)\prod_{j<k}\left(1-\Phi\left(\zeta_{j}\left(c_{i0}\right)\right)\right),\label{eq:cond-p}
\end{equation}
where stochastic function $\zeta_{k}$ is drawn from Gaussian process
$\zeta_{k}\sim GP\left(0,V_{k}\right)$ for $k=1,2,\cdots$. \citet{rodriguez2011nonparametric}
demonstrate the flexibility and computational simplicity of the probit
stick-breaking prior.
\label{paragraph: mglr-intuition}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).
\section{Theoretical Properties\label{sec:Theoretical-Properties}}
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 \citep{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 \citep{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.
\subsection{Identification\label{subsec:Identification}}
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{rem:id-homosk-re}.
\begin{assumption}
\label{assu:(model)}\emph{ (Identification: General Model)}
\end{assumption}
\begin{enumerate}
\item \emph{Model setup: Consider the panel data model in (\ref{eq:general_panel}),}
\begin{enumerate}
\item \emph{$\left(c_{i0},\lambda_{i},\sigma_{i}^{2}\right)$ are i.i.d.\ across
$i$. }
\item \emph{For all $t$, conditional on $\left(y_{it},c_{i,t-1}\right)$,
$x_{it}^{P*}$ is independent of $\left(\lambda_{i},\sigma_{i}^{2}\right)$.}
\item \emph{$\left(x_{i,0:T}^{O},w_{i,0:T}\right)$ are independent of $\left(\lambda_{i},\sigma_{i}^{2}\right)$. }
\item \emph{Conditioning on $c_{i0}$, $\lambda_{i}$ and $\sigma_{i}^{2}$
are independent of each other.}
\item \emph{Let $u_{it}=\sigma_{i}v_{it}$. $v_{it}\sim N\left(0,1\right)$
is i.i.d.\ across $i$ and $t$ and independent of $\left(c_{i,t-1},\lambda_{i},\sigma_{i}^{2}\right)$.}
\end{enumerate}
\item \emph{Identification:}
\begin{enumerate}
\item \emph{The characteristic functions of $\lambda_{i}|c_{i0}$ and $\sigma_{i}^{2}|c_{i0}$
are non-vanishing almost everywhere.}
\item \emph{For all i, $w_{i,0:T-1}$ has full rank $d_{w}$ almost everywhere. }
\item \emph{Let $\tilde{x}_{i,t-1}=x_{i,t-1}-\sum_{s=t+1}^{T}x_{i,s-1}w_{i,s-1}^{\prime}\left(\sum_{s=t+1}^{T}w_{i,s-1}w_{i,s-1}^{\prime}\right)^{-1}w_{i,t-1}$
given by orthogonal forward differencing. Then, the matrix $\mathbb{E}\big[\sum_{t=1}^{T-d_{w}}\tilde{x}_{i,t-1}\tilde{x}_{i,t-1}^{\prime}\big]$
has full rank $d_{x}$.}
\end{enumerate}
\end{enumerate}
\medskip{}
\label{paragraph: cond-corr}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).
\begin{thm}
\label{prop:(Identification)-1} \emph{(Identification: General Model)}
Under Assumption \ref{assu:(model)}, the common parameters $\beta$
and the conditional distribution of individual effects, $f_{\lambda}(\lambda_{i}|c_{i0})$
and $f_{\sigma^{2}}(\sigma_{i}^{2}|c_{i0})$, are all identified.
\end{thm}
\medskip{}
\noindent The argument is similar to \citet{ArellanoBover1995} and
\citet{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 \citep{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 \citep{GalambosSimonelli2004}
to obtain the identification of $f_{\sigma^{2}}$.
\begin{rem}
\label{rem:id-homosk-re}(1) For random coefficients models, Assumption
\ref{assu:(model)}(1-a) is replaced by ``$\left(\lambda_{i},\sigma_{i}^{2}\right)$
are independent of $c_{i0}$ and i.i.d.$\;$across $i$.''
\noindent (2) Under cross-sectional homoskedasticity, we can delete
Assumption \ref{assu:(model)}(1-d), get rid of $\sigma_{i}^{2}$
in conditions 1-a,b,c and $\sigma_{i}^{2}|c_{i0}$ in condition 2-a,
and replace condition 1-e by ``$u_{it}$ is i.i.d.\ across $i$
and $t$ and independent of $\left(c_{i,t-1},\lambda_{i}\right)$.''
\end{rem}
\subsection{Posterior Consistency\label{subsec:Posterior-Consistency}}
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.\label{par: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
\begin{equation}
g\left(\left.D_{i}\right|\vartheta,f\right)=\begin{cases}
\int p\left(\left.D_{i}\right|\vartheta,h_{i}\right)f\left(h_{i}\right)dh_{i}, & \text{if \emph{f} is an unconditional dist.,}\\
\int p\left(\left.\left.D_{i}\right\backslash c_{i0}\right|\vartheta,h_{i}\right)f\left(\left.h_{i}\right|c_{i0}\right)q_{0}\left(c_{i0}\right)dh_{i}, & \text{if \emph{f} is a conditional dist.,}
\end{cases}\label{eq:gen-g}
\end{equation}
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 \citep{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 \citep{PatiDunsonTokdar2013}.
Note that $q_{0}$ is only for theoretical derivation, and there is
no need to estimate it in practice.
\begin{thm}
\label{Thm: general} \emph{(Posterior Consistency: General Semiparametric
Model) }Suppose we have:
\end{thm}
\begin{enumerate}
\item \emph{Individual-specific likelihood:}
\begin{enumerate}
\item \emph{Kullback-Leibler (KL) property: For all $\epsilon>0$,
\[
\Pi\left(\left(\vartheta,f\right):\;D_{KL}\left(g\left(\left.D_{i}\right|\vartheta_{0},f_{0}\right)\parallel g\left(\left.D_{i}\right|\vartheta,f\right)\right)<\epsilon\right)>0.
\]
}
\item \emph{There exists $\delta_{\vartheta}>0$ such that for all $\left\Vert \vartheta_{1}-\vartheta_{2}\right\Vert _{2}<\delta_{\vartheta}$
and $f\in\mathcal{F}$, $\left\Vert g\left(\left.D_{i}\right|\vartheta_{1},f\right)-g\left(\left.D_{i}\right|\vartheta_{2},f\right)\right\Vert _{1}\le C_{g}\left\Vert \vartheta_{1}-\vartheta_{2}\right\Vert _{2},$
for some $C_{g}>0$ not depending on $f$.}
\item \emph{There exists an increasing function $\mathfrak{C}:\;\mathbb{R}_{\ge0}\mapsto\mathbb{R}_{\ge0}$
with $\lim_{x\rightarrow0}\mathfrak{C}\left(x\right)=0$ such that
for all $f\in\mathcal{F}$, $W_{2}\left(f,f_{0}\right)\le\mathfrak{C}\left(\left\Vert g\left(\left.D_{i}\right|\vartheta_{0},f\right)-g\left(\left.D_{i}\right|\vartheta_{0},f_{0}\right)\right\Vert _{1}\right)$.}
\end{enumerate}
\item \emph{Common parameters: There exists an exponentially consistent
sequence of tests $\varphi_{N}\left(D\right)$ for testing $H_{0}:\;\vartheta=\vartheta_{0},\ against\ H_{1}:\;\vartheta\in\Theta^{c},$
i.e.$\;$there exists a constant $C_{\varphi}>0$ such that
\[
(a)\;\mathbb{E}_{\vartheta_{0},f_{0}}\varphi_{N}\left(D\right)=O\left(e^{-C_{\varphi}N}\right),\;and\;(b)\;\sup_{\vartheta\in\Theta^{c},f\in\mathcal{F}}\mathbb{E}_{\vartheta,f}\left[1-\varphi_{N}\left(D\right)\right]=O\left(e^{-C_{\varphi}N}\right),
\]
where $\Theta^{c}\subset\Theta$ and $\vartheta\notin\Theta^{c}$.}
\item \emph{Distribution of individual heterogeneity: Its prior satisfies
a sieve property, i.e.$\;$there exists $\mathcal{F}_{N}\subset\mathcal{F}$
that can be partitioned as $\mathcal{F}_{N}=\cup_{j}\mathcal{F}_{N,j}$
such that, for all $\epsilon>0$,}
\begin{enumerate}
\item \emph{For some $\beta>0$, $\Pi_{f}\left(\mathcal{F}_{N}^{c}\right)=O\left(\exp\left(-\beta N\right)\right)$.}
\item \emph{For some $\gamma>0$, $\sum_{j}\sqrt{\mathcal{N}\left(\epsilon,\mathcal{F}_{N,j}\right)\Pi_{f}\left(\mathcal{F}_{N,j}\right)}=o\left(\exp\left(\left(1-\gamma\right)N\epsilon^{2}\right)\right)$,
where $\mathcal{N}\left(\epsilon,\mathcal{F}_{N,j}\right)$ is the
covering number of $\mathcal{F}_{N,j}$ by balls with radius $\epsilon$
in the $L_{1}$-norm.}
\end{enumerate}
\end{enumerate}
\emph{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.}
\medskip{}
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
\begin{align*}
& \Pi\left(\left.\vartheta\in\Theta_{\delta}^{c}\;\text{or}\;f\in\mathcal{F}_{\epsilon}^{c}\right|D\right)=\Pi_{\vartheta}\left(\left.\vartheta\in\Theta_{\delta}^{c}\right|D\right)+\Pi\left(\left.\vartheta\in\Theta_{\delta}\;\text{and}\;f\in\mathcal{F}_{\epsilon}^{c}\right|D\right)\\
= & \left.\left[\mathbb{P}\left(\vartheta\in\Theta_{\delta}^{c},D\right)+\mathbb{P}\left(\vartheta\in\Theta_{\delta},f\in\mathcal{F}_{\epsilon}^{c},D\right)\right]\right/\mathbb{P}\left(D\right)\\
= & \left.\left[\int_{\Theta_{\delta}^{c}\times\mathcal{F}}R_{N}\left(D,\vartheta,f\right)d\Pi\left(\vartheta,f\right)+\int_{\Theta_{\delta}\times\mathcal{F}_{\epsilon}^{c}}R_{N}\left(D,\vartheta,f\right)d\Pi\left(\vartheta,f\right)\right]\right/\int_{\Theta\times\mathcal{F}}R_{N}\left(D,\vartheta,f\right)d\Pi\left(\vartheta,f\right),
\end{align*}
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 \citep{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.
\label{paragraph: inv-ineq-id}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 \citep{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. \label{paragraph: ptwise-conv}Also,
the result only guarantees pointwise convergence in the space of the
distributions. For uniform forecasting performance in dynamic panel
data models, see \citet{LiuMoonSchorfheide2015}, which considers
an empirical Bayes setup with a nonparametric kernel estimate of the
marginal distribution of data.
\paragraph{Random Coefficients Model.\label{subsec:post-consist-re}}
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{Thm: general}, 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{assu:lag-y-re}(1-e)
below.
\begin{assumption}
\label{assu: (lag-y-y0)} \emph{(Covariates)}
\end{assumption}
\begin{enumerate}
\item \emph{$w_{i,0:T-1}$ is bounded.}
\item \emph{The eigenvalues of $\sum_{t}w_{i,t-1}w_{i,t-1}^{\prime}$ are
no less than some small $m_{w}>0$.}
\item \emph{$x_{i,0:T-1}^{O},$ $x_{i,0:T-1}^{P*}$, and $y_{i0}$ have
finite $4\left(1+\eta^{\prime}\right)$-th moments with $\eta^{\prime}>0$.}
\end{enumerate}
\medskip{}
\noindent 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{Thm: general}(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{Thm: general}(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{Thm: general}(1-a,b).
\begin{assumption}
\label{assu:lag-y-re} \emph{(Distribution of Individual Heterogeneity:
Random Coefficients)}
\end{assumption}
\begin{enumerate}
\item \emph{True distribution }$f_{0}$\emph{:}
\begin{enumerate}
\item \emph{$f_{0}\left(\lambda\right)$ is a continuous density.}
\item \emph{For some $M_{\lambda}>0$, $0<f_{0}\left(\lambda\right)\le M_{\lambda}$
for all $\lambda$.}
\item \emph{$\mathbb{E}_{f_{0}}\left[\log f_{0}\left(\lambda\right)\right]<\infty$.}
\item \emph{$\mathbb{E}_{f_{0}}\left[\log\frac{f_{0}\left(\lambda\right)}{\varphi_{\delta}\left(\lambda\right)}\right]<\infty$,
where $\varphi_{\delta}\left(\lambda\right)=\inf_{\left\Vert \lambda^{\prime}-\lambda\right\Vert _{2}<\delta}f_{0}\left(\lambda'\right)$,
for some $\delta>0$.}
\item \emph{For some $\eta>0$, $\mathbb{E}_{f_{0}}\left[\int\left\Vert \lambda\right\Vert _{2}^{2\left(1+\eta\right)}\right]<\infty$.}
\end{enumerate}
\item \emph{The base distribution of the DPM prior ($G_{0}$) follows a
multivariate-normal-inverse-Wishart distribution, where the degree
of freedom of the inverse Wishart component $\nu_{0}>\max\left(2d_{w},\left(2d_{w}+1\right)\left(d_{w}-1\right)\right)$.}
\end{enumerate}
\medskip{}
\noindent 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 \citet{Canale2017}, condition
2 further ensures the sieve property (Theorem \ref{Thm: general}(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$.
\begin{thm}
\label{prop:(lag-y-re)-1} \emph{(Posterior Consistency: Random Coefficients)}
Suppose we have:
\end{thm}
\begin{enumerate}
\item \emph{Model: Remark \ref{rem:id-homosk-re} for random coefficients
models with cross-sectional homoskedasticity.}
\item \emph{Covariates: $\left(x_{i,0:T},w_{i,0:T}\right)$ satisfies Assumption
\ref{assu: (lag-y-y0)}.}
\item \emph{Common parameters:}
\begin{enumerate}
\item \emph{$\vartheta_{0}$ is in the interior of $\text{supp}\left(\Pi_{\vartheta}\right)$.}
\item \emph{The domain of $\sigma^{2}$ is bounded by $\left[\underline{\sigma}^{2},\;\bar{\sigma}^{2}\right]$
for some $\underline{\sigma}^{2},\bar{\sigma}^{2}>0$.}
\end{enumerate}
\item \emph{Distributions of individual heterogeneity: $f_{0}$ and $\Pi_{f}$
satisfy Assumption \ref{assu:lag-y-re}.}
\end{enumerate}
\emph{Then, the posterior achieves consistency at $\left(\vartheta_{0},f_{0}\right)$.}
\paragraph{Correlated Random Coefficients Model.\label{subsec:post consist cre}}
$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{assu: (lag-y-cre)}(1-e) below.
\begin{assumption}
\label{assu: (lag-y-y0)-1} \emph{(Conditioning set) }$\mathcal{C}$
is compact, and $q_{0}\left(c_{0}\right)>0$ for all $c_{0}\in\mathcal{C}$.
\end{assumption}
\medskip{}
\noindent 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{assu: (lag-y-y0)}(1,3) for random coefficients
models.
\begin{assumption}
\label{assu: (lag-y-cre)} \emph{(Distribution of Individual Heterogeneity:
Correlated Random Coefficients) }
\end{assumption}
\begin{enumerate}
\item \emph{True distribution }$f_{0}$\emph{:}
\begin{enumerate}
\item \emph{$f_{0}\left(\cdot|\cdot\right)$ is jointly continuous in $\left(\lambda,c_{0}\right)$.}
\item \emph{For some $M_{\lambda}>0$, $0<f_{0}\left(\lambda|c_{0}\right)\le M_{\lambda}$
for all $\left(\lambda,c_{0}\right)$.}
\item \emph{$\mathbb{E}_{f_{0},q_{0}}\left[\log f_{0}\left(\lambda|c_{0}\right)\right]<\infty$.}
\item \emph{$\mathbb{E}_{f_{0},q_{0}}\left[\log\frac{f_{0}\left(\lambda|c_{0}\right)}{\varphi_{\delta}\left(\lambda|c_{0}\right)}\right]<\infty$,
where $\varphi_{\delta}\left(\lambda|c_{0}\right)=\inf_{\left\Vert \lambda^{\prime}-\lambda\right\Vert _{2}<\delta}f_{0}\left(\lambda'|c_{0}\right)$,
for some $\delta>0$.}
\item \emph{For some $\eta>0$, $\mathbb{E}_{f_{0},q_{0}}\left[\int\left\Vert \lambda\right\Vert _{2}^{2\left(1+\eta\right)}\right]<\infty$.}
\end{enumerate}
\item \emph{The base distribution of the MGLR}\textsubscript{\emph{x}}\emph{
prior ($G_{0}$) is characterized by }a\emph{ multivariate normal
distribution on $\text{vec}\left(\mu\right)$ and an inverse Wishart
distribution on $\Omega$, where the degree of freedom of the inverse
Wishart component $\nu_{0}>\max\left(2d_{w},\left(2d_{w}+1\right)\left(d_{w}-1\right)\right)$.}
\item \emph{Stick-breaking process: The covariance function for Gaussian
process can be specified as $V_{k}\left(c,\tilde{c}\right)=\tau\exp\left(-A_{k}\left\Vert c-\tilde{c}\right\Vert _{2}^{2}\right),$
where $\tau>0$ is a fixed number.}
\begin{enumerate}
\item \emph{The prior for $A_{k}$ has full support on $\mathbb{R}^{+}$.}
\item \emph{There exist $\beta$, $\gamma>0$ and a sequence $\delta_{N}=O\left(N^{-5/2}\left(\log N\right)^{2}\right)$
such that $\mathbb{P}\left(A_{k}>\delta_{N}\right)\le\exp\left(-N^{-\beta}k^{\left(\beta+2\right)/\gamma}\log k\right)$.}
\item \emph{For the same $\gamma$ as in condition 3-b, there exists an
increasing sequence $r_{N}\rightarrow\infty$ and $\left(r_{N}\right)^{d_{c0}}=o\left(N^{1-\gamma}\left(\log N\right)^{-\left(d_{c_{0}}+1\right)}\right)$
such that $\mathbb{P}\left(A_{k}>r_{N}\right)\le\exp\left(-N\right)$.}
\end{enumerate}
\end{enumerate}
\medskip{}
\noindent These conditions build on \citet{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{assu:lag-y-re}. 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.
\begin{thm}
\label{prop:(lag-y-cre)}\emph{ (Posterior Consistency: Correlated
Random Coefficients) }Suppose we have:
\end{thm}
\begin{enumerate}
\item \emph{Model: Remark \ref{rem:id-homosk-re}(2) for cross-sectional
homoskedastic models.}
\item \emph{Covariates: $\left(x_{i,0:T},w_{i,0:T}\right)$ satisfy Assumptions
\ref{assu: (lag-y-y0)}(2,3) and \ref{assu: (lag-y-y0)-1}.}
\item \emph{Common parameters: Theorem \ref{prop:(lag-y-re)-1}(3).}
\item \emph{Distributions of individual heterogeneity: $f_{0}$ and $\Pi_{f}$
satisfy Assumption \ref{assu: (lag-y-cre)}.}
\end{enumerate}
\emph{Then, the posterior achieves consistency at $\left(\vartheta_{0},f_{0}\right)$.}
\subsection{Density forecasts\label{subsec:theory-dfcst}}
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{prop:dfcst-general}
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.
\begin{thm}
\label{prop:dfcst} \emph{(Density Forecasts: (Correlated) Random
Coefficients with Cross-sectional Homoskedasticity)} Given conditions
in Theorem \ref{prop:(lag-y-re)-1} for random coefficients models
(or conditions in Theorems \ref{prop:(lag-y-cre)} and continuity
of $q_{0}\left(c_{0}\right)$ for correlated random coefficients models),
density forecasts converge to the oracle for all $i$ with $\mathbb{E}_{f_{0}}\left[\left.\left\Vert \lambda\right\Vert _{2}^{2}\right|c_{i0}\right]<\infty$,
i.e.$\;$given $i$, for all $\epsilon>0$, as $N\rightarrow\infty$,
\[
\mathbb{P}\left(\left.W_{2}\left(f_{i,T+1}^{cond},f_{i,T+1}^{oracle}\right)<\epsilon\right|D\right)\rightarrow1,
\]
in probability with respect to the true DGP.
\end{thm}
\medskip{}
\noindent The asymptotic convergence of aggregate-level density forecasts
can then be derived by summing individual-specific forecasts over
different subcategories.
\section{Monte Carlo Simulation\label{sec:Simulation}}
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.
\subsection{Forecast Evaluation and Alternative Predictors\label{subsec:Forecast-Evaluation-Methods}}
The accuracy of the density forecasts is measured by the log predictive
score (LPS) as suggested in \citet{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 \citet{AmisanoGiacomini2007}
(for the LPS) and \citet{timmermann2019comparing} (for panel data,
see their Section 2.6 on general loss functions) to examine the significance
in the LPS difference.
\begin{figure}[!t]
\begin{centering}
\caption{Alternative Predictors\label{fig:alt-predictor}}
\par\end{centering}
\medskip{}
\begin{centering}
\includegraphics[width=0.9\textwidth]{figures/prior_grayscale}
\par\end{centering}
\begin{singlespace}
\noindent \raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{}
For easier illustration, here I consider the baseline model with univariate
$\lambda_{i}$ and homoskedasticity. The black solid and teal dotted
lines represent two draws from each prior (except NP-disc, where the
teal one is also solid). Homog: Because $\lambda^{*}$ is unknown
}\emph{\footnotesize{}ex ante}{\footnotesize{}, the subgraph plots
two vertical lines representing two degenerate distributions with
different locations. Param: The subgraph contains two curves with
different means and variances. NP-disc: See Appendix for a formal
definition of the DP and how it relates to the DPM.}{\footnotesize\par}
\end{singlespace}
\end{figure}
Different predictors can be interpreted as different priors on the
distribution of $\lambda_{i}$. As these priors are distributions
over distributions, Figure \ref{fig:alt-predictor} 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 \citet{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 \emph{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.
\subsection{Baseline Model with Random Effects\label{subsec:baseline-model}}
The specifications are summarized in Table \ref{tab:Simulation-Setup:-Baseline}.
$\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{assu: (lag-y-y0)}(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{assu:lag-y-re}(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.
\begin{table}[!t]
\caption{Simulation Setup: Baseline Model with Random Effects\label{tab:Simulation-Setup:-Baseline}}
\medskip{}
\centering{}
\begin{tabular}{ll}
\hline
\hline Law of motion & $y_{it}=\beta y_{i,t-1}+\lambda_{i}+u_{it},\;u_{it}\sim N\left(0,\sigma^{2}\right)$\tabularnewline
Common parameters & $\beta_{0}=0.8,\;\sigma_{0}^{2}=\frac{1}{4}$\tabularnewline
Initial conditions & $y_{i0}\sim N\left(0,1\right)$\tabularnewline
Sample size & $N=1000,\;T=6$\tabularnewline
\hline
Random Effects: & \tabularnewline
Degenerate & $\lambda_{i}=0$\tabularnewline
Skewed & $\lambda_{i}\sim\frac{1}{9}N\left(2,\frac{1}{2}\right)+\frac{8}{9}N\left(-\frac{1}{4},\frac{1}{2}\right)$,
so $\mathbb{V}\left(\lambda_{i}\right)=1$\tabularnewline
Bimodal & $\lambda_{i}\sim\left(0.35N\left(0,1\right)+0.65N\left(10,1\right)\right)/\sqrt{1+10^{2}\cdot0.35\cdot0.65}$,
so $\mathbb{V}\left(\lambda_{i}\right)=1$\tabularnewline
\hline
\end{tabular}
\end{table}
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{tab:Forecast-Evaluation:-Benchmark} 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.
\begin{table}[!t]
\caption{Density Forecast Evaluation: Baseline Model with Random Effects\label{tab:Forecast-Evaluation:-Benchmark}}
\medskip{}
\begin{centering}
\begin{tabular}{>{\raggedright}p{0.75in}>{\raggedleft}p{0.75in}>{\raggedleft}p{0.75in}>{\raggedleft}p{0.75in}}
\hline
\hline & \multicolumn{1}{r}{Degenerate} & \multicolumn{1}{r}{Skewed} & \multicolumn{1}{r}{Bimodal}\tabularnewline
\hline
\emph{Oracle} & \emph{-725}\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}} & \emph{-798}\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}} & \emph{-766}\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}}\tabularnewline
\hline
Homog & \textbf{-0.2}{\footnotesize{}{*}{*}{*}} & -193{\footnotesize{}{*}{*}{*}} & -424{\footnotesize{}{*}{*}{*}}\tabularnewline
Flat & -102{\footnotesize{}{*}{*}{*}} & -7{\footnotesize{}{*}{*}{*}} & -38{\footnotesize{}{*}{*}{*}}\tabularnewline
Param & -4\textcolor{white}{\footnotesize{}{*}{*}{*}} & -1{\footnotesize{}{*}{*}{*}} & -34{\footnotesize{}{*}{*}{*}}\tabularnewline
NP-disc & \textbf{-0.2}{\footnotesize{}{*}{*}{*}} & -206{\footnotesize{}{*}{*}{*}} & -40{\footnotesize{}{*}{*}{*}}\tabularnewline
NP-R & -4\textcolor{white}{\footnotesize{}{*}{*}{*}} & \textbf{-0.3}\textcolor{white}{\footnotesize{}{*}{*}{*}} & \textbf{-6}\textcolor{white}{\footnotesize{}{*}{*}{*}}\tabularnewline
\hline
\end{tabular}
\par\end{centering}
\medskip{}
\emph{\footnotesize{}Notes:}{\footnotesize{} The density forecasts
are assessed by the LPS and a test combining \citet{AmisanoGiacomini2007}
and \citet{timmermann2019comparing}. For the oracle predictor, the
table reports the exact values of $\text{LPS}\cdot N$ (averaged over
1,000 Monte Carlo samples). For other predictors, the table reports
their differences from the oracle. The tests compare other feasible
predictors with NP-R, with significance levels indicated by {*}: 10\%,
{*}{*}: 5\%, and {*}{*}{*}: 1\%. The entries in bold indicate the
best feasible predictor in each column. }{\footnotesize\par}
\end{table}
To investigate why we obtain better forecasts, Figure \ref{fig:Estimated-:-Benchmark}
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{prop:dfcst}.
\begin{figure}[!t]
\caption{$f_{0}$ vs $\Pi_{f}\left(f\left|y_{1:N,0:T}\right.\right):$ Baseline
Model with Random Effects\label{fig:Estimated-:-Benchmark}}
\medskip{}
\begin{centering}
\begin{tabular}{cccc}
\multicolumn{2}{c}{(a) Skewed} & \multicolumn{2}{c}{(b) Bimodal}\tabularnewline
Param & NP-R & Param & NP-R\tabularnewline
[-0.75ex]\includegraphics[width=0.23\textwidth]{figures/modrun_103_pi_shade_lambda2est4} & \includegraphics[width=0.23\textwidth]{figures/modrun_103_pi_shade_lambda2est6} & \includegraphics[width=0.23\textwidth]{figures/modrun_103_pi_shade_lambda4est4} & \includegraphics[width=0.23\textwidth]{figures/modrun_103_pi_shade_lambda4est6}\tabularnewline
\end{tabular}
\par\end{centering}
\emph{\footnotesize{}Notes:}{\footnotesize{} The subgraphs are constructed
from the estimation results of one of the 1,000 repetitions. The black
solid lines represent the true $\lambda_{i}$ distributions, $f_{0}$.
The teal bands show the posterior distributions of $f$, $\Pi_{f}\left(f\left|y_{1:N,0:T}\right.\right)$.}{\footnotesize\par}
\end{figure}
\subsection{General Model\label{subsec:General-model}}
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{tab:Simulation-Setup:-General}.
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 \citet{dunson2008kernel}
and \citet{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 \citet{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}$.
\begin{table}[!t]
\caption{Simulation Setup: General Model \label{tab:Simulation-Setup:-General}}
\medskip{}
\begin{centering}
\begin{tabular}{l>{\raggedright}p{0.65\textwidth}}
\hline
\hline Law of motion & $y_{it}=\beta y_{i,t-1}+\lambda_{i}^{\prime}w_{i,t-1}+u_{it},\;u_{it}=\sigma_{i}v_{it}$\tabularnewline
Covariates & $w_{i,t-1}=[1,w_{t-1}^{(2)},w_{i,t-1}^{(3)}]^{\prime}$, $w_{t-1}^{(2)}\sim N\left(0,1\right)\mathbf{1}\left(\left|w_{t-1}^{(2)}\right|\le10\right)$,
$w_{i,t-1}^{(3)}\sim\text{Ga}\left(1,1\right)\mathbf{1}\left(w_{i,t-1}^{(3)}\le10\right)$\tabularnewline
Common parameters & $\beta_{0}=0.8$\tabularnewline
Initial conditions & $y_{i0}\sim U\left(0,1\right)$\tabularnewline
Corr.\ random coef.\ & $\lambda_{i}|y_{i0}\sim e^{-2y_{i0}}N\left(y_{i0}v,0.1^{2}vv'\right)+\left(1-e^{-2y_{i0}}\right)N\left(y_{i0}^{4}v,0.2^{2}vv'\right)$,
$v=\left[1,2,-1\right]^{\prime}$\tabularnewline
Cross-sec.\ heterosk.\ & $\sigma_{i}^{2}|y_{i0}\sim\left[0.454\left(y_{i0}+0.5\right)^{2}\cdot\text{IG}\left(51,40\right)+10^{-6}\right]\cdot\mathbf{1}\left(\sigma_{i}^{2}\le10^{6}\right)$\tabularnewline
Sample size & $N=1000,\;T=6$\tabularnewline
\hline
\multicolumn{2}{l}{Innovation distributions:}\tabularnewline
Normal & $v_{it}\sim N\left(0,1\right)$\tabularnewline
Skewed & $v_{it}\sim\frac{1}{9}N\left(2,\frac{1}{2}\right)+\frac{8}{9}N\left(-\frac{1}{4},\frac{1}{2}\right)$\tabularnewline
\hline
\end{tabular}
\par\end{centering}
\medskip{}
\medskip{}
\medskip{}
\begin{centering}
\begin{tabular}{cccc}
\multicolumn{2}{c}{(a) $\left.\lambda_{i1}\right|y_{i0}$} & \multicolumn{2}{c}{(b) $\left.\sigma_{i}^{2}\right|y_{i0}$}\tabularnewline
[-0.25ex]\includegraphics[width=0.23\textwidth]{figures/dgp2_lamb1_joint} & \includegraphics[width=0.23\textwidth]{figures/dgp2_lamb1_cond} & \includegraphics[width=0.23\textwidth]{figures/dgp2_sigma2_joint} & \includegraphics[width=0.23\textwidth]{figures/dgp2_sigma2_cond}\tabularnewline
\end{tabular}
\par\end{centering}
\emph{\footnotesize{}Notes:}{\footnotesize{} In the left two panels,
$\lambda_{i1}$ is the coefficient on $w_{i,t-1}^{\left(1\right)}=1$
and can be interpreted as the heterogeneous intercept. In the second
and fourth panels, the black solid / teal dashed / orange dotted lines
are conditional on $y_{i0}=0.25,\;0.5,\;\text{and }0.75$, respectively.
As $y_{i0}\sim U\left(0,1\right)$, the conditional distribution equals
the joint distribution for all $y_{i0}\in\left[0,1\right]$, i.e.$\;f\left(\left.\lambda_{i1}\right|y_{i0}\right)=f\left(\left.\lambda_{i1}\right|y_{i0}\right)q_{0}\left(y_{i0}\right)=f\left(\lambda_{i1},y_{i0}\right)$.}{\footnotesize\par}
\end{table}
The left two columns of Table \ref{tab:Forecast-Evaluation:-General}
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{tab:Forecast-Evaluation:-General}
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.
\label{paragraph: misspec}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{tab:Forecast-Evaluation:-General},
the relative ranking is the same as the correctly specified case,
and NP-C is still significantly better than the alternatives.
\begin{table}[!t]
\caption{Prior Structures and Density Forecast Evaluation: General Model\label{tab:Forecast-Evaluation:-General}}
\medskip{}
\begin{centering}
\begin{tabular}{ll|ll|rr}
\hline
\hline & & $f_{\lambda}$ & $f_{\sigma^{2}}$ (or $f_{l}$) & \multicolumn{1}{r}{Normal $v_{it}$\textcolor{white}{\footnotesize{}{*}}} & \multicolumn{1}{r}{Skewed $v_{it}$}\tabularnewline
\hline
\emph{Oracle} & & \emph{Known} & \emph{Known} & \emph{-974}\textbf{\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}}} & \emph{-965}\textbf{\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}}}\tabularnewline
\hline
Homog & & $=\delta_{\lambda^{*}}$ & $f_{\sigma^{2}}=\delta_{\sigma^{2*}}$ & -407{\footnotesize{}{*}{*}{*}} & -417{\footnotesize{}{*}{*}{*}}\tabularnewline
Homosk & NP-C & $\sim$ MGLR\textsubscript{x} & $f_{\sigma^{2}}=\delta_{\sigma^{2*}}$ & -134{\footnotesize{}{*}{*}{*}} & -146{\footnotesize{}{*}{*}{*}}\tabularnewline
\hline
Heterosk & Flat & $\propto1$ & $f_{\sigma^{2}}\propto1$ & -384{\footnotesize{}{*}{*}{*}} & -366{\footnotesize{}{*}{*}{*}}\tabularnewline
& Param & $=$ Normal & $f_{\sigma^{2}}=$ IG & -79{\footnotesize{}{*}{*}{*}} & -78{\footnotesize{}{*}{*}{*}}\tabularnewline
& NP-disc & $\sim$ DP & $f_{l}\sim$ DP & -79{\footnotesize{}{*}{*}{*}} & -78{\footnotesize{}{*}{*}{*}}\tabularnewline
& NP-R & $\sim$ DPM & $f_{l}\sim$ DPM & -229{\footnotesize{}{*}{*}{*}} & -224{\footnotesize{}{*}{*}{*}}\tabularnewline
& NP-C & $\sim$ MGLR\textsubscript{x} & $f_{l}\sim$ MGLR\textsubscript{x} & \textbf{-70}\textcolor{white}{\footnotesize{}{*}{*}{*}} & \textbf{-71}\textcolor{white}{\footnotesize{}{*}{*}{*}}\tabularnewline
\hline
\end{tabular}
\par\end{centering}
\medskip{}
\raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{} The prior
structure of Heterosk-Param is detailed in the Appendix. For density
forecast evaluation, see the description in Table \ref{tab:Forecast-Evaluation:-Benchmark}.
Here the tests are conducted with respect to Heterosk-NP-C.}{\footnotesize\par}
\end{table}
\section{Empirical Application: Young Firm Dynamics\label{sec:Empirical-Application}}
Studies have documented that young firm performance is affected by
R\&D and that different firms may react differently \citep{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 \citet{robb2009overview} for further description
of the survey design.
\subsection{Model Specification}
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 \citet{zarutskie2015did}
and \citet{AkcigitKerr2010}, firm performance is measured by employment.
From an economic point of view, young firms make a significant contribution
to employment and job creation \citep{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 \citet{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{tab:descrip}, and the corresponding densities are plotted
in Figure \ref{fig:descrip} 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.
\begin{table}[!t]
\caption{Descriptive Statistics of Observables\label{tab:descrip}}
\medskip{}
\centering{}
\begin{tabular}{lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}lr@{\extracolsep{0pt}.}l}
\hline
\hline & \multicolumn{2}{c}{10\%} & \multicolumn{2}{c}{Mean} & \multicolumn{2}{c}{Med$.$} & \multicolumn{2}{c}{90\%} & \multicolumn{2}{c}{SD} & \multicolumn{2}{c}{Skew$.$} & \multicolumn{2}{c}{Kurt$.$}\tabularnewline
\hline
log emp & 0&69 & 1&59 & 1&39 & 2&20 & 1&02 & 0&59 & 3&42\tabularnewline
R\&D & 0&00 & 0&27 & 0&14 & 0&50 & 0&32 & 1&18 & 3&25\tabularnewline
\hline
\end{tabular}
\end{table}
\subsection{Results}
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{tab:Forecast-Evaluation:app}
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{tab:Forecast-Evaluation:app} 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.
\begin{table}[!t]
\caption{Parameter Estimation and Density Forecast Evaluation: Young Firm Dynamics\label{tab:Forecast-Evaluation:app}}
\medskip{}
\begin{centering}
\begin{tabular}{ll|rc|r}
\hline
\hline & & \multicolumn{2}{c|}{$\beta$} & \multicolumn{1}{c}{LPS{*}N}\tabularnewline
\cline{3-4} \cline{4-4}
& & Mean & SD & \tabularnewline
\hline
\emph{Heterosk} & \emph{NP-C/R} & \emph{0.50} & \emph{0.02} & \textbf{\emph{-195}}\textcolor{white}{\emph{\footnotesize{}{*}{*}{*}}}\tabularnewline
\hline
Homog & & 0.88 & 0.02 & -139{\footnotesize{}{*}{*}{*}}\tabularnewline
Homosk & NP-C & 0.48 & 0.02 & -113{\footnotesize{}{*}{*}{*}}\tabularnewline
\hline
Heterosk & Flat & 0.19 & 0.07 & -134{\footnotesize{}{*}{*}{*}}\tabularnewline
& Param & 0.62 & 0.07 & -63{\footnotesize{}{*}{*}{*}}\tabularnewline
& NP-disc & 0.92 & 0.01 & -88{\footnotesize{}{*}{*}{*}}\tabularnewline
& NP-R & 0.74 & 0.04 & -20{\footnotesize{}{*}{*}}\textcolor{white}{\footnotesize{}{*}}\tabularnewline
& NP-C & 0.53 & 0.03 & -6{\footnotesize{}{*}}\textcolor{white}{\footnotesize{}{*}{*}}\tabularnewline
\hline
\end{tabular}
\par\end{centering}
\medskip{}
\raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{} See the
description of Table \ref{tab:Forecast-Evaluation:-Benchmark} for
density forecast evaluation. Here Heterosk-NP-C/R is the benchmark
for both normalization and significance tests. For Heterosk-NP-C/R,
the table reports the exact values of $\text{LPS}\cdot N$. For other
predictors, the table reports their differences from Heterosk-NP-C/R.}{\footnotesize\par}
\end{table}
Figures \ref{fig:PIT} and \ref{fig:PIT-1} (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 \citep{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.
\begin{figure}[!t]
\begin{centering}
\caption{PIT\label{fig:PIT}}
\par\end{centering}
\medskip{}
\begin{centering}
\begin{tabular}{cc}
Homog & NP-C/R\tabularnewline
[-0.25ex]\includegraphics[width=0.35\textwidth]{figures/modrun4_pitmain_Homog} & \includegraphics[width=0.35\textwidth]{figures/modrun4_pitmain_NPCR}\tabularnewline
\end{tabular}
\par\end{centering}
\raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{} Teal lines
indicate the confidence interval. See Appendix for PITs of all predictors.}{\footnotesize\par}
\end{figure}
Figure \ref{fig:pred-dist-1} 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{fig:pred-dist-1-1}
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{fig:pred-dist} and \ref{fig:pred-dist-2} (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.
\begin{figure}[!t]
\begin{centering}
\caption{Predictive Distributions: Firm-level, 4 Types\label{fig:pred-dist-1}}
\par\end{centering}
\medskip{}
\begin{centering}
\begin{tabular}{cccc}
(a) & (b) & (c) & (d)\tabularnewline
[-0.25ex]\includegraphics[width=0.23\textwidth]{figures/modrun4_pred_dist_firm39} & \includegraphics[width=0.23\textwidth]{figures/modrun4_pred_dist_firm385} & \includegraphics[width=0.23\textwidth]{figures/modrun4_pred_dist_firm371} & \includegraphics[width=0.23\textwidth]{figures/modrun4_pred_dist_firm103}\tabularnewline
\end{tabular}
\par\end{centering}
\raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{} The black
solid (teal dotted) lines are the predictive distributions via the
NP-C/R (Homog).}{\footnotesize\par}
\end{figure}
\begin{figure}[!t]
\begin{centering}
\caption{Predictive Distributions: Aggregated by Sectors\label{fig:pred-dist}}
\par\end{centering}
\medskip{}
\begin{centering}
\begin{tabular}{cc}
Construction & Retail Trade\tabularnewline
[-0.75ex]\includegraphics[width=0.35\textwidth]{figures/modrun4pred_dist_ind_3_main} & \includegraphics[width=0.35\textwidth]{figures/modrun4pred_dist_ind_8_main}\tabularnewline
\end{tabular}
\par\end{centering}
\raggedright{}\emph{\footnotesize{}Notes:}{\footnotesize{} The black
solid (teal dotted) lines are the predictive distributions via the
NP-C/R (Homog). See Appendix for predictive distributions of all sectors.}{\footnotesize\par}
\end{figure}
The latent heterogeneity structure is presented in Figure \ref{fig:joint-dist},
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. \label{paragraph: condcorr_lambdasigma2}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 \emph{$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.
\begin{figure}[!t]
\begin{centering}
\caption{Joint Distributions: $\hat{\lambda}_{i}$ and $y_{i0}$\label{fig:joint-dist}}
\par\end{centering}
\medskip{}
\begin{tabular}{cc}
\rotatebox{90}{\hspace*{2.2cm} \footnotesize{$\hat\lambda_{i1}$}}\hspace{-.5cm} & \includegraphics[bb=0bp 0bp 432bp 432bp,width=0.3\textwidth]{figures/corr41}\tabularnewline
[-3ex] & \footnotesize{Standardized $y_{i0}$}\tabularnewline
\end{tabular}\hspace{-.3cm}
\begin{tabular}{cc}
\rotatebox{90}{\hspace*{2.2cm} \footnotesize{$\hat\lambda_{i2}$}}\hspace{-.5cm} & \includegraphics[bb=0bp 0bp 432bp 432bp,width=0.3\textwidth]{figures/corr42}\tabularnewline
[-3ex] & \footnotesize{Standardized $y_{i0}$}\tabularnewline
\end{tabular}\hspace{-.3cm}
\begin{tabular}{cc}
\rotatebox{90}{\hspace*{2.2cm} \footnotesize{$\hat\lambda_{i2}$}}\hspace{-.5cm} & \includegraphics[bb=0bp 0bp 432bp 432bp,width=0.3\textwidth]{figures/corr12}\tabularnewline
[-3ex] & \footnotesize{$\hat\lambda_{i1}$}\tabularnewline
\end{tabular}\emph{\footnotesize{}}\\
\emph{\footnotesize{}\vspace{0.15cm}
}{\footnotesize\par}
\emph{\footnotesize{}Notes:}{\footnotesize{} $\lambda_{i1}$ is the
heterogeneous intercept, and $\lambda_{i2}$ is the heterogeneous
coefficient on R\&D.}{\footnotesize\par}
\end{figure}
\section{Conclusion\label{sec:Concluding-Remarks}}
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.
\bibliographystyle{ecca}
\bibliography{dp}