rsurGibbs: Gibbs Sampler for Seemingly Unrelated Regressions (SUR)

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/rsurgibbs_rcpp.r

Description

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

Usage

1
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β_i + e_i with i=1,…,m for m regressions
(e(k,1), …, e(k,m))' ~ N(0, Σ) with k=1, …, n

Can be written as a stacked model:
y = Xβ + 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 = ∑ p_i.

β ~ N(betabar, A^{-1})
Σ ~ 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, [email protected].

References

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

See Also

rmultireg

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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))}

Example output

 
Starting Gibbs Sampler for SUR Regression Model
  with  2  regressions
  and   100  observations for each regression
 
Prior Parms: 
betabar
[1] 0 0 0 0 0
A
     [,1] [,2] [,3] [,4] [,5]
[1,] 0.01 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00
[5,] 0.00 0.00 0.00 0.00 0.01
nu =  5
V = 
     [,1] [,2]
[1,]    5    0
[2,]    0    5
 
MCMC parms: 
R=  10  keep=  1  nprint=  100
 
 MCMC Iteration (est time to end - min) 
 Total Time Elapsed: 0.00 
Summary of beta draws
fewer than 100 draws submitted 
Summary of Sigmadraws
fewer than 100 draws submitted 

bayesm documentation built on July 21, 2017, 7:18 p.m.