I have the following model:
$y_{it}=\alpha + x'_{it}\beta_{i} + \epsilon_{it}, \text{ } i=1,2,...,N, \text{ } t=1,2,...,T$ (1)
$\beta_{i}= z'_{i}\gamma+\eta_{i}$ (2)
with $\epsilon_{it} \sim N(0,\sigma_{\epsilon_{i}}^{2})$ and $\eta_{i} \sim N(0,\sigma_{\eta}^{2})$
Also i consider the following prior specifications:
$p(\beta,\gamma) \propto 1 \\ p(\sigma_{\epsilon_{i}}^{2}) \propto \sigma_{\epsilon_{i}}^{-2} \\ p(\sigma_{eta}^{2}) \propto \sigma_{\eta}^{-2}$
which has the following likelihood function:
$p(Y|\theta)=\prod_{i=1}^{N}\int_{\beta_{i}} \left( \prod_{t=1}^{T}\frac{1}{\sigma_{\epsilon_{i}}\sqrt{2\pi}}exp \left( -\frac{1}{2\sigma_{\epsilon_{i}}^{2}}(y_{it}-\alpha-x'_{it}\beta_{i})^{2}\right) \right) \times \frac{1}{\sigma_{\eta}\sqrt{2\pi}}exp \left( -\frac{1}{2{\sigma_{\eta}^{2}}}(\beta_{i}-z'_{i}\gamma)^{2} \right) d\beta_{i}$
I would like to estimate the parameters $\theta=(\alpha,\gamma,\sigma_{\epsilon_{i}}^{2},\sigma_{\eta}^{2})$. So $\theta$ is the parameter vector. Only $\beta_{i}$ is random. So I consider $\beta_{i}$ as the latent variables.
I would like to estimate these by using the following Gibbs sampling schema:
sample $\alpha$ given $\{\beta_{i}\}_{i=1}^{N},\gamma,\sigma_{\epsilon}^{2},\sigma_{\eta}^{2},Y$
sample $\gamma$ given $\{\beta_{i}\}_{i=1}^{N},\alpha,\sigma_{\epsilon}^{2},\sigma_{\eta}^{2},Y$
sample $\sigma_{\epsilon}^{2}$ given $\{\beta_{i}\}_{i=1}^{N},\gamma,\alpha,\sigma_{\eta}^{2},Y$
sample $\sigma_{\eta}^{2}$ given $\{\beta_{i}\}_{i=1}^{N},\gamma,\alpha,\sigma_{\epsilon}^{2},Y$
sample $\{\beta_{i}\}_{i=1}^{N}$ given $\gamma,\alpha,\sigma_{\epsilon}^{2},\sigma_{\eta}^{2},Y$
I know that:
To sample $\alpha$ I re-write (1) as $y_{it} - x'_{it}\beta_{i}=\alpha + \epsilon_{it}$, so the full conditional distribution of $\alpha$ is normal with mean $\hat{\alpha}=\frac{1}{T\times N} \sum_{i=1}^{N}\sum_{t=1}^{T} (y_{it} - x'_{it}\beta_{i})$ and variance $\frac{\sigma_{\epsilon}^{2}}{T \times N}$.
To sample $\sigma_{\epsilon}^{2}$ I consider the regression model in (1) , hence it can be sampled from an inverted Gamma-2 distribution with parameter $ \sum_{i=1}^{N}\sum_{t=1}^{T} (y_{it} - \alpha - x'_{it}\beta_{i})^{2}$ with $N \times T$ degrees of freedom.
To sample $\gamma$ I consider the regression model in (2), so $\gamma$ can be sampled from multivariate normal distribution with mean equal to the standard OLS estimators $\hat{\gamma}=\left(\sum_{i=1}^{N}z_{i}z'_{i}\right)^{-1}(z'_{i}\beta_{i})$ and covariance matrix equal to the standard OLS covariance matrix estimator $\sigma_{\eta}^{2}\left(\sum_{i=1}^{N}z_{i}z'_{i}\right)^{-1}$.
To sample $\sigma_{\eta}^{2}$ again we consider the linear regression in (2), hence $\sigma_{\eta}^{2}$ can be sampled from an inverted Gamma-2 distribution with parameter $\sum_{i=1}^{N}(\beta_{i}-z'_{i}\gamma)^{2}$ and $N$ degrees of freedom.
Lastly, we can sample $\{\beta_{i}\}_{i=1}^{N}$ from a normal distribution with mean \ $\left(\sum_{t=1}^{T}x_{it}^{2}+\left(\frac{\sigma_{\epsilon}}{\sigma_{\eta}}\right)^{2}\right)^{-1}\left(\sum_{t=1}^{T}x_{it}y_{it}+\left(\frac{\sigma_{\epsilon}}{\sigma_{\eta}}\right)^{2}z'_{i}\gamma\right)$ and covariance matrix $\sigma_{\epsilon}^{2}\left(\sum_{t=1}^{T}x_{it}^{2}+\left(\frac{\sigma_{\epsilon}}{\sigma_{\eta}}\right)^{2}\right)^{-1}$.
How can I do this in R ? If you could please give me an example for instance how to sample $\alpha$ I would be very grateful.
Thanks ! :)