bayesm (version 3.0-2)

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

Details

Model: \(y_i = X_i\beta_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 = X\beta + 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

• regdatalist of lists, regdata[[i]]=list(y=yi,X=Xi)

• betabark x 1 prior mean (def: 0)

• Ak 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

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 3. http://www.perossi.org/home/bsm-1

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)%*%U
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)

Mcmc1=list(R=R)

out=rsurGibbs(Data=list(regdata=regdata),Mcmc=Mcmc1)

cat("Summary of beta draws",fill=TRUE)
summary(out$betadraw,tvalues=c(beta1,beta2))
cat("Summary of Sigmadraws",fill=TRUE)
summary(out$Sigmadraw,tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)]))

if(0){
plot(out$betadraw,tvalues=c(beta1,beta2))
}