rivGibbs {bayesm}R Documentation

Gibbs Sampler for Linear "IV" Model


rivGibbs is a Gibbs Sampler for a linear structural equation with an arbitrary number of instruments.


rivGibbs(Data, Prior, Mcmc)



list(y, x, w, z)


list(md, Ad, mbg, Abg, nu, V)


list(R, keep, nprint)


Model and Priors

x = z'δ + e1
y = β*x + w'γ + e2
e1,e2 ~ N(0, Σ)

Note: if intercepts are desired in either equation, include vector of ones in z or w

δ ~ N(md, Ad^{-1})
vec(β,γ) ~ N(mbg, Abg^{-1})
Σ ~ IW(nu, V)

Argument Details

Data = list(y, x, w, z)

y: n x 1 vector of obs on LHS variable in structural equation
x: n x 1 vector of obs on "endogenous" variable in structural equation
w: n x j matrix of obs on "exogenous" variables in the structural equation
z: n x p matrix of obs on instruments

Prior = list(md, Ad, mbg, Abg, nu, V) [optional]

md: p-length prior mean of delta (def: 0)
Ad: p x p PDS prior precision matrix for prior on delta (def: 0.01*I)
mbg: (j+1)-length prior mean vector for prior on beta,gamma (def: 0)
Abg: (j+1)x(j+1) PDS prior precision matrix for prior on beta,gamma (def: 0.01*I)
nu: d.f. parameter for Inverted Wishart prior on Sigma (def: 5)
V: 2 x 2 location matrix for Inverted Wishart prior on Sigma (def: nu*I)

Mcmc = list(R, keep, nprint) [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)


A list containing:


R/keep x p matrix of delta draws


R/keep x 1 vector of beta draws


R/keep x j matrix of gamma draws


R/keep x 4 matrix of Sigma draws – each row is the vector form of Sigma


Rob McCulloch and Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.


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


if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

simIV = function(delta, beta, Sigma, n, z, w, gamma) {
  eps = matrix(rnorm(2*n),ncol=2) %*% chol(Sigma)
  x = z%*%delta + eps[,1]
  y = beta*x + eps[,2] + w%*%gamma
  list(x=as.vector(x), y=as.vector(y)) 
n = 200
p=1 # number of instruments
z = cbind(rep(1,n), matrix(runif(n*p),ncol=p))
w = matrix(1,n,1)
rho = 0.8
Sigma = matrix(c(1,rho,rho,1), ncol=2)
delta = c(1,4)
beta = 0.5
gamma = c(1)
simiv = simIV(delta, beta, Sigma, n, z, w, gamma)

Data1 = list(); Data1$z = z; Data1$w=w; Data1$x=simiv$x; Data1$y=simiv$y
Mcmc1=list(); Mcmc1$R = R; Mcmc1$keep=1

out = rivGibbs(Data=Data1, Mcmc=Mcmc1)

cat("Summary of Beta draws", fill=TRUE)
summary(out$betadraw, tvalues=beta)

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

## plotting examples

[Package bayesm version 3.1-4 Index]