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.
116,895 characters · 20 sections · 0 citation commands
A Semi-Parametric Bayesian Generalized Least Squares Estimator
The Generalized Least Square (gls) estimator is a family of econometric methods that have seen numerous applications in empirical economics. As pointed out by \hyperlink{Wooldridge, 2003}{Wooldridge (2003)}, parametric gls type estimators accommodate a deviation from the assumption that the errors in the model are homoskedastic and have no serial correlation. For example, compared to the ordinary least squares regression model, gls no longer assumes that the covariance matrix of the errors is diagonal with identical diagonal elements.
In a more general setting, each error may be a random vector, which includes some of the most popular applications of gls. For example, the Seemingly Unrelated Regression (sur, \hyperlink{Zellner, 1962}{Zellner, 1962} and \hyperlink{Zellner, 1971}{1971}) has been developed for equation systems, and widely applied to gain efficiency by exploiting the correlation between errors across equations. Similarly, in the analysis of panel data the random effects model (rem) recognizes that there are individual specific, time-invariant features that are unobservable and uncorrelated with the explanatory variables.
However, the parametric gls still maintains the assumption that the error vector of each individual has the same covariance matrix. In reality, however, heterogeneity in error distributions is a major concern in empirical analyses. Such heterogeneity can be caused by observations on individuals or households reflecting variation in features such as the size of the household and the level of income, among others. It is a challenge for analysts who seek efficient estimates and inference with the data to capture the form of the heterogeneity in observations.
The standard Bayesian approach to gls assumes that the error distribution is multivariate normal. Recent developments in Bayesian methods allow the use of prior information to relax this assumption. The Dirichlet prior has been introduced to accommodate heterogeneity in the distributions of both errors (see \hyperlink{Chigira and Shiba, 2015}{Chigira and Shiba, 2015} for an example) and model parameters (\hyperlink{Allenby et al., 1998}{Allenby et al., 1998}) by mixing a fixed number of normal distributions.
A notable drawback of the Dirichlet prior is that the the dimension of the mixing distribution is usually unknown. Bayesian semi-parametric methods introduce more flexibility by letting the data and the prior determine the structure of heterogeneity jointly. The Dirichlet Process ($\mathcal{DP}$) prior\footnote{See \hyperlink{Escobar and West, 1995}{Escobar and West, 1995} and \hyperlink{Escobar and West, 1998}{1998} and \hyperlink{MacEachern, 1998}{MacEachern, 1998} for a reference of the Dirichlet Process prior.} can be used to form a mixing of normal distributions, whose dimension need not be predetermined. In this sense, the use of $\mathcal{DP}$ priors represents a more flexible approach to accommodating heterogeneity than the mixing of a fixed number of normal distributions with the Dirichlet prior.
In the context where heterogeneity in the distributions of errors is a major concern, $\mathcal{DP}$ priors are introduced for the distributional parameters of the errors in the model\footnote{E.g. for each error vector to have its distinct normal distribution, its distributional parameters are a mean vector and a covariance matrix.}. Under this prior, the distributional parameters are put into groups, and assigned a group specific value. As a result, the corresponding errors will have group specific distributions.
A landmark study in this area is \hyperlink{Conley et al., 2008}{Conley et al. (2008)}, who introduce a Bayesian semi-parametric approach to the instrumental variable problem in a two stage least square type model. Due to the endogeneity of some explanatory variables, the errors in the two stages are correlated by construction. Instead of assuming that the joint errors in the two stages have an identical bivariate normal distribution (c.f., \hyperlink{Chao and Phillips, 1998}{Chao and Phillips, 1998}; \hyperlink{Geweke, 1996}{Geweke, 1996}; \hyperlink{Kleibergen and van Dijk, 1998}{Kleibergen and van Dijk, 1998}; \hyperlink{Rossi et al., 2005}{Rossi et al., 2005}), the authors introduce a Dirichlet process prior for the distributional parameters. This provides a semi-parametric version of two stage least squares, where the errors of the two stages jointly follow a non-parametric mixture of normal distributions.
In this paper we focus on relaxing the identical distribution assumption on the errors, but in a different scenario from that of \hyperlink{Conley et al., 2008}{Conley et al. (2008)}. We propose a semi-parametric Bayesian gls that incorporates the $\mathcal{DP}$ prior. The motivation is to incorporate more information in the error distribution by allowing their distributional parameters to differ across observations. The resulting distribution of the error terms will involve a mixture of normal distributions where the number of the normal components is influenced by both the prior and the data. We then introduce two specific cases of semi-parametric Bayesian gls, namely for equation systems and panel data.
The rest of the paper is organized as follow. In Section (ref) we briefly review the literature on the Dirichlet process and its application as a prior for semi-parametric Bayesian estimators. We then consider two special cases of the gls. The sur with the $\mathcal{DP}$ prior (dp-sur) is introduced in Section (ref) with simulation and empirical results. Section (ref) motivates and introduces our semi-parametric Bayesian gls estimator for the random effects model with the $\mathcal{DP}$ prior (dp-rem) for panel data, together with simulation and empirical results. Section (ref) concludes the paper.
In this section we introduce the generic form of Bayesian gls with the $\mathcal{DP}$ prior. In Section (ref) we briefly review the literature in related areas. In Section (ref) we introduce the semi-parametric Bayesian gls.
Bayesian attempts to incorporate heterogeneity in the distributional parameters of the errors in the linear regression model can be traced back to \hyperlink{Geweke, 1993}{Geweke (1993)} with the use of an inverse gamma prior for the variances of the errors. He demonstrates that such a scale mixture of normal distributions is equivalent to the errors having a t-distribution.
Although the model with the t-distributed errors is more flexible than assuming a normal distribution for the errors, this approach depends upon the assumption that the normal distributions are mixed with inverse gamma distributed variances. As pointed out by \hyperlink{Koop, 2003}{Koop (2003)}, relaxing this assumption results in more flexible models, given that the errors are no longer restricted to having a t-distribution. This can be done by using a Dirichlet prior, the conjugate prior of a multinomial distribution, to mix a finite number of normal distributions.
The Dirichlet mixture model has emerged as a widely applied methodology for capturing heterogeneity in both linear and non-linear models, c.f., \hyperlink{Allenby et al., 1998}{Allenby et al. (1998)}, \hyperlink{Li and Tobias, 2011}{Li and Tobias (2011)} and \hyperlink{Chigira and Shiba, 2015}{Chigira and Shiba (2015)}. The main limitation is that it takes a fairly difficult test procedure to determine the “correct” number of mixing components.
In the wake of this limitation of the Dirichlet mixture model, it seems more reasonable to let the data and the prior jointly determine the number of normal components in the mixture. This can be achieved using a Dirichlet prior of infinite dimension, which is the Dirichlet process ($\mathcal{DP}$) introduced by \hyperlink{Ferguson, 1973}{Ferguson (1973)}\footnote{See \hyperlink{Teh, 2011}{Teh (2011)} and \hyperlink{Gershman and Blei, 2012}{Gershman and Blei (2012)} for reviews of the Dirichlet process.}. $\mathcal{DP}$ is the conjugate prior for a non-parametric multinomial distribution of infinite dimensions. The generic form of the $\mathcal{DP}$ can be written as
where $\alpha>0$ is the concentration parameter, and $F_0$ is the base distribution. $F$ is a random distribution that is discrete with probability one\footnote{The level of discreteness is influenced by $\alpha$, the concentration parameter.}.
The $\mathcal{DP}$ is a non-parametric “distribution of distributions” (\hyperlink{Escobar and West, 1995}{Escobar and West, 1995} and \hyperlink{Escobar and West, 1998}{1998}; \hyperlink{MacEachern, 1998}{MacEachern, 1998}), in the sense that a draw, say $F$, from a $\mathcal{DP}$ is a probability distribution itself. Conditional on $n-1$ existing realizations $\{ r_1, r_2, \ldots, r_{n-1} \}$ from $F$, the Chinese restaurant process (\hyperlink{Aldous, 1985}{Aldous, 1985}) provides the predictive probabilities of the $n^{th}$ realization, $r_n$. Due to the fact that $F$ is discrete, the existing realizations will be assigned to groups, where all realizations in the same group take a group specific unique value.
Denote the group id of $r_i$ as $c_i = 1, \ldots, K$, and the unique value of group $c_i$ as $r_{c_i}^*$: if $r_i$ is in group $k$, then $r_i = r_{c_i}^* = r_{k}^*$. The prediction probabilities of $r_n$ is given by
where $n_k$ is the number of realizations that are already in group $k$. \hyperlink{Aldous, 1985}{Aldous (1985)} show that $r_1, r_2, \ldots, r_n$ generated according to the Chinese restaurant process are i.i.d. draws from $F$\footnote{The realisations $r_1, r_2, \ldots, r_n$ generated according to ((ref)) are not independent given that the $n^{th}$ realisation is generated conditioned on the $n-1$ realizations before. However, these realisations are exchangeable, and therefore independent conditional on a distribution $F$.}, i.e.,
A model with a $\mathcal{DP}$ prior on the distribution of parameters is called a $\mathcal{DP}$ mixture model (c.f., \hyperlink{de Carvalho et al., 2013}{de Carvalho et al., 2013}; \hyperlink{Wiesenfarth et al. 2014}{Wiesenfarth et al., 2014}; \hyperlink{Li et al., 2018}{Li et al., 2018} and \hyperlink{Hejblum et al., 2019}{Hejblum et al., 2019}), and is capable of representing very general forms of heterogeneity in the distributions of the observations. The $\mathcal{DP}$ normal mixture model, whose mixture components are normal distributions, can be written as
where $\theta_i$ is the set of parameters of observation $\bm{y}_i$. In the multivariate normal\footnote{Note that $\bm{y}_i$ is a vector.} case, $\theta_i$ consists of the mean vector and covariance matrix, i.e., $\theta_i = \left( \bm{\mu}_i, \bm{\Sigma}_i \right)$.
The posterior probability of $\theta_i$ having the same value as one of the existing $\theta_{-i}$ is
where $\theta^*_{k}$ and $n_{k}$ respectively denote the unique value of group $k$ and the number of observations already in group $k$. $f_{\mathcal{N}} (\cdot)$ denotes the density function of multivariate normal distribution. The posterior probability of $\theta_i$ taking a new value $\theta^*_{new}$ from the base distribution is
where $f_{F_0} \left( \theta_{i} \right)$ is the probability density of the current value $\theta_{i}$ given $F_0$.
In this section we introduce the generic form of the semi-parametric gls estimator, where a $\mathcal{DP}$ prior is introduced on the distributional parameters of the errors. Consider a general linear regression
where $i$ indexes the observation, $\bm{y}_i$ is a $Q \times 1$ vector of dependent variables, $\bm{X_i}$ is a $Q \times K$ matrix of explanatory variables, and $\bm{\beta}$ is a $K \times 1$ vector of coefficients. $\bm{\varepsilon}_i \sim \mathcal{N} \left( \theta_i \right)$ is a $Q \times 1$ vector of errors, where $\theta_i = \left( \bm{\mu}_i, \bm{\Sigma}_i \right)$ denotes the distribution parameters.
To facilitate a more flexible distribution for the errors, where both the mean and covariance matrix of each $\bm{\varepsilon_i}$ is allowed to be different, we do not include a constant in $\bm{X}_i$. Our semi-parametric gls estimator introduces a $\mathcal{DP}$ prior on the distribution of $\theta_i$, resulting in the errors having a non-parametric mixture of normal distributions. The hierarchical prior for $\theta_i $ is
where $\alpha$ is the concentration parameter, and $F_0$ is the base distribution of the $\mathcal{DP}$, respectively.
Due to the discreteness of $F$ under the $\mathcal{DP}$ prior, the values of some $\theta_i$'s will be the same, thus putting them into the same group. The distribution parameters of $\bm{\varepsilon}_i$ can then be written as
where $c_i$ denotes the group id of $\bm{\varepsilon}_i$, and the superscript $*$ denotes the group-specific values of parameters. If $c_i = c_j$, $i, j \in \{ 1, 2, \ldots, N \}$, $\bm{\varepsilon}_i$ and $\bm{\varepsilon}_j$ share the same group id and parameters, i.e., $\bm{\mu}_{c_i}^* = \bm{\mu}_{c_j}^*$ and $\bm{\Sigma}_{c_i}^* = \bm{\Sigma}_{c_j}^*$. Such “grouping” characteristic can help to reveal the structure of the unobserved heterogeneity in the data.
A straightforward choice of the base distribution is the normal-inverse Wishart distribution, i.e., $\theta \sim \mathcal{NIW} (\nu_0, \bm{W}_0, \bm{\lambda}_0, \kappa_0)$, or equivalently,
where $\nu_0$, $\bm{W}_0$, $\bm{\lambda}_0$ and $\kappa_0$ are the hyper-parameters of the normal-inverse Wishart distribution. As the normal-inverse Wishart distribution is the conjugate prior for the parameters of a multivariate normal distribution, such a base distribution simplifies the evaluation of the posterior probability of $\theta_i$ taking a new value in ((ref)). Having specified the $\mathcal{DP}$ prior, we now outline the process of drawing from the posterior distributions of the parameters.
In order to take draws for the parameters in the model, the following Gibbs sampler is employed. We need to draw the distributional parameters $\Theta = \{ \theta_{c_i}^* \}_{i=1}^{N} = \{ \bm{\mu}_{c_i}^*, \bm{\Sigma}_{c_i}^* \}_{i=1}^{N}$, the regression parameters $\bm{\beta}$ and the concentration parameter $\alpha$ of the $\mathcal{DP}$ prior, i.e.,
We begin with drawing the distributional parameters of the errors. For this we view the residuals $\bm{e}_i$ ($i = 1, \ldots, N$) as the “observed” errors. According to ((ref)) in the Chinese restaurant process, the probability of distributional parameters $\theta_i$ being assigned to group $c_j$ is
Similarly by ((ref)), the probability of $\theta_i$ being assigned to a new group, i.e., taking a new value drawn from $F_0$ is
Normalizing $\tilde{p} = \{ \tilde{p}_1, \ldots, \tilde{p}_N \}$ gives us a multinomial probability vector
A draw can then be taken within $1, \ldots, N$ from the corresponding multinomial distribution to update the group membership. If the multinomial draw is $j \neq i$, $\theta_i$ is assigned to the group $c_j$, while if the draw is $i$, $\theta_i$ is assigned to a new group.
It should be mentioned that due to the choice of normal-inverse Wishart distribution for $F_0$, the integral in ((ref)) is
where $Q$ is the dimension of the error term vector $\bm{\varepsilon}_i$, $\Gamma_Q(\cdot)$ is the multivariate gamma function, and
As for the value of the hyper-parameters, we follow \hyperlink{Conley et al., 2008}{Conley et al. (2008)} and set $\nu_0 = Q + 0.004$, $\bm{W}_0 = 0.17 \bm{I}$, $\bm{\lambda}_0 = \bm{0}$ and $\kappa_0 = 0.016$.
Once the group memberships of all the distributional parameters $\theta_i$ have been updated, the unique values of parameters for each group are redrawn (as suggested by \hyperlink{Escobar and West, 1998}{Escobar and West, 1998}). Let $N_k$ and $\bar{\bm{e}}_k$ be the number of residuals assigned to group $k$ and their sample mean, respectively. The unique values of this group $\theta_{k}^* = (\bm{\mu}_k^*, \bm{\Sigma}_k^*)$ is redrawn from a normal-inverse Wishart distribution denoted by $\mathcal{NIW} (\nu_k, \bm{W}_k, \bm{\lambda}_k, \kappa_k)$, whose parameters are
where $\bm{S}_k = \sum_{c_i = k}{\left( \bm{e}_i - \bar{\bm{e}}_k \right) \left( \bm{e}_i - \bar{\bm{e}}_k \right)'}$.
After draws have been taken for the distributional parameters, we move on to the second part of the mcmc algorithm in (ref), which draws from the posterior of the regression parameters $\bm{\beta}$. The likelihood of $\bm{\beta}$ in ((ref)) is
We specify a normal prior for $\bm{\beta}$, i.e.,
where $\bm{b}_0$ and $\bm{V}_0$ denote the prior mean and covariance matrix of $\bm{\beta}$, respectively. The posterior of $\bm{\beta}$ is
where
and
For the hyper-parameters we specify $\bm{b}_0 = \bm{0}$ and $\bm{V}_0 = 1000 \bm{I}_K$, in order to prevent a prior that is overly informative.
Finally, to make draws of the concentration parameter $\alpha$, we adopt the prior introduced by \hyperlink{Conley et al., 2008}{Conley et al. (2008)}, which is
where $\alpha_{min}$ and $\alpha_{max}$ are the pre-set lower and upper bounds of $\alpha$. Larger $\alpha$ leads to more groups being generated on average, i.e., the $\mathcal{DP}$ being less discrete. \hyperlink{Antoniak, 1974}{Antoniak (1974)} gives the distribution of $K$, the number of groups, conditioned on $\alpha$. The corresponding posterior of $\alpha$ is
We set $\alpha_{min}$ to 0.1083 such that the mode of $K$ is 1, and set $\alpha_{max}$ to 1.834 so that the mode of $K$ is 5% of the sample size\footnote{To test the sensitivity to the prior of $\alpha$, the hyper-parameter $\alpha_{max}$ has been adjusted so that the mode of $K$ is 10% and 50% of the sample size, respectively. In our experiments the results are not sensitive to these changes in $\alpha_{max}$.}. Following the suggestion of \hyperlink{Conley et al., 2008}{Conley et al. (2008)}, we set $\tau$ to 0.8.
In this section we introduce how the $\mathcal{DP}$ prior is incorporated with the sur for equation systems. Consider a system of $M$ equations
where $\bm{y}_m = \left[y_{mi}\right]_{m=1}^{M}$, $\bm{\varepsilon}_m = \left[\varepsilon_{mi}\right]_{m=1}^{M}$ are $N \times 1$ vectors of dependent variables and errors, respectively. $\bm{X}_m$ is an $N \times K_m$ matrix of explanatory variables. Following our assumptions in (ref), there are no constants in the equations. $\bm{\beta}_m$ is a $K_m \times 1$ vector of parameters.
In the presence of correlation between errors across the equations there exists an efficiency gain by utilising a system estimator. The sur (\hyperlink{Zellner, 1962}{Zellner, 1962}) was introduced for this task. The equations are stacked in the following way
or simply
In comparison with the Bayesian ols estimator that usually assumes $\varepsilon_{mi} \overset{iid}{\sim} \mathcal{N}(\mu_m, \sigma_m^2)$ for all $m$, the errors $\bm{\varepsilon}_i$ are now identically multivariate normally distributed, i.e., $\bm{\varepsilon}_i \overset{iid}{\sim} \mathcal{N}(\bm{\mu}, \bm{\Sigma})$. The covariance matrix of $\bm{\varepsilon}$ is then
where ”$\otimes$" stands for the Kronecker product.
The gls type estimators (sur in this case) utilize the information in the covariance matrix to transform the data, so that the transformed errors are homoskedastic with no serial correlation. Such transformation is reflected in the likelihood as in ((ref)).\footnote{However, the covariance matrix $\bm{\Omega}$ has a specific form as in the sur in ((ref)), instead of the general, positive definite symmetric form of a covariance matrix.} The prior and posterior of the parameters can be defined similarly to equations ((ref)) to ((ref)).
Although sur accounts for the cross-equation correlation of errors, as \hyperlink{Wooldridge, 2003}{Wooldridge (2003)} has noted, the errors are assumed to be identically distributed. Moreover, unlike the frequentist gls estimator, this distribution is usually assumed to be normal in Bayesian methods. In this section we propose a new dp-sur estimator that makes no a priori assumptions on the family of distribution of the errors. If we allow each observation $i$ to have its own distributional parameters, flexibility of the error distribution will lead to identification problems given cross sectional data. A compromise is to assign the observations into groups with the $\mathcal{DP}$ prior as in ((ref)), in which case the distribution parameters of $\bm{\varepsilon}_i$ of the sur in ((ref)) is
The main difference from the parametric Bayesian sur is that the covariance matrix of each observation is now given in ((ref)), which allows each group of observations to have its own unique values for the parameters. Then draws can be taken from the posteriors of parameters according to the Gibbs sampler described in Section (ref).
In this section we conduct simulation experiments designed to evaluate and compare the performances of three estimators. The first one is our semi-parametric dp-sur in Section (ref). The second is a Bayesian sur, where the errors have a parametric mixture of normal distributions with a Dirichlet prior (abbreviated as dir-sur hereafter) on the distributional parameters of the errors. The third is a parametric Bayesian sur where the errors have a normal distribution (abbreviated as nor-sur hereafter).
With not loss of generality, all simulation experiments are based on a two-equation system. Two explanatory variables are included in each equation, which are drawn from normal distributions with the following parameters
and
We set the values of the parameters to $\bm{\beta}_1 = (1, -2)^{\prime}$ and $\bm{\beta}_2 = (-1, 2)^{\prime}$.
We generate errors from two types of distributions to demonstrate our dp-sur. The first is the multivariate log-normal distribution, i.e., $\ln \bm{\varepsilon_i} \sim \mathcal{N} \left( \bm{0}, \bm{\Sigma} \right)$. The log-normal distribution is fat tailed with a positive mean, thus suitable for examining the performance of our semi-parametric dp-sur. For the covariance matrix of the bivariate normal distribution of $\ln \bm{\varepsilon}_i$, we specify its form as
Without loss of generality, we let the variances of $\varepsilon_{1i}$ and $\varepsilon_{2i}$ be identical, and fix the correlation between them at 0.5. In order to explore the performances of the estimators when extreme values are generated with different probabilities, we set $\sigma^2$ equal to three values: 1, 1.5 and 2.
The second type of error distribution is designed to demonstrate the performance of our dp-sur when the errors have multi-modal distributions. In order to adjust the heaviness of the tails of the distributions as well, we employ mixtures of multivariate Student t distributions, which are scale mixtures of multivariate normal distributions (\hyperlink{Andrews and Mallows, 1974}{Andrews and Mallows, 1974}). To avoid unnecessary complexity, we utilise a mixture of two multivariate t distributions with the same degrees of freedom ($df$) and scale matrix, but with different means. More specifically, the scale matrix shared by both multivariate t distributions is
Four $df$ are specified to adjust the heaviness of the tails of the mixture distribution, namely 2\footnote{We avoid the df of 1 because the multivariate t distribution has no well defined expectation under this circumstance.}, 4, 6 and $\infty$, with the tails becoming less heavy. Note that when the $df$ is $\infty$, the mixture distribution reduces to a mixture of multivariate normal distributions. The mean vectors of the two multivariate t are $[-1, -1]'$ and $[4, 4]'$. They are chosen to ensure that the mean vectors of the two mixture components are located far enough from one another, so that the mixture distributions are bi-modal for all four $df$. The two multivariate t distributions are allocated mixture weights of 0.4 and 0.6, respectively, resulting in the mixture distributions being asymmetric.
Three sample sizes are chosen for the simulation experiments: 100, 200, and 300. For each sample size, 100 samples are generated. In the following section we will report the average posterior means, standard deviations and mean squared errors over the 100 samples.
In this section we present the results of the simulations for the three estimators dp-sur, dir-sur and nor-sur.
In Table (ref) we report the posterior means, standard deviations and the mean squared errors of the parameters averaged over the 100 samples, estimated with the three estimators when the errors are multivariate log-normal. The columns are divided into three blocks for the three sample sizes, each containing the results of the three estimators. The rows are also divided into three blocks for the three values of the variance $\sigma^2$. Each block contains the posterior mean, standard deviation and mean squared errors of the parameters. It should be noted that the dir-sur uses a fixed number of mixture components, which in our simulation studies is set to the posterior mode of the number of clusters given by the dp-sur averaged over the 100 samples.
One may see from Table (ref) that all three estimators give posterior means that are very close to the true values of the parameters. This indicates that all three estimators provide good point estimators for the parameters.
Regarding the posterior standard deviations, it is not surprising that all three estimators provide smaller posterior deviations when the sample size increases for each of the three $\sigma^2$. Table (ref) also shows that the dp-sur posterior standard deviations are uniformly smaller than the dir-sur ones, which are in turn smaller than the nor-sur posterior standard deviations. It is worth noticing that when $\sigma^2$ gets larger, the ratio of the dp-sur posterior standard deviations to the dir-sur become smaller for all sample sizes. This indicates that through a non-parametric mixture of normal distributions, our dp-sur achieves increasingly less dispersed posteriors than the \textsc{dir-sur} that employs a parametric mixture when the “true” distribution of the errors is more skewed and heavy tailed.
The same phenomena are also observed comparing the dp-sur posterior standard deviations to those of the nor-sur. In addition, the posterior standard deviations of the dp-sur are rather similar across the three $\sigma^2$ for the same sample size, which demonstrates its capability to fit more skewed and heavy tailed error distributions without having to make the posterior more dispersed. In contrast, both parametric estimators, i.e., the dir-sur and the nor-sur, record considerable larger posterior standard deviations when $\sigma^2$ increases.
The mean squared errors (mse) are one of the most important measures for the performance of estimators. Results in Table (ref) show that the mse of all three estimators decrease as sample size increases for all values of $\sigma^2$. We also see from Table (ref) that the dp-sur again dominates the two parametric estimators dir-sur and nor-sur, giving smaller mse under all circumstances. Similar to the behaviour of the posterior standard deviations\footnote{This is not surprising, for the mean squared error is the sum of the squared bias of the estimator and its variance. As all three estimators in our simulations achieved posterior means close to the true values of the parameters, the differences in their \textsc{mse} are mostly driven by their posterior standard deviations.}, for all sample sizes the ratios of the \textsc{dp-sur} \textsc{mse} to the \textsc{dir-sur} and \textsc{nor-sur} ones become smaller as $\sigma^2$ increases. Moreover, the \textsc{mse} of the \textsc{dp-sur} remain similar when $\sigma^2$ becomes larger, while those of the parametric \textsc{dir-sur} and \textsc{nor-sur} both increase considerably. This again indicates that our semi-parametric \textsc{dp-sur} has superior performance when the “true” distribution of the errors is fat tailed.
The results with the mixed multivariate t errors are presented in Table (ref). The four horizontal blocks represent different $df$, i.e., 2, 4, 6 and $\infty$. They show that the posterior means given by all three estimators are close to the truths, indicating that they all give good point estimators for the parameters.
The posterior standard deviations given by our semi-parametric dp-sur are smaller than those of the dir-sur and nor-sur under all circumstances. This demonstrates that our dp-sur produces posteriors that are considerably more concentrated than the two parametric estimators. It is worth noticing that for each sample size, the ratios of the dp-sur posterior standard deviations to those of the dir-sur and \textsc{nor-sur} are the smallest when $df = 2$, and increases with larger $df$. This is due to the tails of the multivariate t distributions mixed being heavier when $df$ is smaller. The same situation is observed with respect to the advantage of \textsc{dir-sur} over \textsc{nor-sur}. In addition, the \textsc{dir-sur} and \textsc{nor-sur} posterior standard deviations are almost the same when $df = \infty$, while \textsc{dp-sur} posterior standard deviations are still smaller, which indicates that the parametric mixture model experiences difficulties in identifying the heterogeneity in the error distribution.
Similar to the posterior standard deviations, the mse of the dp-sur are uniformly smaller than the dir-sur and nor-sur ones. This indicates that our semi-parametric dp-sur outperforms the parametric estimators when the error distribution is asymmetric and bi-modal. In addition, the advantages regarding mse of the \textsc{dp-sur} over the \textsc{dir-sur}/\textsc{nor-sur} are the largest when the $df$ is 2, and becomes smaller when the $df$ increases. That is, the semi-parametric \textsc{dp-sur} proves more efficient when the tails of the mixed multivariate t distribution are heavier. Moreover, one can observe that the \textsc{mse} of the \textsc{dir-sur} and the \textsc{nor-sur} are similar when $df = \infty$, while the \textsc{dp-sur} still has smaller \textsc{mse}. It demonstrates that the performance advantage of the \textsc{dir-sur} using a parametric mixture of normal distributions over the \textsc{nor-sur} diminishes faster than that of the \textsc{dp-sur} as the tails of the error distribution become less heavy.
In this section we apply our dp-sur estimator to the demand for factors of production with a generalized Leontief cost function (\hyperlink{Diewert, 1971}{Diewert, 1971}). The model is an equation system with as many equations as there are factors. We do not impose symmetry or homogeneity restrictions to make the model more general.
The dataset we use are from \hyperlink{Malikov et al. 2016}{Malikov et al. (2016)}, which contains 799 observations on 285 large U.S. banks in 2002, 2004 and 2006. The data includes quantities and prices of the inputs including labour, physical assets and borrowed funding, and the quantity of output, which is the loans made by a bank.
The equation system for the demands for factors is
where $L$, $A$ and $F$ denote the quantity of labour, physical assets and borrowed funds, respectively; $T$ denotes the trend variable; $Y$ denotes output, and $P_k$ is the price of factor $k$, with $k \in \{L, A, F\}$. For the errors we assume that conditional on the explanatory variables, $\left( \varepsilon_{Li}, \varepsilon_{Ai}, \varepsilon_{Fi} \right) \sim \mathcal{N} \left(\bm{\mu}_i, \bm{\Sigma}_i \right)$ where $\bm{\mu}_i$ and $\bm{\Sigma_i}$ are respectively the mean vector and covariance matrix of individual $i$. We allow the errors to be correlated across the three equations in the system for any particular individual, i.e., $\text{cov} \left( \varepsilon_{ki}, \varepsilon_{si} \right) \neq 0$, with $k, s \in \{L, A, F\}$ indexing equations.
With the generalized Leontief cost function, the cross price elasticities of the factors are given by
The own price elasticities are
As the main interest in a demand system for factors of production is the price elasticities, we report results regarding the price elasticities of the three factors. There are nine elasticities within the three-input system, which are indexed by $\{L, A, F \}$ as stated before. In Table (ref) we present the posterior means, standard deviations and p-values of the elasticities by the three estimators: dp-sur, dir-sur and nor-sur. In addition, we also demonstrate the performances of the three estimators with several figures.
From Table (ref) one can see that the posterior means given by the three estimators are of the same signs and relatively close in magnitude for all nine elasticities. Among them, the posterior means of all three own-price elasticities are negative, with labour being the most elastic, followed by physical assets, and funding is the lease elastic input, which is expected with the banking industry. The cross-price elasticities demonstrate that labour and assets are substitutes, while assets and funding are complements. The relationship between labour and funding is not very clear: the demand for labour will decrease when the price for funding rises as indicated by $E_{LF}$, but $E_{FL}$ shows that the demand for funding will increase when labour becomes more expensive. However, it should be noted that both $E_{LF}$ and $E_{FL}$ are rather small in magnitude, showing the both are inelastic to changes in the price of the other.
The posterior standard deviations of our semi-parametric dp-sur are always smaller than the dir-sur or nor-sur ones. This difference is particularly significant with the cross-price elasticities regarding funding. Such a difference contributes to the fact that only the dp-sur posterior p-values are below the usually chosen significance level of 0.05 for $E_{FL}$ and $E_{FA}$.
To better illustrate the performances of the three estimators, we present the 95% highest posterior density intervals (hpdi's) of the elasticities in the demand system of U.S. banks for factors in Figure (ref). We observe that the hpdi's of dp-sur are the narrowest intervals for all nine elasticities. For instance, in the case of the elasticity of funding w.r.t. the asset price (e_fa), only our dp-sur shows significance at the 5% level.
We also present the histograms of $E_{FA}$ to further demonstrate the performance of the three estimators in Figure (ref). One may see that only the dp-sur gives such an hpdi that excludes 0. From the histogram one may also see that the posterior given by our dp-sur is left skewed with the non-parametric mixture. In fact, the posterior mode of the number of clusters is 3, showing that heterogeneity is detected in the error distributions.
In Figure (ref) we present the posterior means of predictive densities of the residuals obtained by the dp-sur and nor-sur. One may see that with the semi-parametric dp-sur, the predictive density of the residuals is bi-modal, which is not captured by the nor-sur.
The gls has also seen numerous applications with panel data models, and in particular the random effects model (rem). In a panel with $N$ units and $T$ time periods, the error of each unit\footnote{We use the term “unit” to denote the cross section here. In practice it can be households, firms, countries or persons.} is a $T \times 1$\footnote{That is, the $Q$ in the generic semi-parametric gls in Section (ref) is $T$ in this context. We use $T$ here following panel data protocols.} vector. We will relax the assumption of parametric Bayesian gls for the rem (\hyperlink{Koop, 2003}{Koop, 2003}) that the error vectors for all units have the same distribution. In this section we propose a semi-parametric Bayesian approach by introducing $\mathcal{DP}$ priors on the distributional parameters of the random effects and idiosyncratic errors. We follow the same approach as in the dp-sur method in terms of applying the $\mathcal{DP}$ prior on the distributional parameters.
Consider the following panel data model
where $i$ and $t$ index the cross section and time series dimensions of the data, respectively. $y_{it}$ is the dependent variable, $\bm{x}_{it}$ denotes the explanatory variables, and $\bm{\beta}$ is a conformable vector of parameters. $u_i$ is the time-invariant unobservable of unit $i$, and $\eta_{it}$ the idiosyncratic error term. $\varepsilon_{it} = u_i + \eta_{it}$ is the composite error.
In Bayesian methods the difference between the fixed and random effects lies in the choice of prior for the individual effects $u_i$: the fixed effects model assumes a non-hierarchical prior for $u_i$; and the rem includes a hierarchical prior. The prior for $u_i$ in the rem can be written as
where $\mu_u$ and $\sigma_u^2$ are the mean and variance\footnote{Note that the distributional parameters of $u_i$ and $\eta_{it}$ are often assumed to be random, and have their own priors. However, for the moment we leave them fixed for the sake of simplicity.} of $u_i$, respectively. Assuming $\eta_{it} \overset{iid}{\sim} \mathcal{N} (\mu_{\eta}, \sigma_{\eta}^2)$, the likelihood of $\bm{\beta}$ marginalized over $u_i$ in the Bayesian rem may be written as
where $\bm{X}_i = \left[x_{1it}, \ldots, x_{Kit}\right]_{t=1}^{T}$ is a $T \times K$ matrix of explanatory variables, and $\bm{y}_i = [y_{it}]_{t=1}^{T}$ is a $T \times 1$ vector of dependent variables. $\mu = \mu_u + \mu_{\eta}$ is the mean of the composite error $\varepsilon_{it}$, and $\iota_T$ is a $T \times 1$ vector of ones. $\bm{\Sigma}$ is the covariance matrix of the $T \times 1$ composite error vector $\bm{\varepsilon}_i = [\varepsilon_{i1}, \varepsilon_{i2}, \ldots, \varepsilon_{iT} ]'$. Assuming the usual strict exogeneity in rem, $\bm{\Sigma}$ is
In this paper our primary focus is the provision of more efficient inference by exploiting information in the heterogeneous distributions of unobservables. Thus, we maintain the approach of gls, and focus on heterogeneity in the distributional parameters of the unobservables instead of introducing heterogeneity to the distribution of model parameters themselves\footnote{ \hyperlink{Kleinman and Ibrahim, 1998}{Kleinman and Ibrahim (1998)} and \hyperlink{Kyung et al., 2010}{Kyung et al. (2010)} studied the heterogeneity in the model parameters across the units ($\bm{\beta_i}$). In contrast to our case, the $\mathcal{DP}$ prior was put on the parameters themselves in these papers.}. In this sense, our method is in the same spirit as the literature pioneered by \hyperlink{Conley et al., 2008}{Conley et al. (2008)}. We relax the identical distribution assumptions for both $\eta_{it}$ and $u_i$ by introducing two independent $\mathcal{DP}$ priors on their distributional parameters\footnote{\hyperlink{Hirano, 2002}{Hirano (2002)} presents a semi-parametric autoregressive panel data model with individual effects. A $\mathcal{DP}$ prior is introduced on distributional parameters of idiosyncratic errors. Our estimator, although not considering dynamic panels, allows both the idiosyncratic errors and the individual effects to be heterogeneous in their distributions by introducing two independent $\mathcal{DP}$ priors for them, which leads to more flexibility.}.
The $\mathcal{DP}$ prior for $\theta_{\eta, it} = ( \mu_{\eta, it}, \sigma^2_{\eta, it})$, the parameters of the idiosyncratic error $\eta_{it}$, is
where $\alpha_{\eta}$ and $G_0$ denote the concentration parameter and base distribution of the $\mathcal{DP}$ prior, respectively. As the prior is introduced across all idiosyncratic errors, the grouping of distributional parameters $\theta_{\eta, it}$ are not restricted to either within the cross section or within the time series dimensions. This means $\eta_{it}$ and $\eta_{is}$ can be allocated to different groups $c_{it}$ and $c_{is}$ such that they have different distributions with parameters $\theta_{\eta, c_{it}}^* = (\mu_{\eta, c_{it}}^*, \sigma_{\eta, c_{it}}^{2*})$ and $\theta_{\eta, c_{is}}^* = (\mu_{\eta, c_{is}}^*, \sigma_{\eta, c_{is}}^{2*})$, respectively. Similarly, $\eta_{it}$ and $\eta_{jt}$ can be in the same group, i.e., $c_{it} = c_{jt}$, making them identically distributed with parameters $\theta_{\eta, c_{it}}^* = \theta_{\eta, c_{jt}}^*$.
We write the $\mathcal{DP}$ prior for $\theta_{u, i} = (\mu_{u, i}, \sigma^2_{u, i})$, the parameters of the individual effects $u_i$'s, as
where $\alpha_u$ is the concentration parameter, and $F_0$ is the base distribution of the $\mathcal{DP}$ prior. Given that we maintain the assumption of strictly exogeneity of rem, it is reasonable to introduce a $\mathcal{DP}$ prior independent of $G$\footnote{For two mixtures of normal distributions, it is possible to introduce, e.g., a Hierarchical Dirichlet Process prior (\hyperlink{Teh et al., 2005}{Teh et al., 2005}) instead of two independent $\mathcal{DP}$ when dependence between the two priors is necessary. However, as we are considering the individual effects $u_i$ and idiosyncratic errors $\eta_{it}$ in the rem framework, it is plausible to introduce two independent $\mathcal{DP}$ priors in our case.} for the distributional parameters $\theta_{u, i}$.
The $\mathcal{DP}$ prior on the distributional parameters of individual effects $u_i$ generates groupings across the $N$ units. As a result, if $u_i$ and $u_j$ are in different groups $c_i$ and $c_j$, their parameters will take the values $\theta_{u, c_i}^* = (\mu_{u, c_i}^*, \sigma_{u, c_i}^{2*})$ and $\theta_{u, c_j}^* = (\mu_{u, c_j}^*, \sigma_{u, c_j}^{2*})$, respectively. This relaxes the rem assumption that the individual effects are identically distributed, as $u_i$ and $u_j$ are allowed to have different distributional parameters.
The mean vector of the composite error vector $\bm{\varepsilon}_i = u_i \iota + \bm{\eta}_i$ is then given by
where $\bm{\mu}_{\bm{\eta}_i}$ is the mean vector of $\bm{\eta}_i$. The covariance matrix of $\bm{\varepsilon}_i$ is
where $\bm{\Sigma}_{\bm{\eta}_i} = \mathbf{diag} \left(\sigma_{\eta, c_{i1}}^{2*}, \ldots, \sigma_{\eta, c_{iT}}^{2*} \right)$ is the diagonal covariance matrix of $\bm{\eta}_i$, where the diagonal elements are the variances of the idiosyncratic errors of unit $i$ over all time periods.
The mcmc algorithm utilised to draw from the posterior distributions is similar to that described in Section (ref) with a few minor changes due to the particular form of the dp-rem. The most significant change is that now the individual effects $\bm{U} = \{u_i\}_{i=1}^{N}$ must also be drawn. In addition, there are two sets of distributional parameters, i.e., $\Theta_u$ and $\Theta_\eta$ for $u_i$ and $\eta_{it}$, respectively. Similarly, there are concentration parameters for the two independent $\mathcal{DP}$ priors, namely $\alpha_u$ and $\alpha_{\eta}$. Accordingly, the Gibbs sampler consists of
Among the parameters, the drawing of the concentration parameters $\alpha_u$ and $\alpha_{\eta}$ is exactly the same as described in Section (ref), because the two $\mathcal{DP}$ priors in our dp-rem are independent. The drawing of $\Theta_u$ and $\Theta_\eta$ are also similar to the procedure described in ((ref)) to ((ref)). It should be noted that the normal-inverse Wishart distribution in ((ref)) can still serve as the base distributions $G_0$ and $F_0$. However, as both $u_i$ and $\eta_{it}$ are scalars, they become the normal-inverse gamma distribution with both $\bm{\lambda}_0$ and $\bm{W}_0$ being scalars consequently.
The drawing of $\bm{U}$ is based on the fact that for a given $i$, the hierarchical prior for $u_i$ is
which is the conjugate prior for the normal likelihood, leading to the posterior of $u_i$ being normal as well. The posterior variance of $u_i$ is
The posterior mean of $u_i$ is
As for the regression parameter $\bm{\beta}$, its likelihood marginalized over ${u_i}$ is given by
which differs from ((ref)) only in that $\bm{\mu}_i$, the mean of the composite error $\bm{\varepsilon}_i$ is included. Compared with the marginal likelihood of the parametric Bayesian rem in ((ref)), the covariance matrix of the composite error vector $\bm{\varepsilon_i}$ is allowed to be different for each unit $i$ in the panel.
Given a conjugate normal prior for $\bm{\beta}$, i.e.,
where $\bm{b}_0$ and $\bm{V}_0$ respectively denote the prior mean and covariance matrix of $\bm{\beta}$, the posterior covariance matrix of $\bm{\beta}$, marginalized over $u_i$, is given by
The posterior mean vector is
A modified version of our dp-rem can be introduced in the spirit of the correlated rem introduced by \hyperlink{Mundlak, 1978}{Mundlak (1978)} and further discussed by \hyperlink{Chamberlain, 1982}{Chamberlain (1982)}, \hyperlink{Wooldridge, 2005}{Wooldridge (2005)}, \hyperlink{Murtazashvili and Wooldridge, 2008}{Murtazashvili and Wooldridge (2008)} and \hyperlink{Wooldridge, 2019}{Wooldridge (2019)}. This model offers a middle ground between the fixed and random effects models by allowing the individual effects to be correlated with $\bm{X}_i$ in a specific manner. This is achieved by specifying the individual effects as a linear function of the within unit means of the explanatory variables, i.e.,
where $\bar{x}_{ki} = 1/T \sum_{t=1}^{T}{x_{kit}}$. Then a $\mathcal{DP}$ prior as in ((ref)) can be applied to $u_i$, and the analyses of the parameters are similar.
We carry out a series of simulation experiments to demonstrate the performance of our dp-rem. Similar to the simulations for the dp-sur, our dp-rem is compared with two other estimators. The first one is a Bayesian rem with a Dirichlet prior on the distributional parameters of $u_i$ and $\eta_{it}$ (dir-rem hereafter), where $u_i$ and $\eta_{it}$ have parametric mixtures of normal distributions. The second one assumes that $u_i$ and $\eta_{it}$ are normal distributed (nor-rem hereafter).
The simulations are carried out using the following specification
where the explanatory variables are generated from
with
The individual effects $u_i$ and the idiosyncratic errors $\eta_{it}$ are independently generated. As with the simulation experiments for the dp-sur, we specify two types of distributions: a log-normal distribution and a bi-modal distributions. For the log-normal distribution, the variances of $\ln u_i$ and $\ln \eta_{it}$ are chosen to be equal and denoted by $\sigma^2$, which are set to three values: 1, 1.5 and 2, in order to explore the performance of our dp-rem under different circumstances.
To create a bi-modal distribution, we mix two non-central univariate t distributions with non-centrality parameters -1 and 4. Like in the experiment design in Section (ref), the $df$ are set to 2, 4, 6 and $\infty$ to adjust the heaviness of the tails. The weights assigned to the two non-central t distributions are 0.4 and 0.6, so that the mixture distribution is asymmetric.
To explore the performance of our dp-rem estimator with different sample sizes, we fix the number of periods in the panel at 3 and set the number of cross sections to 100, 200 and 300. 100 samples are generated for each sample size, and the results reported in this section are the average over all the samples.
Table (ref) presents the results of the simulations with log-normal distributed individual effects $u_i$ and idiosyncratic errors $\eta_{it}$. For the posterior means we observe that the three estimators all perform well. When $\sigma^2 = 2$, we observe the greatest differences for the dir-rem and nor-rem posterior means from the true values. It shows that the two parametric methods are slightly inferior as point estimators when the distributions of $u_i$ and $\eta_{it}$ in ((ref)) are more skewed and heavy tailed. This can be caused by the distribution of the composite errors being a convolution of two log-normal distributions, which are fat tailed. The two parametric estimators dir-rem and nor-rem struggle to give good point estimates for the parameters in this case, while our semi-parametric dp-rem uniformly performs well.
The posterior standard deviations of all three estimators decrease as the sample size increases for all three values of $\sigma^2$. Among them, the dp-rem posterior standard deviations are always smaller than dir-rem and nor-rem ones. This indicates that the semi-parametric dp-rem provides posterior distributions that are less dispersed than the two parametric methods, as a result of the non-parametric mixture of normal distributions. It can also be seen that the posterior standard deviations of our dp-rem remains similar when $\sigma^2$ increases for any given sample size. In comparison, the posterior standard deviations of the two parametric estimators increase considerably when $\sigma^2$ gets larger.\footnote{This result was also observed in the simulations for equation systems as in Table (ref).}
The mse for the dp-rem are always smaller than their dir-rem and nor-rem counterparts. The superior performance of the dp-rem demonstrates the advantage of the non-parametric mixture with the $\mathcal{DP}$ prior over the parametric mixture using the Dirichlet prior, and the homogeneous normal errors. As with the posterior standard deviations, the mse of the \textsc{dp-rem} remain similar when $\sigma^2$ increases, while the \textsc{dir-rem} and \textsc{nor-rem} ones increase considerably. This further demonstrates the superiority of our semi-parametric \textsc{dp-rem} when the individual effects and idiosyncratic errors have fat tailed distributions.
Table (ref) reports the results with mixed t distributed individual effects and idiosyncratic errors. It can be seen that for all three estimators their posterior means are close to the truth in all settings. That is, all three estimators perform well as point estimators.
The posterior standard deviations of our dp-rem are the smallest in all scenarios, indicating that its posteriors are more concentrated than the dir-rem and nor-rem ones. We also observe that for a given sample size the advantages of the dp-rem over the parametric dir-rem and nor-rem with respect to posterior standard deviations are the largest when $df = 2$, and become smaller as $df$ increases. This is the result of the tails of the t distributions mixed in the error distributions becoming less heavy. In addition, the \textsc{dir-rem} gives almost the same posterior standard deviations as the \textsc{nor-rem} when $df$ are 4, 6 and $\infty$. This demonstrates that the advantages of the \textsc{dp-rem} over the \textsc{dir-rem} and \textsc{nor-rem} fall at a slower rate than that of the \textsc{dir-rem} over the \textsc{nor-rem} when $df$ increases.
Our dp-rem dominates the two parametric estimators regarding mse under all circumstances. Similar to the posterior standard deviations, the advantages of the dp-rem over the dir-rem and nor-rem in terms of mse fall as the $df$ increases, as the tails of the mixed t distributions become less heavy. It should be noted that the \textsc{mse} of \textsc{dir-rem} and \textsc{nor-rem} are almost identical when $df = 4$. That is, the advantage of the \textsc{dir-rem} over the \textsc{nor-rem} has diminished. In contrast, the advantage of the \textsc{dp-rem} over the two parametric estimators is still present. This demonstrates that the parametric \textsc{dir-rem} and \textsc{nor-rem} are not as good as the semi-parametric \textsc{dp-rem} at identifying the heterogeneity in error distributions when the tails become less heavy.
As a demonstration of the dp-rem with real data, in this section we present the results regarding a model for wages of U.S. workers. The data are from \hyperlink{Cornwell and Rupert, 1988}{Cornwell and Rupert (1988)} with 595 individuals over a period of 7 years, from 1976 to 1982. In fact, it allows us to demonstrate the correlated rem in (ref), as outlined below. The model is given by
where the dependent variable is the logarithms of wages, and the explanatory variables are experience in years ($Exp$), dummies for marriage status ($Mar$) and the individual being female ($Fem$), as well as the years of education ($Edu$). As there are strong reasons to suspect that the unobserved individual effect $v_i$ is correlated with the explanatory variables due to omitted variables such as ability and motivation, the correlated rem is selected. It specifies $v_i$ as
where the sample averages of experience ($AExp$) and marriage status ($AMar$) of individual $i$ are included, as they are the two time variant variables in ((ref)). The posterior means, standard deviations and p-values of the parameters with the three estimators dp-rem, dir-rem and nor-rem are presented in Table (ref). We also demonstrate their respective performances with a number of figures.
From Table (ref) we observe that the posterior means are similar for all parameters. Among the explanatory variables, experience and education are positively correlated with workers' wages, while being married and female negatively influences wages. The semi-parametric dp-rem gives the smallest posterior standard deviations among all three estimators. As the posterior modes of the numbers of clusters for $u_i$ and $\varepsilon_i$ are 2 and 6, respectively, heterogeneity is detected in both the individual effects' and the idiosyncratic errors' distributions. The posterior p-values of the dp-rem are also the smallest for all parameters. It should be noted that both parameters in (ref) have posterior p-values less than 5%, indicating that $v_i$ are correlated with the explanatory variables, providing support for the application of the correlated rem.
In Figure (ref) we present the 95% hpdi's of the three estimators for all regression parameters in (ref). One can observe that the semi-parametric dp-rem yields the narrowest 95% hpdi's for all parameters, followed by the dir-rem with a parametric mixture. The parametric nor-rem has the widest hpdi's for all the parameters, though for some parameters its \textsc{hpdi}'s are rather similar to those of \textsc{dir-rem}, e.g. with the parameter $\beta_4$ for education.
Figure (ref) presents the histograms of the posterior draws for the education parameter ($\beta_4$). Though all three intervals exclude zero, the dp-rem one is the shortest, as a result of the efficiency gain from exploring the heterogeneity in the error distributions.
Figure (ref) presents the predictive densities of the fitted individual effects and residuals obtained by the semi-parametric dp-rem and the parametric nor-rem, respectively. One can see that both the fitted individual effects and the residuals have multi-modal distributions with our dp-rem, which is not captured by the nor-rem. Thus, the normal distribution assumption could cause potential losses in efficiency.
In this paper we address the potential violation of the assumptions made by parametric Bayesian gls estimators that the errors are homogeneous regarding their distributions. Such assumptions are likely to be problematic in reality particularly when micro data are used, as the features of individuals or households are likely to lead to observations having different distributions. We present a semi-parametric Bayesian gls where the error distribution is a non-parametric mixture of normal distributions by introducing a Dirichlet process prior on the distributional parameters of the errors. The number of normal components is decided jointly by the data and the prior in such a mixture, which is able to cover a large variety of distributions. The errors are grouped by the $\mathcal{DP}$ prior, with those in the same group having the same distributional parameters and thus the same distribution. Two specific cases of the semi-parametric Bayesian gls are then introduced, which are the sur for equation systems and the rem for panel data.
Our dp-sur and dp-rem methods are demonstrated with a series of simulation experiments consisting of two scenarios, where the errors have a log-normal distribution and a mixture of t distributions that are bi-modal and asymmetric, respectively. When the errors have a log-normal distribution, which is fat tailed, our semi-parametric gls estimator gives smaller posterior standard deviations, as well as smaller mean squared errors in all settings. Such advantages over the parametric estimators are greater when the variances of the errors are larger, leading to heavier tails of the error distributions. When the errors are bi-modal as a mixture of t distributions, our semi-parametric gls estimators also out-perform the parametric estimators with respect to posterior standard errors and mean squared errors in all scenarios. Such advantages are larger when the degrees of freedom of the t mixture components are smaller, i.e., the tails of the mixture are heavier.
We apply our dp-sur method to the demands for production factors with the generalized Leontief cost function using a dataset of the U.S. banking industry. Heterogeneity is detected in the sample by our semi-parametric estimator. The dp-sur posterior standard deviations are smaller than the dir-sur ones using a parametric mixture of normal distributions, as well as the nor-sur ones for all the demand elasticities. In addition, the posterior p-values of the dp-sur are also smaller for all the parameters.
Our dp-rem is also applied to a study U.S. workers' wages, where there is a strong reason to suspect that the unobserved individual effects are correlated with the explanatory variables due to omitted variables describing individual features such as abilities. The correlated rem is then estimated. The dp-rem detects heterogeneity in the distributions of both the individual effects and the idiosyncratic errors with this sample. Our dp-rem obtains smaller posterior standard deviations and p-values than the parametric dir-rem and nor-rem for all parameters.
We are grateful to the useful comments received from Debopam Bhattacharya, Xiaohong Chen, Gernot Doppelhofer, Oliver Linton and Justin Tobias.