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
where
represents
Note: Error terms have non-zero means.
DO NOT include intercepts in the or
matrices.
This is different from
rivGibbs
which requires intercepts to be included explicitly.
where and
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.
is the natural conjugate prior for
:
and
These parameters are collected together in the list .
It is highly recommended that you use the default settings for these hyper-parameters.
uniform[alim[1], alimb[2]]
dim(data)-1 + exp(z)
uniform[dim(data)-1+nulim[1], nulim[2]]
uniform[vlim[1], vlim[2]]
Argument Details
Data = list(y, x, w, z)
y: | vector of obs on LHS variable in structural equation |
x: | vector of obs on "endogenous" variable in structural equation |
w: | matrix of obs on "exogenous" variables in the structural equation |
z: | matrix of obs on instruments
|
Prior = list(md, Ad, mbg, Abg, lambda, Prioralpha, lambda_hyper)
[optional]
md: | -length prior mean of delta (def: 0) |
Ad: | PDS prior precision matrix for prior on delta (def: 0.01*I) |
mbg: | -length prior mean vector for prior on beta,gamma (def: 0) |
Abg: | 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: | matrix that reports the probability that each draw came from a particular component (here, a one-column matrix of 1s) |
zdraw: | matrix that indicates which component each draw is assigned to (here, null) |
compdraw: | A list of lists of 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 |
|
betadraw |
|
alphadraw |
|
Istardraw |
|
gammadraw |
|
nmix |
a list containing: |
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
}