Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,873 characters · 30 sections · 54 citation commands
Forecasting with Bayesian Grouped Random Effects in Panel Data
{5pt} {7pt}
JEL CLASSIFICATION: C11, C14, C23, C53, G31
KEY WORDS: Panel Data; Grouped Heterogeneity; Random Effects; Dirichlet Process; Set Forecast; Density Forecast; Investment
\thispagestyle{empty} \setcounter{page}{0}
With the increasing availability of panel data, many works have examined and demonstrated its central role in the empirical research throughout the social and business sciences. Analysis of panel data has various edges over that of pure cross-sectional or time-series data. The most important one is that the panel data provide researchers with a flexible way to model both heterogeneity among individuals, firms, regions, and countries and possible structural changes over time. Apart from the principal role in the model estimation, it is interesting and essential to study their relevance for forecasting. Among novel methods emerged recently, the latent group structure in the heterogeneity attracts wide attention. In this paper, we allow for grouped patterns of unobserved heterogeneity in the dynamics panel data models and evaluate whether this latent structure improves the predictive performance in an extensive collection of short time series.
In the dynamics panel data model, it is common to assume that each cross-sectional unit has unique intercept. This assumption introduces a large number of parameters that become a burden in estimation. In models that have as many parameters as individual units, fixed effects estimators are known to suffer from the “incidental parameters" problem neyman1948, which can bring about significant biases in estimates of common parameters. This problem becomes severe in short panels even if the number of units goes to infinity chamberlain1980,nickell1981, and the fixed-effects themselves are often poorly estimated. An unreliable estimate leads to concerns about the predictive power of panel data models as inaccurate estimates affect forecasts in all aspects.
To address this issue\footnote{Another important strand of literature implements generalized method of moments (GMM) methods to eliminate bias, see arellano1991, arellano1995 and blundell1998. Though successfully solved the “incidental parameters" problem, this set of methods doesn't allow for any latent group structure.}, econometricians attempt to reduce the number of unknown parameters by dividing units into a finite number of groups. The premise of this idea is that units in the same group share the unit-specific parameters. Previous works include bonhomme2015, ando2016, su2016, bester2016, su2019, bonhomme2019, and cheng2019. Moreover, finite mixture model provide a well-known probabilistic approach to model-based clustering mcnicholas2010,fruhwirth2011a. With a finite number of groups, econometricians could avoid “incidental parameter” problem under several particular assumptions and derive consistent estimators for the common parameters.
However, the convenience of the group structure does not come without any cost. The number of groups is an unknown but fixed quantity, and the need to specify the number in advance is deemed one of the significant drawbacks of applying these methods in a clustering context. Many methods have been suggested to estimate the optimal number a posteriori from the data such as BIC keribin2000,bonhomme2015, marginal likelihoods fruhwirth2004, or the integrated classification likelihood biernacki2000. Bayesian approaches sometimes pursue a similar strategy, often adding the DIC to the list of model choice criteria, e.g., celeux2006 and kim2019. If both $N$ and $T$ are large enough, the information criterion could select the true group structure. However, with a short time span, these criteria might fail to achieve their goal. As noted in bonhomme2015, the choice of the number of groups is crucial to estimation and inference for model parameters. Misspecification in group number forces the algorithm to consider incorrect group membership. We will later show that it is the information criteria that substantially affects the performance of Grouped Fixed Effects (GFE) estimator proposed by BM.
The contributions of this paper are fourfold. First, closely following kim2019 and liu2020, we develop a posterior sampling algorithm that addresses the nonparametric estimation of latent grouped effects and proposes Bayesian Grouped Random Effects (BGRE) estimator. The number of groups is treated as an unknown parameter that is estimated jointly with the component-specific parameters under the assumption that group membership remains constant over time. Instead of using the Finite Mixture model, which needs to preset the number of groups, we use Dirichlet Process (DP) prior, in particular the stick-breaking prior, that allows for infinite potential groups. The entire posterior sampler builds upon the blocked Gibbs sampling\footnote{Unlike the Pólya urn Gibbs sampler escobar1995, blocked Gibbs sampler approach avoids marginalizing over the prior and thus allows for direct sampling of the nonparametric posterior, leading to computational and inferential advantages.} proposed by ishwaran2001.
Second, we leverage the researcher's prior knowledge of the latent group structure to improve the estimation and forecasting. In particular, we summarize and incorporate the information of subjective group structure in the prior distribution of the membership probabilities. If the subjective prior on the group structure is more precise than the random guess, even with incorrect presumed number of groups, we show that including it in the prior improves the performance of the BGRE estimators as it guides the group membership estimates.
Third, we explore the potential link between the proposed BGRE estimators and unsupervised machine learning method. Theoretically, we show that our block Gibbs sampler for the BGRE estimator is closely related to the Kmeans algorithm macqueen1967 under certain assumptions. In particular, both algorithms assign units to the closest centroid when forming the clusters and recalculate the means of the new cluster afterward. To compare the performance of clustering, we modify our algorithm to incorporate Kmeans and construct a two-step BGRE estimator where individuals are clustered in the first step using Kmeans, and the group-specific heterogeneity is estimated in the second step. In the simulation section, we document that our BGRE estimators dominate the two-step GRE estimator in terms of the performance of both clustering and forecasting. We also find that the two-step GRE estimators with Kmeans algorithm severely underestimate the number of groups under all data generating processes, whereas BGRE estimators deliver accurate estimates.
Last but not least, we examine the performance of BGRE estimators using various sets of simulated data and real data. The Monte Carlo study presents that grouped heterogeneity brings gains in estimating group structure and one-step ahead point, set, and density forecasting relative to commonly used predictors with different parametric priors on individual effects. In particular, our estimators outperform BM's Grouped Fixed Effects (GFE) estimator in various settings of the data generating process. The better performance is primarily due to the accurate estimate of the group structure. Regarding other predictors, we show that failing to model group structure and to pool information across units severely deteriorates the results for both estimation and forecasting. Finally, we use our method to forecast the investment rate across a broad range of firms. The BGRE estimators offer better performance than the standard panel data models in forecasting. This reveals that incorporating the latent group structure provides a great amount of flexibility and improves the predictive power of the underlying panel data model.
Our paper relates to several branches of the literature. Our work is closely related to BM, KW, and bonhomme2019. All of these three papers aim to estimate the unobserved heterogeneity in a linear dynamic panel data model and develop statistical inference methods. BM estimates the parameters of the model using the GFE estimator that minimizes the least-squares criterion for all possible groupings of the cross-sectional units. They jointly estimate the individual types and the model's parameter given the number of groups and perform model selection afterward. One the other hand, BLM modify this method and split the procedure into two steps with Kmeans clustering algorithm is used in the first step. From Bayesian's point of view, KW proposes a full Bayesian estimator that simultaneously estimates the group structure and parameters. Unfortunately, none of these works examine the potential forecasting gain when considering the group structure. Our work will fill in this gap and examine the performance of the BGRE estimator in various scenarios.
The paper that most related to ours is liu2020. She considers a linear dynamic panel data model and implements the finite mixture model to estimate the underlying distribution of unit-specific intercepts. However, the assumptions of the underlying model in liu2020 are different from ours. She assumes fully heterogeneous intercepts in a panel data model, which amounts to the standard random-effects model. The finite mixture model in her setting serves as a method to pool information across units. In our work, we specify a group-specific intercept in the population level, and the semiparametric method is the critical ingredient to deliver not only the estimates but also the group structure. Moreover, her main object of interest is to construct individual-specific density forecasts for a panel, while our work includes point, set, and density forecasts and, most importantly, evaluates the performance of group clustering.
This paper also relates to the literature on nonparametric Bayesian approach in the group structure estimation problem. We model the unknown distribution of the heterogeneous coefficient (including grouped intercept and innovation variance) as the Dirichlet Process of Normals with potential infinite groups. The idea of sampling from the Dirichlet Process Model has been widely used by a number of authors including escobar1995, neal2000, ishwaran2001, molitor2010, yau2011, hastie2015, liverani2015, liu2019, and liu2020. To make the infinite-dimensional problem operable, our blocked Gibbs sampler relies on the slice sampling described by walker2007, a more efficient version was later proposed by kalli2011.
We proceed as follows. In section (ref), we present the specification of a linear dynamic panel data model and discuss the construction and evaluation of point, set, and density forecasts. Section (ref) provides details on nonparametric Bayesian priors and subjective group structure priors. It also documents the connection to the Kmeans algorithm. We conduct various Monte Carlo experiments in section (ref) to examine the performance of the proposed estimator in a controlled environment in the light of point, density, and set forecasts. We also examine the performance of a few variants of the BGRE estimator. In section (ref), we conduct empirical analysis in which we forecast the investment rate across firms. Finally, we conclude in section (ref). A description of the data sets, additional empirical results, and derivations are relegated to the appendix.
We consider a panel with observations for cross-sectional units $i=1, \ldots, N$ in periods $t=1, \ldots, T$. Given a panel data set $\{ (y_{it}, x_{it})\}$, a simple linear dynamic panel data model with grouped patterns of heterogeneity takes the following form:
where $x_{i t}$ are a $p \times 1$ vector of exogenous variables, they are uncorrelated with $\varepsilon_{i t}$ but is allowed to be arbitrarily correlated with $\alpha_{g_{i} t}$. $\alpha_{g_{i} t}$ denote the time-varying group-specific heterogeneity. The subscript $g_{i} \in \{1,..., K \}$ is the group membership variable with unknown and unconstrained $K$. $y_{it-1}$ is the lagged outcome variable. $\rho$ is the homogeneous AR(1) parameters that are common for all cross-sectional units, and $\beta_{i}$ is a $p \times 1$ vector of unit-specific slope coefficients. $\varepsilon_{i t}$ is the idiosyncratic error term featured by zero mean and grouped heteroskedasticity $\sigma_{g_{i}}^{2}$, with cross-sectional homoskedasticity being a special case where $\sigma_{g_{i} }^{2}=\sigma^{2}$. This setting leads to a heterogeneous panel with group pattern modeled through $\alpha_{g_{i} t}$ and $\sigma_{g_{i}}^{2}$.
By stacking all observations for unit $i$, we get an aggregated model:
where $y_{i} = \left[y_{i 1}, y_{i 2}, \ldots, y_{i T}\right]^{\prime}$, $y_{i,-1} = \left[y_{i 0}, y_{i 1}, \ldots, y_{i T-1}\right]^{\prime}$, $T \times 1$, $x_{i} = \left[x_{i 1}, x_{i 2}, \ldots, x_{i T}\right]^{\prime}$, $\alpha_{g_{i}} = \left[\alpha_{g_{i} 1}, \alpha_{g_{i} 2}, \ldots, \alpha_{g_{i} T}\right]^{\prime}$, $\varepsilon_{i} = \left[\varepsilon_{i 1}, \varepsilon_{i 2}, \ldots, \varepsilon_{i T}\right]^{\prime}$, $\boldsymbol{\Sigma}_{g_{i}} = \sigma_{g_{i}}^{2} \mathbf{I}_T$. To indicate the component from which each observation stems, we introduce a group membership variable $G = \left[ g_{1}, \ldots, g_{N}\right]$ taking values in $\{1, \ldots, K\}^{N}$. Define a set of unit that belongs to group $k$: $C_k = \left \{i \in \{1,2,...,N\} | g_i = k \right \}$. Let $|C_k|$ denote the cardinality of the set $C_k$.
Following sun2005, lin2012 and BM, we assume that individual group membership does not vary over time. In addition, for any group $i \neq j$, we assume that they have different paths of random effects, e.g., $\alpha_{i} \neq \alpha_{j}$, and no single unit can simultaneously belong to these two groups: $C_{i} \bigcap C_{j} = \emptyset$.
The main goal of this paper is to estimate the grouped random effects $\alpha_{g_i}$, common parameter $\rho$, hetergenous coeffecients $\beta_{i}$ and group membership $G$ using full sample and provide the point, set, and density forecasts of $y_{i t+h}$ for each unit $i$. Throughout this paper, we focus on the one-step ahead forecast where $h=1$. For the multiple-step forecast, the procedure can be extended by iterating $y_{i T+h}$ in accordance with ((ref)) given the estimate of parameters and realizations of data.
Our goal is to generate one-step ahead forecasts of $y_{i, T+1}$ for $i = 1,...,N$ conditional on the history of observations,
and newly available exogenous variables $x_{i T+1}$ at $T+1$. For illustration purpose, we drop $X$ and $x_{i T+1}$ from notations but we always condition on these exogenous variables.
The posterior predictive distribution for unit $i$ is given by
where $\Theta$ is a vector of parameters $\Theta = \{\rho, \beta_i, \alpha_{g_i}, \Sigma_{g_i}, g_i \}$. This density is the posterior expectation of the following function,
which is invariant to relabeling the components of the mixture. Therefore, given $M^*$ posterior draws, the density estimated from the MCMC draws is
Therefore, we can draw samples from $\hat{p}(y_{i T+1} | Y)$ by simulating ((ref)) forward conditional on the posterior draws of $\Theta $ and observations.
We evaluate the point forecasts via the Root Mean Square Forecast Error (RMSFE) under the quadratic loss function averaged across units. Let $\hat{y}_{i T+1} $ represent the predicted value conditional on the observed data up to period $T$, the loss function is written as
where $y_{i, T+1}$ is the realization at $T+1$ and $\hat{\varepsilon}_{iT+1}$ denote the forecast error.
The optimal posterior forecast under quadratic loss function is obtain by minimizing the posterior risk,
This implies optimal posterior forecast is the posterior mean,
We construct set forecasts $CS_{i T+1}$ from the posterior predictive distribution of each unit. In particular, we adopt a Bayesian approach and report the highest posterior density interval (HPDI), which is the narrowest connected interval with coverage probability of $1-\alpha$. Put differently, it requires that the probability of $y_{i T+h} \in CS_{iT+1}$ conditional on having observed the history $Y$ is at least $1-\alpha$, i.e.,
and this interval is the shortest among all possible single connected candidate sets. Let $\delta^{l}$ be the lower bound and $\delta_{u}$ be the upper bound, then $CS_{i T+1} = \left[\delta_{i}^{l}, \delta_{i}^{u}\right]$.
The assessment of set forecasts in simulation studies and empirical applications is based on two metrics: (1) the cross-sectional coverage frequency,
and (2) the average length of the sets $C_{i T+1}$,
To compare the performance of density forecast for various estimators, we examine the continuous ranked probability score (CRPS) across units. The CRPS is frequently used to assess the respective accuracy of two probabilistic forecasting models. It is a quadratic measure of the difference between the forecast cumulative distribution function (CDF), $F^{T+1}_i(y)$, and the empirical CDF of the observation with the formula as follows,
where $y_{i T+h}$ is the realization at $T+1$.
Moveover, we report another metric called the average log predictive scores (LPS) to assess the performance of the density forecast from the view of the probability distribution function (PDF). As suggested in geweke2010, the LPS for a panel reads as,
To evaluate the statistical superiority of pooling within $K$ clusters, we report the bias, standard deviation, average length of 95% credible set, and frequentist coverage of the posterior mean estimate of $\rho$ across Monte Carlo repetitions. For the random effects $\alpha$, we only present the average bias as it may not be of interest for most empirical analysis.
To estimate the number of groups, we derive a point estimator from its posterior distribution, typically, the posterior mean, which is consistent with a quadratic loss function. In the empirical analysis, we also consider the posterior mode suggested by malsiner2016, which is equal to the most frequent number of non-empty components visited during MCMC sampling. These approaches constitute an automatic and straightforward strategy to estimate the unknown number of groups without using model selection criteria or marginal likelihoods.
Until now, our main focus is the group structure for the intercepts $\alpha_i$ while $\rho$ and $\beta_i$ are left unchanged as in a standard panel data model. Here, we can easily extend our model and allow for joint group-specific heterogeneity in $\alpha$, $\beta$ and $\rho$. Then, the extended model is written as,
where $\widecheck{x}_{it} = [1 \; y_{it-1} \; x'_{it}]' $, and $\theta_{g_{i}} = [\alpha_{g_{i}} \; \rho_{g_{i}} \; \beta_{g_{i}}^{\prime}]'$.
From a group $k$, with a joint conjugate prior for parameters $\theta_{k}$, we modify our block Gibbs Sampler to draw $\alpha_k$, $\rho_k$ and $\beta_k$ simultaneously from their joint posterior distribution. The detailed derivation of the posterior distributions are presented in Appendix (ref).
In this section, we provide details in Bayesian analysis. In Section (ref), we document the specification of the prior distribution for all parameters, including the auxiliary variable in the random coefficient model, and the subjective group prior if econometricians have prior knowledge on group structure. Section (ref) outlines the posterior sampler and the proposed algorithm is shown in the Appendix (ref). Finally, we provide preliminary thoughts on the connection between our Bayesian method and unsupervised machine learning methods in Section (ref).
Two sets of specifications of prior distributions are considered in this section. In the first prior specification, we concentrate on random effects models and implement a full Bayesian analysis. In addition, we specify a hyperprior for the distribution of unobserved heterogeneity and then construct a joint posterior for the coefficients of this hyperprior as well as the actual unit-specific and common coefficients. In the second specification, econometricians could provide useful information on the latent group structure and incorporate it in the prior.
In this paper, we focus on the random coefficients model where heterogeneous parameters $\alpha_{g_{i} t}$ and $\sigma_{g_i}$ are independent and are assumed to be independent of the initial value of each unit $y_{i0}$. The specification can be extended to correlated random coefficients model by modeling the joint distribution of heterogeneous parameters and initial values $y_{i0}$.
A typical choice in the nonparametric Bayesian literature is the Dirichlet Process (DP) prior or stick-breaking prior. With group probabilities $\pi_k$ and parameter in prior: (mean, variance) = $(\mu_{\alpha}, \Sigma_{\alpha})$, a draw of $\alpha_{g_{i} t}$ from the DP prior could be viewed as a mixture of point mess with the probability mass function,
where $\delta_{x}$ denotes the Dirac-delta function concentrated at $x$, each $\alpha_{k t}$ is drawn from a normal distribution and $K$ is unknown. $\mu_{\alpha}$ are set to the OLS estimate of $\alpha$ assuming $K = 1$ and $ \Sigma_{\alpha}$ equals $200 \times \hat{\Sigma}_{\alpha}$ where $\hat{\Sigma}_{\alpha}$ is the standard deviation of the OLS estimator. In the same fashion, we can define the DP prior for grouped heteroskedasticity $\sigma^2_{g_{i}}$ given identical group probabilities $\pi_k$:
where each component is drawn from inverse-Gamma distribution.
Put together, the posterior draws of grouped related coefficients can be characterized by a grouped triplet $\{\pi_k, \alpha_{k}, \sigma^2_{k}\}$ for $k = 1,2,..$, and $\alpha_k = [\alpha_{k1}, \alpha_{k2},..., \alpha_{kT}]'$. Importantly, the distributions of both $\alpha_{k}$ and $\sigma^2_{k}$ are discrete, because draws can only take the values in the set $\left\{(\alpha_{k}, \sigma^2_{k}): k \in \mathbb{Z}^{+}\right\}$. This nonparametric nature makes the Dirichlet Process prior an ideal choice for clustering problems especially when the distinct number of clusters is unknown beforehand. The group parameters $(\alpha_{k}, \sigma^2_{k})$ are assumed to follow the base distribution $B_0$ which is an independent (non-conjugate) Multivariate Normal-Inverse-Gamma (IMNIG) distribution.
On the other hand, the group probability is formalized through an infinite-dimensional stick-breaking prior governed by the concentration parameter $a$,
where $\xi_k$, which are called stick lengths, are independent random variables drawn from the beta distribution $Beta(1, a)$. This construction can be viewed as a stick-breaking procedure, where at each step, we independently and randomly break the leftover of a stick of unit length and assign the length of this break to the current value of $\pi_{k}$. The smaller $a$ is, the less of the stick will be left for subsequent values (on average), yielding more concentrated distributions.
The concentration parameter $a$ specifies how strong this discretization is. As $a \to 0$, the realizations are all concentrated at a single value, while when $a \to \infty$, the realizations become continuous-valued from its based distribution. escobar1995 shows that the number of estimated groups under a DP prior is sensitive to $a$, which indicates that a data-driven estimate should be more reasonable. To determine how discrete we want and how many groups are needed given the data, it is convenient to treat $a$ as a parameter under the nonparametric Bayesian framework. Put differently, we can set up a relatively general hyperprior for $a \sim Gamma \left(0.4, 10\right)$, and update it based on the observations. This step generates a posterior estimate of $a$, which implicitly chooses the optimal $K$ without re-estimating the models with different numbers of groups.
Finally, the prior distribution for the common parameter $\rho$ is chosen to be a normal distribution,
The prior of heterogeneous parameter $\beta_i$ follows,
To sum up, in the random coefficients model, we specify the Dirichlet Process priors for group random effects $\alpha_{g_i t}$ and heteroskedasticity $\sigma^2_{g_{i}}$, a stick-breaking process for group probabilities $\pi_k$, a hyperprior for the concentration parameter $a$ and a normal prior for the common parameter $\rho$ and heterogeneous parameter $\beta_i$.
Frequently, researchers could provide a group structure on all or at least part of the units based on personal expertise and the nature of individuals. For example, firms coming from the same industry may share a similar growth pattern with relative high probability; countries having the same level of development form comparable fiscal policies. Though this presumed group structure might be subjective or purely based on theoretical analysis, it is still valuable to integrate this group information as it guides estimation when it enters the algorithm via a prior.
To integrate the prior knowledge with our model, we introduce the prior for the membership probability $\omega_{i} = [\omega_{i1}, \omega_{i2},..., \omega_{iK^p}]'$ for all unit $i = 1,..,N$, where $K^p$ is the preset number of groups. Namely, before estimating the group membership, the researcher assigns each unit to different groups with a set of subjective group-specific probabilities, and these probabilities will enter the algorithm through a prior distribution for $\omega_{i}$. We name this prior for $\omega_{i}$ as Subjective Group Probability (SGP) Prior. In practice, one could provide a table (for example, Table (ref)) documenting the subjective group probability of a unit falling into a specific group.
It is worth noting that such a prior design is rather general and flexible, which is able to account for different scenarios. A special case of the SGP prior is that, for each $i = 1,2,..., N$, the researcher is fairly confident in her knowledge and sets one of $\omega_{ik}$ to 100%. This is equivalent to the case where the researcher exactly partitions $N$ units into $K^p$ predetermined groups.
Building on the prior for the random effects model in Section (ref), we allow for incorporating the researchers' prior knowledge while inheriting the feature of reallocating units and changing the number of groups along the MCMC sampling. These flexible features enable the block Gibbs sampler to automatically correct and update the imprecise subjective prior, especially when $K^p$ doesn't match the true number of groups.
To incorporate these subjective group probabilities, it is important to choose a proper prior for $\omega$. Dirichlet distribution is an applicable candidate among assorted densities since it is the conjugate prior of the multinomial distribution, which facilitates direct sampling, and, most importantly, provides a natural channel to integrate subjective group probability.
To see this, we use the simplest case where the number of potential groups ($K^*$) in an iteration equals the presumed $K^p$. Let $\omega_{i} = [\omega_{i1}, \omega_{i2},..., \omega_{iK^p}]$ be the vector of group-specific probability for unit $i$, we set the prior density for $\omega_{i}$ as an unsymmetric Dirichlet distribution,
where $a_{ik}$ are concentration parameters and strictly positive. Conditional on $\omega_{i}$, the group membership $g_i$ is assumed to be drawn from a multinomial distribution,
It can be shown posterior probability of $\omega_{i}$ given $g_i$ is also a Dirichlet distribution with modified hyperparameters: $Dir (a_{i1} + \mathbf{1}(g_i = 1), a_{i2}+ \mathbf{1}(g_i = 2),..., a_{iK^p}+ \mathbf{1}(g_i = K^p))$. Hence we can direct sample $\omega_{i}$ from it posterior distribution.
Another important property of Dirichlet distribution that enables itself to be the most suitable prior is that we can tie our prior probability directly with the expected value of $\omega_{ik}$,
To integrate the researcher's prior knowledge, one only need to deliberately choose a set of $\{a_{ik}\}$ such that the expected probability matches her subjective probability on groups.
Since the membership probabilities $\omega_{ik}$ are updated based on observations, and we allow for reallocating units and changing the number of groups along the MCMC sampling, a revision of our block Gibbs sampler is needed to adjust for such changes. The details of the new algorithm are presented in Appendix (ref). In practice, we can restrict $\sum_{i = k}^{K^p} a_{ik}$ to be 1 so that $a_{ik}$ represents both the subjective group probability for unit $i$ belonging to group $k$ and the prior mean of $\omega_{ik}$.
Draws from the joint posterior distribution can be obtained by using blocked Gibbs sampling. The proposed algorithm is based on ishwaran2001 and walker2007. Though the algorithm in ishwaran2001 has been widely used for sampling stick-breaking priors, it alone can't fulfill our need for estimating the number of groups without any predetermined level or upper bound since it requires a finite-dimensional prior and truncation. To avoid approximation and a predetermined number of groups $K$, we implement slice-sampling proposed by walker2007 and modify the framework of ishwaran2001 with additional posterior sampling steps. Using the conjugate priors specified in Section (ref), each parameter is directly drawn from its posterior distribution.
In the Appendix (ref), we provide detailed derivations for the conditional distributions over which the Gibbs sampler iterates. We focus on the time-varying grouped random effects model with grouped heteroskedasticity, which is the most sophisticated specification. Other specifications can be estimated by merely ignoring time effects in $\alpha$'s or shutting down the heteroskedasticity.
One feature of our proposed block Gibbs sampler is that it partitions $N$ units into $G$ groups, and, at the same time, generates posterior draws for parameters. This Gibbs sampler and our BGRE estimator inevitably remind us of one of the most popular clustering algorithms in the area of unsupervised machine learning: the Kmeans algorithm. Indeed, the Kmeans algorithm plays a crucial role in BM and BLM, who estimate the grouped fixed-effects from the frequentists' point of view. In this subsection, we seek to illustrate the similarity and connection between our block Gibbs sampler and the Kmeans algorithm in the limit.
We start with the Kmeans clustering algorithm. Given a set of observations $\left(\mathbf{z}_{1}, \mathbf{z}_{2}, \ldots, \mathbf{z}_{N}\right),$ where each observation contains the dependent variables and covariates, $(y_i, x_i')$. Kmeans clustering aims to partition the $N$ observations into $K$ sets so as to minimize the within-cluster sum of squares,
The algorithm alternates between reassigning points to clusters and recomputing the means. For the assignment step, one computes the squared Euclidean distance from each point to each cluster mean, and then assign each observation to the cluster with the nearest mean. The update step of the algorithm recalculates centroid for observations assigned to each cluster and updates $ \boldsymbol{\mu}_{k}$ for all $k$.
According to the block Gibbs sampler, we assign unit $i$ to group $k$ conditional on the draws of other parameters (Eq. ((ref))) with probability,
where $c_{ik} = (2\pi)^{-\frac{T}{2}} \Sigma_{k}^{-\frac{1}{2}} \mathbf{1} (u_i < \pi_k)$, $\tilde{y}_i = y_i - \rho y_{-1,i} - x_i \beta_i $, and $y_{-1,i}$ are the lagged values of $y_i$. If we assume homoskedasticity, i.e., $\Sigma_{k} = \Sigma$ for all $k$, then in the limit as $\Sigma \rightarrow 0,$ the value of $p(g_i = k | \rho, \beta, \alpha, \Sigma, G^{(i)}, Y,X)$ approaches zero for all $k$ except for the one corresponding to the smallest weighted distance $(\tilde{y}_i - \alpha_k)' \Sigma_{k}^{-1} (\tilde{y}_i - \alpha_k)$. In this case, this step is akin to the assignment step of Kmeans but using a weighted Euclidean distance. Then, conditional on newly estimated group membership, we update the group random effects $\alpha_k$ through Bayesian linear regression only using the units of group $k$. This step exactly recalculates the means of the new clusters, establishing the equivalence of the update step.
Having this similarity in mind, it is natural to include the Kmeans algorithm in our Monte Carlo experiment and explore its performance relative to our BGRE estimator in terms of accuracy of clustering. Notably, following BLM, we construct a 2-step GRE estimator equipped with K-mean algorithm in the first step. The performance of this 2-step estimator is assessed in Section (ref).
In this section, we conducted Monte Carlo simulation experiments to examine the performance of various Grouped Random Effects (GRE) estimators under different data generating processes (DGPs) and prior assumptions. These DGPs differ in whether the random effects are time-invariant or time-varying and whether to introduce heterogeneity in the variance of innovations. Such designs allow us to examine not only how our approach performs under DGPs with particular features, but also the reliability of appropriately estimating the number of clusters.
We consider a setting with sample size $N = 100$, and time span $T = 11$. The last observation of each unit forms the hold-out sample for evaluation as we focus on one-step ahead forecasts. A similar framework can be applied to multiple-step-ahead forecasts by iterating from period to period. We set the true numbers of groups $K^0 = 4$. Given $N$ and $K^0$, we partition the entire sample into $K^0$ balanced blocks with $N/K_0$ units in each block\footnote{If $N/K_0$ is not an integer, use $\lfloor N/K_0 \rfloor$ for group 1,2,..,$K_0-1$, the last group contains the residual units.}. For each DGP, 100 datasets are generated, and we run the block Gibbs samplers for each data set with $M = 10,000$ iterations after a burn-in of 5,000 draws.
The Monte Carlo simulation is based on the dynamic panel data model in ((ref)), in which we suppress the exogenous predictors $x_{i t}$ for simplicity. In short, we consider four linear dynamic DGPs in this section. DGP1 and DGP2 involve time-invariant random effects while time-varying random effects are allowed in DGP3. Moreover, DGP1 and DGP3 consider homoskedasticity, but DGP2 has heteroskedastic innovations. DGP4 is the standard panel data model without a group structure. Throughout these DGPs, the random effects $\alpha_{k}$ and idiosyncratic error $\varepsilon_{i t}$ are standard normal distributed, independent across $i$, $k$, and $t$, and mutually independent. $\varepsilon_{i t}$ is independent of all regressors. The data are simulated according to the following DGPs:
{\bf DGP1:} Time-invariant grouped random effects, homoskedasticity {\bf(Grp Ti-Homo)}. \\ This DGP is the most naive panel data model with group pattern in the random effect.
with $\rho = 0.7$, $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, 0.8 \right)$ and $\alpha_k \stackrel{iid}{\sim} N \left(k, 0.5^2 \right)$ for $k = 1,2,...,K^0$.
{\bf DGP2:} Time-invariant grouped random effects, heteroskedasticity {\bf(Grp Ti-Hetero)}. \\ This DGP aims to incorporate heteroskedasticity, which leads to a slightly complicated process that is hard to estimate.
with $\rho = 0.7$, $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, \sigma_{g_{i}}^{2}\right)$ where $\sigma_k^2 = 1.5\left( 1 - \frac{k-1}{K^0} \right)^2$, and $\alpha_k \stackrel{iid}{\sim} N \left(k, 0.5^2 \right)$ for $k = 1,2,...,K^0$.
{\bf DGP3:} Time-varying grouped random effects, homoskedasticity {\bf(Grp Tv-Homo)}. \\ So far, we have focused on time-invariant models. But when estimated on real data, it's reasonable to believe the random effect could a have time pattern. Hence, we introduce various time-varying patterns of the random effects while keeping the assumption of homoskedasticity to avoid over-parameterization.
with $\rho = 0.7$, $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, 1 \right)$ and $\alpha_{it} \stackrel{iid}{\sim} N \left(\underline{\alpha}_{g_i t}, 0.5^2 \right)$ where $\underline{\alpha}_{g_i t}$ varies across periods and groups as depicted in Figure (ref). To enrich the patterns of time-varying random effects, we construct 4 different paths. Group 1 has a constant mean for $\alpha_{i}$. The means for $\alpha_{i}$ in group 2 are also constant but experience a structure change at $T=5$. Group 3 and group 4 are equipped with monotonically increasing/decreasing means.
{\bf DGP4:} Time-invariant random effects, homoskedasticity, no group structure {\bf(Std Ti-Homo)}. This is the standard panel data model with unit-specific random effects and identical variance for the innovations.
with $\rho = 0.7$, $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, 0.8 \right)$ and $\alpha_{i} \stackrel{iid}{\sim} N \left(0, 0.5^2 \right)$.
We consider four types of BGRE estimators that differ on the assumptions made on the random effects (RE) and the variance of errors: (1) time-invariant grouped RE with homoskedasticity (Ti-Homo); (2) time-invariant grouped RE with heteroskedasticity (Ti-Hetero); (3) time-varying grouped RE with homoskedasticity (Tv-Homo); (4) time-varying grouped RE with heteroskedasticity (Tv-Hetero)\footnote{For the Tv-Homo and Tv-Hetero estimator, as we allow for time effects in $\alpha_i$, we use the most recent $\alpha_{iT}$ to make one-step ahead prediciton. This is equivalent to assume the law of motion of $\alpha_{it}$ is a random walk. Modeling the trend of $\alpha_{it}$ would result in a more accurate forecast, but this is beyond the Scope of this paper.}. For instance, Ti-Homo estimator assumes the true model to have time-invariant grouped random effects, and the variance of error terms to be constant across units. Besides the results we will show below, we modify the DGPs and conduct more experiments using (a) larger variance of error terms $\sigma^2$, (b) shorter time span, and (c) different true number of groups $K^0$. These additional results are available in Appendix (ref).
Regarding alternative estimators, we consider the following Bayesian estimators that have different prior assumptions on the random effects $\alpha_i$. (1) Bayesian pooled estimator (Pooled): $\alpha_i$ is treated as a common parameter as $\rho$ does, this means all units share the same prior level of $\alpha_i$; (2) flat-prior estimator (Flat): assume $p(\alpha_i) \; \propto \; 1$, this amounts to draw samples from a posterior whose mode is the MLE estimate. Given the estimate of common parameter, there is no pooling across units, $\alpha_i$'s are estimated only using their own history; (3) Parametric-prior estimator (Param): assume $\alpha_i \sim N(\mu, \pi^2)$, where a Normal-Inverse-Gamma hyperprior is further imposed on $(\mu, \pi^{2})$\footnote{The Normal-Inverse-Gamma hyperprior for $(\mu, \pi^{2})$ used in the Monte Carlo simulation is as follow: $\mu | \pi^{2} \sim N(m, v\pi^{2})$ with $m$ equating to the pooled OLS estimator of $\alpha_i$ and $v = 1$; $\pi^{2}$ follows $IG(\nu_{\pi}/2, \delta_{\pi}/2)$ with $\nu_{\pi} = 6$ and $\delta_{\pi} = 4$.}, this prior be thought of as a limit case of the DP prior when the concentration parameter $a \rightarrow 0$, so there is only one cluster, and $\left(\mu, \pi^{2}\right)$ are directly drawn from the base distribution.
Table (ref) shows the estimate comparison among alternative predictors. For the DGP1, Ti-Homo and Ti-Hetero estimator are the best in every aspect. This is as expected since they correctly model the time-invariant random effects. Among these two estimators, when we allow for group‐level heteroskedasticity, the optimal number of groups decreases as Ti-Hetero underestimates the number of groups. The coverage probability, however, is not well-controlled, both of which are below the nominal coverage of 0.95. The Flat estimator also has good performance in terms of RMSE of $\rho$. Nonetheless, its coverage probability is relatively low: only 23% of credible sets successfully contain the true values. This is due to the relatively large biases for $\alpha_{i}$. The rest predictors are considerably worse. This implies that completely ignoring group structure (Pooled, Param) results in a significantly inferior estimate, so does wrongly modeling time-varying random effects (Tv-Homo, Tv-Hetero). Regarding the performance of clustering, Ti-Homo slightly underestimates the number of groups with an average $K$ equals to 3.60 while the truth is 4.
For the case of DGP2, we keep time-invariant random effects while assuming heteroskedasticity. Tv-Homo and Tv-Hetero are still dominating. Tv-Hetero generates the best results with an accurate estimate of the number of groups as it correctly models heteroskedasticity which in turn improves the estimation efficiency. The Flat estimator closely follows them, and the rest are worse. Regarding DGP3, when time-varying random effects are introduced in the model, Tv-Homo and Tv-Hetero estimator yield the best performance and the estimated average $K$s are close to the truth. But the biases are arguably low for these two estimators in sacrifice for small standard deviation and short credible intervals. It is worth noting that, unlike Ti-Homo in DGP1 and Ti-Hetero in DGP2, though correctly specified, the bias for $\alpha_{i}$ is still comparatively high. This is because, for simplicity, we don't model the law of motion for $\alpha_{it}$ and simply assume $\alpha_{iT+1} = \alpha_{iT}$, which results in large bias in $\alpha_i$. As regards the DGP4 that doesn't have a group structure, the Flat estimator is the best since it doesn't pool cross-sectional information but estimate the unit-specific random effects. All of the BGRE estimators have almost identical performances with the estimated number of groups close or equal to 1. Since the Pooled and Param estimators assume no group structure, both have similar estimates as the BGRE estimators.
Table (ref) reports the predictive performance of a range of parametric forecasts. For the DGP1, the best forecasts are generated by the Ti-Homo estimator, as it is correctly specified in this environment. It has the smallest RMSFE, the shortest average length of the credible set, correct coverage probability, the largest LPS, and the smallest CRPS. Although allowing for heteroskedasticity along with the time-invariant random effects, Ti-Hetero generates a perfect point forecast as well as Ti-Homo. But Ti-Hetero introduces uncertainty revealed by a slightly wider credible set and worse density forecast. Moreover, estimators involving time-varying random effects (Tv-Homo and Tv-Hetero) worsen the forecast. Finally, incorrectly imposing no latent group pattern substantially deteriorates the predictive performance in all aspects.
For the DGP2, Ti-Hetero is expected to have an edge, and indeed, it dominates the remaining alternatives. Ti-Homo performs as great as Ti-Hetero in point forecast. This is because, from the previous section, we know that Ti-Homo generates accurate point estimates for both $\rho$ and $\alpha_{i}$. But Ti-Homo fails in the set forecast and density forecast, which illustrates the importance of modeling heteroskedasticity. Again, the rest estimators suffer apart from the Flat estimator.
In the DGP3, Ti-Homo and Ti-Hetero are doing badly by not capturing the time effects in $\alpha_{g_{i}}$. Tv-Homo and Tv-Hetero are the best, beating the rest by a large margin, and equally accurate in this setup. The coverage probability for these two estimators is slightly lower than that of Ti-Homo and Ti-Hetero in part due to uncertainty introduced by more parameters of interest. Pooled, Flat, and Param estimator neglect both group structure and time-varying random effects, hence generating poor forecasts.
Regarding DGP4, all estimators beside Param have comparable forecasts as all of them deem no group structure in this environment (all estimated $K$s are close to 1). However, Param suffers from high variance, and it always generates the widest credible interval and worse density forecast as in other DGP's.
This section conducts four sets of Monte Carlo simulation experiments aiming to examine the variants of BGRE estimator: (1) the GFE estimator proposed by BM, (2) a two-step GRE estimator with Kmeans, (3) the BGRE estimator with the subjective group prior, and (4) the BGRE estimator imposed with true $K^0$. The main text ignores part of the estimators we consider in the previous section and focuses on the correctly specified estimator for each DGP.
In addition to four DGPs specified in Section (ref), we design three new DGPs that inherit the main features from DGP1, DGP2, and DGP 3, including balanced group structure and unit variance structure. However, instead of drawing $\alpha_{g_i}$ from the normal, we choose to use a constant $\alpha_k$ for each group. In this way, we could focus on the clustering result rather than repetitions to average out the randomness brought by random effects. We impose $K^0 = 4$ and use this number throughout the Gibbs sampling.
{\bf DGP5:} Time-invariant grouped fixed effects, homoskedasticity.
with $\rho = 0.7$, $\alpha_k = k$ for $k = 1,2,...,K^0$ and $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, 1 \right)$.
{\bf DGP6:} Time-invariant grouped fixed effects, heteroskedasticity.
with $\rho = 0.7$, $\alpha_k = k$ for $k = 1,2,...,K^0$ and $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, \sigma_{g_{i}}^{2}\right)$ with $\sigma_k^2 = 1.5\left( 1 - \frac{k-1}{K^0} \right)^2$.
{\bf DGP7:} Time-varying grouped fixed effects, homoskedasticity.
with $\rho = 0.7$, $\varepsilon_{i t} \stackrel{iid}{\sim} N \left(0, 1 \right)$ and $\alpha_{g_i t} = \underline{\alpha}_{g_i t}$ where $\underline{\alpha}_{g_i t}$ varies across periods and groups as depicted in Figure (ref).
In this experiment, we compare our BGRE estimator with BM's GFE estimator . In particular, we assess the performance of point forecast in $y$ and the accuracy of coefficient estimates (group random effects $\alpha$ and common parameter $\rho$).
To compare the performance of estimators, we use the default numerical setting\footnote{The default settings are as follow: (1) Number of groups = 4; Number of covariates = 1; Standard errors: 0 (no standard errors). (2) For algorithm 0, Number of simulations = 100; (3) For algorithm 1, Number of simulations = 10, Number of neighbors = 10 , Number of steps = 10.} in BM. It is worth noticing that BM relies on information criteria to ex post select the optimal number of groups. Hence we consider at most 10 groups and estimate the number of groups $K$ in accordance with the following Akaike information criterion (AIC)\footnote{We also tried the alternative choice $\hat{\sigma}^{2} \frac{k T + N + p + 1}{N T} \ln (N T)$ for the penalty. This corresponds to the default BIC used in bonhomme2015. We found that, in this case, BIC selected the smallest possible number of groups for all DGPs, i.e., no group structure, whereas the truth is $K^{0}=4$. Moreover, other forms of BIC could always select the largest $K$ as well. Due to the inaccurate estimate of group structure and substantially poor performance, we don't show the result with the default BIC.}:
where $\widehat{\sigma}^{2}$ is a consistent estimate of the variance of $\varepsilon_{i t}$:
The results are shown in the Table (ref). As BM proposes two algorithms\footnote{These two algorithms could generate different estimates as shown in the case of DGP6 below.}, we present the results for four versions of the GFE estimator. The first two estimators ($\text{GFE}_{a0}$ and $\text{GFE}_{a1}$) equip with the Iterative and Variable Neighborhood Search algorithm, respectively. We impose the true number of groups $K^0$ for the other two estimators ($\text{GFE}^0_{a0}$ and $\text{GFE}^0_{a1}$), i.e., we don't perform model selection but choose $\hat{K} = K^0 = 4$ directly.
For the DGP5 and DGP6, Ti-Homo and Ti-Hetero estimator outperform the GFE estimators in all aspects. The GFE estimator's poor performance is mainly due to the incorrect estimate of the number of groups. This also emphasizes that an inaccurate estimate of group structure would deteriorate both estimates and point forecasts. It is worthy noting that, even imposing the true number of groups, GFE estimators ($\text{GFE}^0_{a0}$ and $\text{GFE}^0_{a1}$) are still straggling and worse than their counterparts with $\hat{K}$ selected by the information criterion. $\text{GFE}^0_{a0}$ and $\text{GFE}^0_{a1}$ generate relatively high bias for both $\alpha_{i}$ and $\rho$. In the case of DGP7, Tv-Homo and Tv-Hetero still have better performance than the GFE estimators whose models are chosen by the information criterion. Furthermore, Tv-Homo and Tv-Hetero perform only marginally worse than $\text{GFE}^0_{a0}$ and $\text{GFE}^0_{a1}$ estimators that are imposed with the true $K^0$. This is because they overestimate the number of groups in some posterior draws, and hence, on average, the posterior mean forecast and estimate are slightly off.
Moreover, the accuracy of the GFE estimator is profoundly affected by the choice of the information criteria. We implement several information criteria proposed by bai2002 in this Monte Carlo experiment and find that there is no single criterion that consistently selects the correct number of groups nor close to the truth. As the GFE estimator is designed for the time-varying model, we finally select the AIC mentioned above, which chooses a model that is closest to the true model in time-varying DGPs. But this deliberately selected AIC fails to improve the performance of the GFE estimator in DGP7. Once we switch to other forms of AIC or BIC, these results will not hold anymore. These facts also emphasize the importance of not relying on ex-post model selection and the superiority of our Bayesian estimators.
In this section, we explore whether subjective group structure improves the accuracy of forecast and group clustering. We conduct two Monte Carlo experiments corresponding to SGP Prior defined in Section (ref).
We consider five scenarios, each of which differs in the structure of subjective prior probability and hence the prior group probability $\pi$. The exact specification is characterized in Table (ref). The first three scenarios set the preset number of groups as the truth $K^0$, whereas subjective group probabilities are different in levels. In the scenario 1, the researcher is pretty confident about her decision and assigns 100% to the right group for each unit, which amounts to knowing the true membership. However, she is entirely uninformed and cluster units with even probability for each group (i.e., $\omega_{ik} = 1/K^0$ for $\forall i,k$) in the scenario 3. Scenario 2 is an intermediate case where one is less confident in her knowledge and correctly assigns a unit to its group with the prior probability of 70% (other groups equally split the remaining 30%). For the scenario 4 and 5, the number of groups is different from the truth. We assume the researcher divides all units into $K \ne K^0$ even groups with the prior probability of 100%\footnote{The last two scenarios aim to evaluate the performance of SGP prior when the number of groups is wrong. Instead of randomly assigning a unit to each groups with a set of probability, we assume the econometrician is confident on her prior and set 100% for a particular group. In particular, we assign the first $N/K$ units into group 1, the next $N/K$ units into group 2, and so on. We also run other designs for scenario 4 and 5 with different prior probabilities. The results show that, as long as a certain amount of units are correctly clustered into groups, the performances of SPG-RE estimators are slightly better than those of the BGRE estimator.}.
We conduct the simulation experiment for the SGP Prior, where we examine the impact of subjective group probability prior on the performance of estimate and forecast under DGP3 (time-varying random effects model). To adopt the SGP prior, we revise the BGRE estimator and construct the new Bayesian estimator under the assumptions: (1) time-varying random effects and (2) heteroskedasticity, which corresponds to the most general case of a panel data model. We name it as Subjective Group Probability Random Effect (SGP-RE) estimator. Given different specifications of SGP prior, we report the performance of five SGP-RE estimators relative to the Tv-Hetero estimator (benchmark model).
Figure (ref) and (ref) depict the relative performance of various SGP-RE estimators against the Tv-Hetero estimator\footnote{The full results are presented in the Appendix (ref).}. Remember that the prior knowledge is the most accurate in scenario 1, where the researcher is equivalent to know the true group structure. In this regard, a clear gain in estimate emerges as the RMSE for $\rho$, bias for $\rho$ and $\alpha_i$ generated by SGP-RE1 (100% confidence) decreases by more than 35% relative to the Tv-Hetero estimator. As we move from scenario 1 to the other scenarios, the prior information becomes less accurate. SGP-RE2 (70% confidence) beats the benchmark with moderate improvement on the RMSE ($\sim$22%) and the bias ($\sim$32%) while SGP-RE3 (uninformed) underperforms the benchmark by a huge margin. The poor performance of SGP-RE3 is not surprising. Though it correctly specifies $K = K^0$, the uninformative prior forces the algorithm to consider other incorrect groups with a large chance (= $1 - \frac{1}{K^0}$), and hence deteriorates the performance of both estimates and forecasts. In terms of the one-step ahead forecast, SGP-RE1 leads the rest by scoring the lowest values for each metrics (and the highest LPS), closely followed by SGP-RE2. SGP-RE3 is suffering in terms of point and set forecast. As for scenario 4 and 5, despite the flexibility of allowing for changing $a$ and $K^p$ along MCMC sampling, the incorrect specifications of the group number deteriorate the estimation, both of which fail to deliver reliable estimates for $\rho$ and $\alpha$. Nonetheless, such prior structures help point and density forecasts. Both estimators beat the benchmark and generate comparable LPS to SGP-RE1 and SGP-RE2. This valuable improvement mainly results from the fact the algorithm could exploit the prior knowledge on group structure that successfully partitions merely a fraction of units and adapt the number of groups accordingly.
In short, regarding the overall performance of the SGP prior, the best case is that the researcher has a relative confident prior and knows the true number of groups. In this case, the SGP-RE estimator dominates the Tv-Hetero estimator from every angle. However, in practice, we rarely come up with such a precise prior due to the incomplete understanding of the population of the data. Instead, we might be less confident on our knowledge or even specify more/fewer groups than the truth. Under this circumstance, the SGP-RE estimator could still deliver a better density forecast because of the great exploration of the prior information and the adaptive scheme featured by our Bayesian method.
In this section, we compare our BGRE estimator with a two-step GRE estimator, where units are clustered into groups in the first step using Kmeans algorithm, and the model is then estimated in the second step with group-specific heterogeneity. Unlike BLM, we implement the Bayesian framework in the second step to echo other Bayesian estimators presented in the previous section. This two-step procedure allows us to examine the clustering accuracy of Kmeans relative to our full Bayesian estimate as two-step GRE estimators can be viewed as a special case of the GRE estimator with group membership determined by Kmeans. The optimal number of clusters in Kmeans is selected by the average silhouette method.
To avoid cluttering the tables in the main text, we depict the selected results\footnote{The full results are presented in the Appendix (ref).} for estimates and forecast in Figure (ref) and (ref). Each bar represents the performance of the two-Step GRE estimators against the performance of the original GRE estimator. Above zero indicates the 2-step estimator underperforms the benchmark while a 2-step estimator is better when its bars show negative values. The main models are correctly specified for each DGP, i.e., Ti-Homo for DGP 1, Ti-Hetero for DGP 2, and Tv-Homo for DGP 3.
Figure (ref) presents the point estimates for each DGP. We document the root mean squared forecast error, absolute bias, standard deviation, and the average length of the 95% credible set for the common parameter $\rho$, while the metric for the random effects is the absolute bias. According to these measures, the two-Step GRE estimators perform worse than the Bayesian GRE estimator as they introduce much higher bias in the estimate of $\rho$ (hence larger RMSE for $\rho$) and $\alpha_i$. It is worth noting that the estimator equipped with Kmeans doesn't affect the standard deviation and the average length of 95% credible set of $\rho$.
The inferior performance of the 2-step GRE estimator is due to the inaccurate estimate of group structure. Table (ref) reports the estimated number of groups from the two-step GRE estimators with Kmeans and the BGRE estimator. Regarding the performance of clustering, the Kmeans algorithm severely underestimates the number of groups as it prefers much less groups, while the true number is 4. Meanwhile, our BGRE approach accurately estimates the number of groups, though slightly underestimated in the DGP 1.
Figure (ref) shows the point, set, and density forecast for each DGP. As Kmeans fails to estimate the group structure, none of the 2-step GRE estimators outperform the GRE estimators. Namely, Kmeans clustering doesn't help make a more accurate forecast, and instead, it generates a much higher forecast bias and standard deviation.
In this section, we illustrate the use of Bayesian Grouped Random Effects estimators in a cross-firm study. We revisit the investment regression and use a different version of the dynamic grouped panel model to forecast the investment rate\footnote{The investment rate for a firm in a particular year is defined as the fraction of capital expenditures in property, plant, and equipment in terms of the beginning-of-year capital stock.} for a panel of firms in all industries. Instead of using the traditional Tobin's Q-type investment regression, we implement a new scheme proposed by gala2019, who directly estimates the corporate investment rate without Tobin's Q. Again, our main focus is the one-step ahead point, set and density forecast. Due to space limitations, we only report forecast results for the most recent year in the main text. Summary statistics, and additional details of implementation are stored in Appendix (ref) and (ref).
We consider a general model with grouped latent heterogeneity in $\alpha_{i}$. Following hsiao1997 and gala2019, the investment equation is specified as,
where capital stock, $K_{it}$ is defined as net property, plant and equipment; $I_{it}$ is capital investment; $CF_{it}$, is a liquidity variable defined as cash flow minus dividends; $Y_{t}$ is the end-of-year sales; $\varepsilon_{it}$ are the normally distributed error terms. The subscript $i$ denotes companies, and $t$ denotes time. Unlike the commonly specified investment equation using Tobin's Q, the additional terms, including the natural logarithm of lagged capital and sales-to-capital ratio, are based on the regression proposed by gala2019. The lagged values of the investment rate are included as explanatory variables to avoid endogeneity problems.
As we focus on forecasting, we can relax a few assumptions to achieve better predictive performance. These assumptions include time-invariant random effects $\alpha_{g_i}$, homoskedasticity in $\sigma_{i}$ and homogeneous coefficients for all dependent variables ($\beta_{i \cdot}$ = $\beta_{i}$). Table (ref) summarizes the estimators and their properties we consider in this section. The implementation of time-invariant RE and homoskedasticity is similar to the one in the previous section, i.e., construct four versions of BGRE estimator: Ti-Homo, Ti-Hetero, Tv-Homo, and Tv-Hetero. Despite the fact that the homogeneous slopes have been frequently rejected in empirical studies of estimates and inference, such an assumption could provide potential improvement in forecasts.
The individual company data are obtain from COMPUSTAT Annual database. To account for potential structural breaks and the advanced speed of capital accumulation in the recent decades, our sample is composed of a balanced panel of firms for the years 2000 to 2019, that includes firms from all industries with no missing value in accounting data.
We keep only firm-years that have non-missing information required to construct the primary variables of interest, namely: capital stock $K$, investment $I$, liquidity $CF$, and sales revenues $Y$. The further details of constructing the sample can be found in the Appendix (ref). The final sample comprises 337 firms and the observations for each firm is 20.
To examine the performance of various estimators with limited observations, we choose to use a rolling window of 15 years. In this sense, we create five balanced panels which end in years 2014, ..., 2018 ($t=T$), respectively. The observations in the next year ($t=T+1$) are reserved for pseudo-out-of-sample forecast evaluation. For illustrative purposes, we will present the results for the year 2019 in the remainder of this section. The full results are presents in Appendix (ref).
We begin the empirical analysis by comparing the performance of point, set, and density forecast for the last panel (in-sample periord: 2005 - 2018). We aim to forecast the investment rate in 2019. We consider all the model specifications depicted in the Table (ref) and their performance is presented in the Table (ref). Throughout the analysis, the Flat estimator serves as the benchmark as it essentially assumes individual effects. In Table (ref), the third column shows the RMSFE for the one-step ahead forecast. For the panel considered in the table, we first notice that the best model is the Tv-Hetero (time-varying random effects, heteroskedasticity) in homogeneous coefficients specification. It outperforms the benchmark -- Flat estimator by 25%. Ti-Hetero also delivers accurate point forecast, which suggests time effects provide merely marginal improvement. Under heterogeneous coefficients specification, for the BGRE estimators, though all of them beat the Flat estimator, their RMSFEs are relatively larger. This may arise from the fact that heteroskedasticity alone can capture a great amount of individual effects, thus imposing heterogeneous coefficients in $\beta_{i}$ may overfit the model and lead to poor forecasts.
The fourth column documents the average number of latent groups in $\alpha_{i}$. Most of our BGRE estimators deem a group structure with more than six underlying components. And as we will discuss later, this rich group structure is the cornerstone of the accurate and flexible density forecast. In Figure (ref), we present the posterior distribution of the number of groups for those BGRE estimators that have more than two groups. Regardless of the predictive performance, most estimators agree on the number of groups, with the posterior mode ranging from 6 to 7. On the other hand, the SIC code, which categorizes companies into the industries by their business activities, suggests that there are ten different industries in our sample. This indicates that our block Gibbs sampling algorithm reshuffles the default group structure and optimally pools firms from several sectors.
The fifth and sixth columns present the average coverage rate with the nominal coverage probability of 95% and the average length of the 95% credible set. In general, the coverage rates generated by the homogeneous coefficients specifications are substantially larger than the sets obtained from the models with heterogeneous coefficients and are closer to the nominal coverage probability of 95%. However, the homogeneous setting doesn't improve the average length of the credible set. Indeed, the decrease in the average length is evident for the heterogeneous coefficient models, and it becomes even more pronounced once we impose heteroskedasticity. This is because the sizeable cross-sectional variation in the posterior predictive distributions leaves plenty of room for heterogeneous and heteroskedastic models to shorten the credible set while maintaining the coverage rate in a reasonable range.
The last two columns in Table (ref) depict the performance of the density forecast. Consistent with the point forecast, the Ti-Hetero and Tv-Hetero models under homogeneous coefficients specification have comparable performance and dominate the rest with larger LPS and smaller CRPS. The Ti-Hetero has the largest LPS while Tv-Hetero scores the lowest CRPS. These facts emphasize that incorporating homogeneous coefficients and heteroskedasticity is crucial for density forecasting while time effects are not important.
The results for estimation and forecast in other years are presented in Appendix (ref). In short, the result for 2019 is representative, as most conclusions discussed above also apply for other years. Although no single estimator consistently dominates the rest across the years, at least one of our BGRE estimators always offers the best performance and beats the standard panel data models.
To further investigate the posterior predictive density, we plot the densities of the investment rate generated by the Tv-Hetero and Flat estimator under homogeneous specification in Figure (ref). Both posterior predictive densities have similar posterior means while imposing the grouped random effects remarkably sharpens the density around the mean. The reason is that Tv-Hetero estimator leverages latent group structure and pools the information of firms that share great similarities while the Flat estimator treats individual firm separately and makes a prediction based on limited observations.
Figure (ref) further aggregates the predictive density over industries. Comparing Tv-Hetero and Flat estimator across industries, several observations stand out. First, Tv-Hetero predictive densities tend to be more concentrated in each industry. This is in line with Figure (ref). Second, there is substantial heterogeneity in density forecasts across sectors. While pooling the forecast for all firms yielding a well-behaved bell shape, the posterior predictive densities for the individual industries are in various shapes. This would pose difficulties for the standard estimator to forecast the future and call for a flexible model specification. Third, the Flat estimator is not flexible enough to portray the potential non-normal predictive distribution due to the restrictive normality assumption. In this case, our BGRE estimator, especially the Tv-Hetero estimator, manages to depict the skewed trimodal distribution for the Retail Trade division and bimodality in the Finance, Insurance and Real Estate divison via combining different groups.
This paper studies the estimation and prediction of a dynamic panel data model with latent grouped random effects. We adopt a nonparametric Bayesian approach to identify coefficients and group membership in the random effects simultaneously. This approach avoids the severe issue introduced by the ex-post model selection and allows us to incorporate any forms of prior knowledge on group structure. In Monte Carlo experiments, we show that the BGRE estimators have the edge over standard Bayesian estimators. Regarding clustering, the BGRE estimators generate comparable performance with the Kmeans algorithm. Our empirical application to investment rates across firms reveals that the estimated latent group structure provides a great amount of flexibility and improves point, set, and density forecasts.
The present work raises interesting issues for further research. First, it may be appealing to consider group structures in the AR(1) parameters and heterogeneous coefficients. This would allow us further to reduce the complexity of a panel data model and may improve predictive performance. Second, more clever attempts could be made to incorporate the subjective prior group structure. Our proposed method summarizes prior information in the prior of the membership probability, which can be further improved to overcome its limitation. Third, our method can be extended to nonstationary panels, where panel units and co-integrating relationships may possess latent group structures. Four, the assumption that an individual cannot change its group identity during the whole sampling period can be relaxed in the next step, leading to an even more flexible specification.
\setstretch{1}
\setstretch{1.3}