rivDP {bayesm}

Linear "IV" Model with DP Process Prior for Errors

Description

rivDP is a Gibbs Sampler for a linear structural equation with an arbitrary number of instruments. rivDP uses a mixture-of-normals for the structural and reduced form equations implemented with a Dirichlet Process prior.

Usage

rivDP(Data, Prior, Mcmc)

Arguments

Data

list(y, x, w, z)

Prior

list(md, Ad, mbg, Abg, lambda, Prioralpha, lambda_hyper)

Mcmc

list(R, keep, nprint, maxuniq, SCALE, gridsize)

Details

Model and Priors

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 non-zero 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-(alpha-alpha_{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 hyper-parameters.

\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]]

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, 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 Details

nmix is a list with 3 components. Several functions in the bayesm package that involve a Dirichlet Process or mixture-of-normals 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 one-column 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 inner-most lists has 2 elemens: a vector of draws for mu and a matrix of draws for the Cholesky root of Sigma.

Value

A list containing:

deltadraw

R/keep x p array of delta draws

betadraw

R/keep x 1 vector of beta draws

alphadraw

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

Istardraw

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

gammadraw

R/keep x j array of gamma draws

nmix

a list containing: probdraw, zdraw, compdraw (see “nmix Details” section)

Author(s)

Peter Rossi, Anderson School, UCLA, [e-mail].

References

For further discussion, see "A Semi-Parametric Bayesian Approach to the Instrumental Variable Problem," by Conley, Hansen, McCulloch and Rossi, Journal of Econometrics (2008).

See also, Chapter 4, Bayesian Non- and Semi-parametric Methods and Applications by Peter Rossi.

See Also

rivGibbs

Examples

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

## simulate scaled log-normal 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
}