rivDP {bayesm}R Documentation

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δ+e1x = z'\delta + e1
y=βx+wγ+e2y = \beta*x + w'\gamma + e2
e1,e2e1,e2 \sim N(θi)N(\theta_{i}) where θi\theta_{i} represents μi,Σi\mu_{i}, \Sigma_{i}

Note: Error terms have non-zero means. DO NOT include intercepts in the zz or ww matrices. This is different from rivGibbs which requires intercepts to be included explicitly.

δ\delta \sim N(md,Ad1)N(md, Ad^{-1})
vec(β,γ)vec(\beta, \gamma) \sim N(mbg,Abg1)N(mbg, Abg^{-1})
θi\theta_{i} \sim GG
GG \sim DP(alpha,G0)DP(alpha, G_0)

alphaalpha \sim (1(alphaalphamin)/(alphamaxalphamin))power(1-(alpha-alpha_{min})/(alpha_{max}-alpha{min}))^{power}
where alphaminalpha_{min} and alphamaxalpha_{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.

G0G_0 is the natural conjugate prior for (μ,Σ)(\mu,\Sigma): Σ\Sigma \sim IW(nu,vI)IW(nu, vI) and μΣ\mu|\Sigma \sim N(0,Σ(x)a1)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.

λ(a,nu,v):\lambda(a, nu, v):
aa \sim uniform[alim[1], alimb[2]]
nunu \sim dim(data)-1 + exp(z)
zz \sim uniform[dim(data)-1+nulim[1], nulim[2]]
vv \sim uniform[vlim[1], vlim[2]]

Argument Details

Data = list(y, x, w, z)

y: nx1n x 1 vector of obs on LHS variable in structural equation
x: nx1n x 1 vector of obs on "endogenous" variable in structural equation
w: nxjn x j matrix of obs on "exogenous" variables in the structural equation
z: nxpn x p matrix of obs on instruments

Prior = list(md, Ad, mbg, Abg, lambda, Prioralpha, lambda_hyper) [optional]

md: pp-length prior mean of delta (def: 0)
Ad: pxpp x p PDS prior precision matrix for prior on delta (def: 0.01*I)
mbg: (j+1)(j+1)-length prior mean vector for prior on beta,gamma (def: 0)
Abg: (j+1)x(j+1)(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: ncompxR/keepncomp 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/keepxnobsR/keep x nobs matrix that indicates which component each draw is assigned to (here, null)
compdraw: A list of R/keepR/keep lists of ncompncomp 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/keepxpR/keep x p array of delta draws

betadraw

R/keepx1R/keep x 1 vector of beta draws

alphadraw

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

Istardraw

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

gammadraw

R/keepxjR/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, perossichi@gmail.com.

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
}

[Package bayesm version 3.1-6 Index]