rsurGibbs {bayesm}

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)

Details

Model and Priors

y_i = X_i\beta_i + e_i with i=1,\ldots,m for m regressions
(e(k,1), \ldots, e(k,m))' \sim N(0, \Sigma) with k=1, \ldots, n

Can be written as a stacked model:
y = X\beta + e where y is a nobs*m vector and p = length(beta) = sum(length(beta_i))

Note: must have the same number of observations (n) in each equation but can have a different number of X variables (p_i) for each equation where p = \sum p_i.

\beta \sim N(betabar, A^{-1})
\Sigma \sim IW(nu,V)

Argument Details

Data = list(regdata)

regdata:

list of lists, regdata[[i]] = list(y=y_i, X=X_i), where y_i is n x 1 and X_i is n x p_i

Prior = list(betabar, A, nu, V) [optional]

betabar:	p x 1 prior mean (def: 0)
A:	p x p prior precision matrix (def: 0.01*I)
nu:	d.f. parameter for Inverted Wishart prior (def: m+3)
V:	m x m scale parameter for Inverted Wishart prior (def: nu*I)

Mcmc = list(R, keep) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint:	print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)

Value

A list containing:

betadraw

R x p matrix of betadraws

Sigmadraw

R x (m*m) array of Sigma draws

Author(s)

Peter Rossi, Anderson School, UCLA, [e-mail].

References

For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

rmultireg

Examples

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10}
set.seed(66)

## simulate data from SUR
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(0.5, 0.2, 0.2, 0.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)

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

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

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