HMMpanelRE  R Documentation 
HMMpanelRE generates a sample from the posterior distribution of the hidden Markov randomeffects model discussed in Park (2011). The code works for panel data with the same starting point. The sampling of panel parameters is based on Algorithm 2 of Chib and Carlin (1999). This model uses a multivariate Normal prior for the fixed effects parameters and varying individual effects, an InverseWishart prior on the randomeffects parameters, an InverseGamma prior on the residual error variance, and Beta prior for transition probabilities. The user supplies data and priors, and a sample from the posterior distribution is returned as an mcmc object, which can be subsequently analyzed with functions provided in the coda package.
HMMpanelRE(
subject.id,
time.id,
y,
X,
W,
m = 1,
mcmc = 1000,
burnin = 1000,
thin = 1,
verbose = 0,
b0 = 0,
B0 = 0.001,
c0 = 0.001,
d0 = 0.001,
r0,
R0,
a = NULL,
b = NULL,
seed = NA,
beta.start = NA,
sigma2.start = NA,
D.start = NA,
P.start = NA,
marginal.likelihood = c("none", "Chib95"),
...
)
subject.id 
A numeric vector indicating the group number. It should start from 1. 
time.id 
A numeric vector indicating the time unit. It should start from 1. 
y 
The dependent variable 
X 
The model matrix of the fixedeffects 
W 
The model matrix of the randomeffects. W should be a subset of X. 
m 
The number of changepoints. 
mcmc 
The number of MCMC iterations after burnin. 
burnin 
The number of burnin iterations for the sampler. 
thin 
The thinning interval used in the simulation. The number of MCMC iterations must be divisible by this value. 
verbose 
A switch which determines whether or not the
progress of the sampler is printed to the screen. If

b0 
The prior mean of 
B0 
The prior precision of 
c0 

d0 

r0 
The shape parameter for the InverseWishart prior on variance matrix for the random effects. Set r=q for an uninformative prior where q is the number of random effects 
R0 
The scale matrix for the InverseWishart prior on variance matrix for the random effects. This must be a square qdimension matrix. Use plausible variance regarding random effects for the diagonal of R. 
a 

b 

seed 
The seed for the random number generator. If NA, current R system seed is used. 
beta.start 
The starting values for the beta vector. This can either be a scalar or a column vector with dimension equal to the number of betas. The default value of NA will use draws from the Uniform distribution with the same boundary with the data as the starting value. If this is a scalar, that value will serve as the starting value mean for all of the betas. When there is no covariate, the log value of means should be used. 
sigma2.start 
The starting values for 
D.start 
The starting values for the beta vector. This can either be a scalar or a column vector with dimension equal to the number of betas. The default value of NA will use draws from the Uniform distribution with the same boundary with the data as the starting value. If this is a scalar, that value will serve as the starting value mean for all of the betas. When there is no covariate, the log value of means should be used. 
P.start 
The starting values for the transition matrix. A
user should provide a square matrix with dimension equal to the
number of states. By default, draws from the 
marginal.likelihood 
How should the marginal likelihood be
calculated? Options are: 
... 
further arguments to be passed 
HMMpanelRE
simulates from the randomeffect hidden Markov
panel model introduced by Park (2011).
The model takes the following form:
y_i = X_i \beta_m + W_i
b_i + \varepsilon_i\;\; m = 1, \ldots, M
Where each group i
have
k_i
observations. Randomeffects parameters are assumed to
be timevarying at the system level:
b_i \sim
\mathcal{N}_q(0, D_m)
\varepsilon_i \sim \mathcal{N}(0,
\sigma^2_m I_{k_i})
And the errors: We assume standard, conjugate priors:
\beta
\sim \mathcal{N}_p(b0, B0)
And:
\sigma^{2} \sim
\mathcal{IG}amma(c0/2, d0/2)
And:
D \sim
\mathcal{IW}ishart(r0, R0)
See Chib and Carlin (1999) for more details.
And:
p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M
Where M
is the number of states.
NOTE: We do not provide default parameters for the priors on
the precision matrix for the random effects. When fitting one of
these models, it is of utmost importance to choose a prior that
reflects your prior beliefs about the random effects. Using the
dwish
and rwish
functions might be useful in choosing
these values.
An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda
package. The object contains an attribute prob.state
storage matrix that contains the probability of state_i
for
each period, and the logmarginal likelihood of the model
(logmarglike
).
Jong Hee Park, 2012. “Unified Method for Dynamic and CrossSectional Heterogeneity: Introducing Hidden Markov Panel Models.” American Journal of Political Science.56: 10401054. <doi: 10.1111/j.15405907.2012.00590.x>
Siddhartha Chib. 1998. “Estimation and comparison of multiple changepoint models.” Journal of Econometrics. 86: 221241. <doi: 10.1016/S03044076(97)001152>
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. “MCMCpack: Markov Chain Monte Carlo in R.”, Journal of Statistical Software. 42(9): 121. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v042.i09")}.
## Not run:
## data generating
set.seed(1977)
Q < 3
true.beta1 < c(1, 1, 1) ; true.beta2 < c(1, 1, 1)
true.sigma2 < c(2, 5); true.D1 < diag(.5, Q); true.D2 < diag(2.5, Q)
N=30; T=100;
NT < N*T
x1 < runif(NT, 1, 2)
x2 < runif(NT, 1, 2)
X < cbind(1, x1, x2); W < X; y < rep(NA, NT)
## true break numbers are one and at the center
break.point = rep(T/2, N); break.sigma=c(rep(1, N));
break.list < rep(1, N)
id < rep(1:N, each=NT/N)
K < ncol(X);
ruler < c(1:T)
## compute the weight for the break
W.mat < matrix(NA, T, N)
for (i in 1:N){
W.mat[, i] < pnorm((rulerbreak.point[i])/break.sigma[i])
}
Weight < as.vector(W.mat)
## data generating by weighting two means and variances
j = 1
for (i in 1:N){
Xi < X[j:(j+T1), ]
Wi < W[j:(j+T1), ]
true.V1 < true.sigma2[1]*diag(T) + Wi%*%true.D1%*%t(Wi)
true.V2 < true.sigma2[2]*diag(T) + Wi%*%true.D2%*%t(Wi)
true.mean1 < Xi%*%true.beta1
true.mean2 < Xi%*%true.beta2
weight < Weight[j:(j+T1)]
y[j:(j+T1)] < (1weight)*true.mean1 + (1weight)*chol(true.V1)%*%rnorm(T) +
weight*true.mean2 + weight*chol(true.V2)%*%rnorm(T)
j < j + T
}
## model fitting
subject.id < c(rep(1:N, each=T))
time.id < c(rep(1:T, N))
## model fitting
G < 100
b0 < rep(0, K) ; B0 < solve(diag(100, K))
c0 < 2; d0 < 2
r0 < 5; R0 < diag(c(1, 0.1, 0.1))
subject.id < c(rep(1:N, each=T))
time.id < c(rep(1:T, N))
out1 < HMMpanelRE(subject.id, time.id, y, X, W, m=1,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=b0, B0=B0, c0=c0, d0=d0, r0=r0, R0=R0)
## latent state changes
plotState(out1)
## print mcmc output
summary(out1)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.