bayesm (version 2.0-9)

rsurGibbs: Gibbs Sampler for Seemingly Unrelated Regressions (SUR)

Description

rsurGibbs implements a Gibbs Sampler to draw from the posterior of the Seemingly Unrelated Regression (SUR) Model of Zellner

Usage

rsurGibbs(Data, Prior, Mcmc)

Arguments

Data

list(regdata)

Prior

list(betabar,A, nu, V)

Mcmc

list(R,keep)

Value

• list of MCMC draws
• betadraw
  R x k array of betadraws
• Sigmadraw
  R x (m*m) array of Sigma draws

concept

• bayes
• Gibbs Sampler
• regression
• SUR model
• Seemingly Unrelated Regression
• MCMC

Details

Model: $y_i = X_ibeta_i + e_i$. i=1,...,m. m regressions. (e(1,k), ..., e(m,k)) $\sim$ $N(0,Sigma)$. k=1, ..., nobs. We can also write as the stacked model: $y = Xbeta + e$ where y is a nobs*m long vector and k=length(beta)=sum(length(betai)). Note: we must have the same number of observations in each equation but we can have different numbers of X variables Priors: $beta$ $\sim$ $N(betabar,A^{-1})$. $Sigma$ $\sim$ $IW(nu,V)$. List arguments contain

• regdata

{list of lists, regdata[[i]]=list(y=yi,X=Xi)} betabar{k x 1 prior mean (def: 0)} A{k x k prior precision matrix (def: .01I)} nu{ d.f. parm for Inverted Wishart prior (def: m+3)} V{ scale parm for Inverted Wishart prior (def: nu*I)} R{ number of MCMC draws } keep{ thinning parameter - keep every keepth draw }

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 3. http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rmultireg

Examples

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10}
##
## simulate data from SUR
set.seed(66)
beta1=c(1,2)
beta2=c(1,-1,-2)
nobs=100
nreg=2
iota=c(rep(1,nobs))
X1=cbind(iota,runif(nobs))
X2=cbind(iota,runif(nobs),runif(nobs))
Sigma=matrix(c(.5,.2,.2,.5),ncol=2)
U=chol(Sigma)
E=matrix(rnorm(2*nobs),ncol=2)y1=X1%*%beta1+E[,1]
y2=X2%*%beta2+E[,2]
##
## run Gibbs Sampler
regdata=NULL
regdata[[1]]=list(y=y1,X=X1)
regdata[[2]]=list(y=y2,X=X2)
Mcmc=list(R=R)
out=rsurGibbs(Data=list(regdata=regdata),Mcmc=Mcmc)
##
## summarize GS output
if(R < 100) {begin=1} else {begin=100}
end=Mcmc$R
cat("betadraws ",fill=TRUE)
mat=apply(out$betadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(c(beta1,beta2),mat); rownames(mat)[1]="beta"; print(mat)
cat("Sigmadraws ",fill=TRUE)
mat=apply(out$Sigmadraw[begin:end,],2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="sigma"; print(mat)