rivGibbs {bayesm}R Documentation

Gibbs Sampler for Linear "IV" Model

Description

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

Usage

rivGibbs(Data, Prior, Mcmc)

Arguments

Data

list(y, x, w, z)

Prior

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

Mcmc

list(R, keep, nprint)

Details

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)

Value

A list containing:

deltadraw

R/keep x p matrix of delta draws

betadraw

R/keep x 1 vector of beta draws

gammadraw

R/keep x j matrix of gamma draws

Sigmadraw

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

Author(s)

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

References

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

Examples

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

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
if(0){plot(out$betadraw)}

[Package bayesm version 3.1-4 Index]