Gibbs sampling in R

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 ! :)

In terms of writing this in R, here is an example I found: http://www.stat.purdue.edu/~zhanghao/MAS/handout/gibbsBayesian.pdf

You could then iterate through different levels of α and select the level yielding the best results.

• Thank you for your answer Hobbes. As far as I understand from the code, you (or the author of this document) initialise $\alpha$ and $\beta$ but not $\tau$. I guess this is because if you have calculated all the previous ones you can also compute $\tau$. But I would like to ask you, how could you set the type of flat prior for $\sigma_{\eta}^{2}$ (the $\sigma_{\eta}^{-2}$ ). ( For the $\beta$ and the $\gamma$ I know that you can just as prior a uniform distribution. Thanks ! Jun 21 '16 at 9:04
• This is not my work, just an interesting implementation I found. τ was initialized to 1 in this example. As for σ2η, this is, I believe, wrapped into the rgamma function and is defined by the shape, alpha + n/2. So I think it depends on your distribution and the α you decide to use. I hope I'm not misinterpreting it, this is an interesting problem but beyond my normal expertise. Jun 21 '16 at 15:59

Here is the code I wrote to answer my question. It might not be the most efficient one but it works. Sharing is caring :)

niter<-10000
N<-nrow(Y)

T<-10 # I will take into consideration until t=10 to estimate my parametres and then I will forecast the rest values t=11,12,... etc
result <- matrix(0,nrow = niter, ncol = 3) #Here I store 3 parameters: alpha, sigma_epsilon, sigma_eta. I do not store gamma here
result_betas<-matrix(0,nrow=niter,ncol=N) #the betas
result_gammas<-matrix(0,nrow=niter,ncol=nrow(Z)) #here I store the gammas

alpha<-0
gamma<-0

sigma_epsilon_squared<-1
sigma_eta_squared<-1
beta<-matrix(0, nrow=N, ncol=1)
beta2<-0

for (i in 1:niter)
{
mean_alpha<-0
for (obs in 1:N)
{
for (time in 1:T)
{
mean_alpha<-mean_alpha+Y[obs,time] - beta[obs]*X[obs,time]
}
}
alpha<-rnorm(1, mean = (1/(T*N))* mean_alpha, sd = sqrt((sigma_epsilon_squared)/(T*N)))

param_sigma_epsilon_squared<-0
for (obs in 1:N)
{
for (time in 1:T)
{
param_sigma_epsilon_squared<-param_sigma_epsilon_squared + (Y[obs,time] - alpha - beta[obs]*X[obs,time])^2
}
}
sigma_epsilon_squared <- rinvgamma(1, (N*T)/2, param_sigma_epsilon_squared/2 )

mean_gamma<-0 #this is a vector
mean_gamma1<-0 #this is a matrix
mean_gamma2<-0
sigma_gamma<-0
for(obs in 1:N)
{
mean_gamma1<-mean_gamma1 + Z[,obs]%*%t(Z[,obs]) #this is (k x k)
mean_gamma2<-mean_gamma2 + Z[,obs]*beta[obs]

#mean_gamma<-mean_gamma + 1/(Z[obs,]%*%t(Z[obs,])) * (Z[obs,]*beta[obs,1])
#sigma_gamma<-sigma_gamma+1/(Z[obs,]%*%t(Z[obs,]))
}
mean_gamma<-solve(mean_gamma1)%*%mean_gamma2
sigma_gamma<-solve(mean_gamma1)
gamma<-mvrnorm(n = 1, mu=mean_gamma, Sigma=sqrt(sigma_eta_squared)*sigma_gamma)

result_gammas[i,]<-gamma

mean_sigma_eta_squared<-0
for (obs in 1:N)
{
mean_sigma_eta_squared<- mean_sigma_eta_squared + (beta[obs] - t(Z[,obs])%*%gamma)^2
}
sigma_eta_squared<-rinvgamma(1, N/2, mean_sigma_eta_squared/2 )

result[i,]<-c(alpha,sigma_epsilon_squared,sigma_eta_squared)

mean_beta<-0
mean_beta1<-0 #this is scalar
mean_beta2<-0 #this is scalar
sigma_beta<-0
for (obs in 1:N)
{
for (time in 1:T)
{
mean_beta1<-mean_beta1 + X[obs,time]^2 + (sigma_epsilon_squared/sigma_eta_squared)
mean_beta2<-mean_beta2 + X[obs,time]*Y[obs,time] + (sigma_epsilon_squared/sigma_eta_squared)*(t(Z[,obs])%*%gamma)

}
mean_beta<-(1/mean_beta1)*mean_beta2
sigma_beta<-1/mean_beta1
beta2<-rnorm(1, mean = mean_beta, sd = sqrt((sigma_epsilon_squared)*sigma_beta))

beta[obs]<-beta2
result_betas[i,obs]<-beta2
}
if(i%%10==0)
{
print(i)
}

}