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'\delta + e1
y = \beta*x + w'\gamma + e2
e1,e2
\sim
N(\theta_{i})
where \theta_{i}
represents \mu_{i}, \Sigma_{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.
\delta
\sim
N(md, Ad^{1})
vec(\beta, \gamma)
\sim
N(mbg, Abg^{1})
\theta_{i}
\sim
G
G
\sim
DP(alpha, G_0)
alpha
\sim
(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 (\mu,\Sigma)
:
\Sigma
\sim
IW(nu, vI)
and \mu\Sigma
\sim
N(0, \Sigma(x) a^{1})
These parameters are collected together in the list \lambda
.
It is highly recommended that you use the default settings for these hyperparameters.
\lambda(a, nu, v):
a
\sim
uniform[alim[1], alimb[2]]
nu
\sim
dim(data)1 + exp(z)
z
\sim
uniform[dim(data)1+nulim[1], nulim[2]]
v
\sim
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:  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) 
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:
deltadraw 

betadraw 

alphadraw 

Istardraw 

gammadraw 

nmix 
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
}