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


2 Answers 2


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.

  • $\begingroup$ 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 ! $\endgroup$
    – quant
    Jun 21, 2016 at 9:04
  • $\begingroup$ 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. $\endgroup$
    – Hobbes
    Jun 21, 2016 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 :)


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


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

for (i in 1:niter) 
  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)))

  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
  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])
  gamma<-mvrnorm(n = 1, mu=mean_gamma, Sigma=sqrt(sigma_eta_squared)*sigma_gamma)


  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 )


  mean_beta1<-0 #this is scalar
  mean_beta2<-0 #this is scalar
  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)

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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.