Gibbs sampling in R

Question

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

Answer 1

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.

Answer 2

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)
  }

}