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'\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 |
|
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
}