Linear "IV" Model with DP Process Prior for Errors
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.
rivDP(Data, Prior, Mcmc)
Data |
list(y, x, w, z) |
Prior |
list(md, Ad, mbg, Abg, lambda, Prioralpha, lambda_hyper) |
Mcmc |
list(R, keep, nprint, maxuniq, SCALE, gridsize) |
x = z'δ + e1
y = β*x + w'γ + e2
e1,e2 ~ N(θ_{i}) where θ_{i} represents μ_{i}, Σ_{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.
δ ~ N(md, Ad^{-1})
vec(β, γ) ~ N(mbg, Abg^{-1})
θ_{i} ~ G
G ~ DP(alpha, G_0)
alpha ~ (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 (μ,Σ):
Σ ~ IW(nu, vI) and μ|Σ ~ N(0, Σ(x) a^{-1})
These parameters are collected together in the list λ.
It is highly recommended that you use the default settings for these hyper-parameters.
λ(a, nu, v):
a ~ uniform[alim[1], alimb[2]]
nu ~ dim(data)-1 + exp(z)
z ~ uniform[dim(data)-1+nulim[1], nulim[2]]
v ~ uniform[vlim[1], vlim[2]]
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 Detailsnmix 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.
|
A list containing:
|
R/keep x p array of delta draws |
|
R/keep x 1 vector of beta draws |
|
R/keep x 1 vector of draws of Dirichlet Process tightness parameter |
|
R/keep x 1 vector of draws of the number of unique normal components |
|
R/keep x j array of gamma draws |
|
a list containing: |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
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.
rivGibbs
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
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.