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)
Data |
list(y, x, w, z) |
Prior |
list(md, Ad, mbg, Abg, nu, V) |
Mcmc |
list(R, keep, nprint) |
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)
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:
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 |
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.
http://www.perossi.org/home/bsm-1
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)}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.