rivDP {bayesm}  R Documentation 
rivDP
is a Gibbs Sampler for a linear structural equation with an arbitrary number of instruments. rivDP
uses a mixtureofnormals for the structural and reduced form equations implemented with a Dirichlet Process prior.
rivDP(Data, Prior, Mcmc)
Data 
list(y, x, w, z) 
Prior 
list(md, Ad, mbg, Abg, lambda, Prioralpha, lambda_hyper) 
Mcmc 
list(R, keep, nprint, maxuniq, SCALE, gridsize) 
x = z'δ + e1
y = β*x + w'γ + e2
e1,e2 ~ N(θ_{i}) where θ_{i} represents μ_{i}, Σ_{i}
Note: Error terms have nonzero means.
DO NOT include intercepts in the z or w matrices.
This is different from rivGibbs
which requires intercepts to be included explicitly.
δ ~ N(md, Ad^{1})
vec(β, γ) ~ N(mbg, Abg^{1})
θ_{i} ~ G
G ~ DP(alpha, G_0)
alpha ~ (1(alphaalpha_{min})/(alpha_{max}alpha{min}))^{power}
where alpha_{min} and alpha_{max} are set using the arguments in the reference below.
It is highly recommended that you use the default values for the hyperparameters of the prior on alpha.
G_0 is the natural conjugate prior for (μ,Σ):
Σ ~ IW(nu, vI) and μΣ ~ N(0, Σ(x) a^{1})
These parameters are collected together in the list λ.
It is highly recommended that you use the default settings for these hyperparameters.
λ(a, nu, v):
a ~ uniform[alim[1], alimb[2]]
nu ~ dim(data)1 + exp(z)
z ~ uniform[dim(data)1+nulim[1], nulim[2]]
v ~ uniform[vlim[1], vlim[2]]
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, lambda, Prioralpha, lambda_hyper)
[optional]
md:  plength 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) 
Prioralpha:  list(Istarmin, Istarmax, power) 
$Istarmin:  is expected number of components at lower bound of support of alpha (def: 1) 
$Istarmax:  is expected number of components at upper bound of support of alpha (def: floor(0.1*length(y)) ) 
$power:  is the power parameter for alpha prior (def: 0.8) 
lambda_hyper:  list(alim, nulim, vlim) 
$alim:  defines support of a distribution (def: c(0.01, 10) ) 
$nulim:  defines support of nu distribution (def: c(0.01, 3) ) 
$vlim:  defines support of v distribution (def: c(0.1, 4) )

Mcmc = list(R, keep, nprint, maxuniq, SCALE, gridsize)
[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) 
maxuniq:  storage constraint on the number of unique components (def: 200) 
SCALE:  scale data (def: TRUE ) 
gridsize:  gridsize parameter for alpha draws (def: 20) 
nmix
Detailsnmix
is a list with 3 components. Several functions in the bayesm
package that involve a Dirichlet Process or mixtureofnormals return nmix
. Across these functions, a common structure is used for nmix
in order to utilize generic summary and plotting functions.
probdraw:  ncomp x R/keep matrix that reports the probability that each draw came from a particular component (here, a onecolumn matrix of 1s) 
zdraw:  R/keep x nobs matrix that indicates which component each draw is assigned to (here, null) 
compdraw:  A list of R/keep lists of ncomp lists. Each of the innermost lists has 2 elemens: a vector of draws for mu and a matrix of draws for the Cholesky root of Sigma .

A list containing:

R/keep x p array of delta draws 

R/keep x 1 vector of beta draws 

R/keep x 1 vector of draws of Dirichlet Process tightness parameter 

R/keep x 1 vector of draws of the number of unique normal components 

R/keep x j array of gamma draws 

a list containing: 
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see "A SemiParametric Bayesian Approach to the Instrumental Variable Problem," by Conley, Hansen, McCulloch and Rossi, Journal of Econometrics (2008).
See also, Chapter 4, Bayesian Non and Semiparametric Methods and Applications by Peter Rossi.
rivGibbs
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} set.seed(66) ## simulate scaled lognormal errors and run k = 10 delta = 1.5 Sigma = matrix(c(1, 0.6, 0.6, 1), ncol=2) N = 1000 tbeta = 4 scalefactor = 0.6 root = chol(scalefactor*Sigma) mu = c(1,1) ## compute interquartile ranges ninterq = qnorm(0.75)  qnorm(0.25) error = matrix(rnorm(100000*2), ncol=2)%*%root error = t(t(error)+mu) Err = t(t(exp(error))exp(mu+0.5*scalefactor*diag(Sigma))) lnNinterq = quantile(Err[,1], prob=0.75)  quantile(Err[,1], prob=0.25) ## simulate data error = matrix(rnorm(N*2), ncol=2)%*%root error = t(t(error)+mu) Err = t(t(exp(error))exp(mu+0.5*scalefactor*diag(Sigma))) ## scale appropriately Err[,1] = Err[,1]*ninterq/lnNinterq Err[,2] = Err[,2]*ninterq/lnNinterq z = matrix(runif(k*N), ncol=k) x = z%*%(delta*c(rep(1,k))) + Err[,1] y = x*tbeta + Err[,2] ## specify data input and mcmc parameters Data = list(); Data$z = z Data$x = x Data$y = y Mcmc = list() Mcmc$maxuniq = 100 Mcmc$R = R end = Mcmc$R out = rivDP(Data=Data, Mcmc=Mcmc) cat("Summary of Beta draws", fill=TRUE) summary(out$betadraw, tvalues=tbeta) ## plotting examples if(0){ plot(out$betadraw, tvalues=tbeta) plot(out$nmix) # plot "fitted" density of the errors }