rivDP: 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)

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)

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)

Author

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.

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`.

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.

Examples

Run this code

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
}

Run the code above in your browser using DataLab

`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

`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)`)

`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`.