HMMpanelFE | R Documentation |
HMMpanelFE generates a sample from the posterior distribution of the fixed-effects model with varying individual effects model discussed in Park (2011). The code works for both balanced and unbalanced panel data as long as there is no missing data in the middle of each group. This model uses a multivariate Normal prior for the fixed effects parameters and varying individual effects, an Inverse-Gamma 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.
HMMpanelFE(
subject.id,
y,
X,
m,
mcmc = 1000,
burnin = 1000,
thin = 1,
verbose = 0,
b0 = 0,
B0 = 0.001,
c0 = 0.001,
d0 = 0.001,
delta0 = 0,
Delta0 = 0.001,
a = NULL,
b = NULL,
seed = NA,
...
)
subject.id |
A numeric vector indicating the group number. It should start from 1. |
y |
The response variable. |
X |
The model matrix excluding the constant. |
m |
A vector of break numbers for each subject in the panel. |
mcmc |
The number of MCMC iterations after burn-in. |
burnin |
The number of burn-in 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 |
|
delta0 |
The prior mean of |
Delta0 |
The prior precision of |
a |
|
b |
|
seed |
The seed for the random number generator. If NA, current R system seed is used. |
... |
further arguments to be passed |
HMMpanelFE
simulates from the fixed-effect hidden Markov
pbject level:
\varepsilon_{it} \sim \mathcal{N}(\alpha_{im},
\sigma^2_{im})
We assume standard, semi-conjugate priors:
\beta \sim
\mathcal{N}(b_0,B_0^{-1})
And:
\sigma^{-2} \sim
\mathcal{G}amma(c_0/2, d_0/2)
And:
\alpha \sim
\mathcal{N}(delta_0,Delta_0^{-1})
\beta
, \alpha
and
\sigma^{-2}
are assumed a priori independent.
And:
p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M
Where M
is the number of states.
OLS estimates are used for starting 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 sigma
storage
matrix that contains time-varying residual variance, an attribute
state
storage matrix that contains posterior samples of
hidden states, and an attribute delta
storage matrix
containing time-varying intercepts.
Jong Hee Park, 2012. “Unified Method for Dynamic and Cross-Sectional Heterogeneity: Introducing Hidden Markov Panel Models.” American Journal of Political Science.56: 1040-1054. <doi: 10.1111/j.1540-5907.2012.00590.x>
Siddhartha Chib. 1998. “Estimation and comparison of multiple change-point models.” Journal of Econometrics. 86: 221-241. <doi: 10.1016/S0304-4076(97)00115-2>
Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. “MCMCpack: Markov Chain Monte Carlo in R.”, Journal of Statistical Software. 42(9): 1-21. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v042.i09")}.
## Not run:
## data generating
set.seed(1974)
N <- 30
T <- 80
NT <- N*T
## true parameter values
true.beta <- c(1, 1)
true.sigma <- 3
x1 <- rnorm(NT)
x2 <- runif(NT, 2, 4)
## group-specific breaks
break.point = rep(T/2, N); break.sigma=c(rep(1, N));
break.list <- rep(1, N)
X <- as.matrix(cbind(x1, x2), NT, );
y <- rep(NA, NT)
id <- rep(1:N, each=NT/N)
K <- ncol(X);
true.beta <- as.matrix(true.beta, K, 1)
## compute the break probability
ruler <- c(1:T)
W.mat <- matrix(NA, T, N)
for (i in 1:N){
W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
}
Weight <- as.vector(W.mat)
## draw time-varying individual effects and sample y
j = 1
true.sigma.alpha <- 30
true.alpha1 <- true.alpha2 <- rep(NA, N)
for (i in 1:N){
Xi <- X[j:(j+T-1), ]
true.mean <- Xi %*% true.beta
weight <- Weight[j:(j+T-1)]
true.alpha1[i] <- rnorm(1, 0, true.sigma.alpha)
true.alpha2[i] <- -1*true.alpha1[i]
y[j:(j+T-1)] <- ((1-weight)*true.mean + (1-weight)*rnorm(T, 0, true.sigma) +
(1-weight)*true.alpha1[i]) +
(weight*true.mean + weight*rnorm(T, 0, true.sigma) + weight*true.alpha2[i])
j <- j + T
}
## extract the standardized residuals from the OLS with fixed-effects
FEols <- lm(y ~ X + as.factor(id) -1 )
resid.all <- rstandard(FEols)
time.id <- rep(1:80, N)
## model fitting
G <- 100
BF <- testpanelSubjectBreak(subject.id=id, time.id=time.id,
resid= resid.all, max.break=3, minimum = 10,
mcmc=G, burnin = G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, Time = time.id)
## get the estimated break numbers
estimated.breaks <- make.breaklist(BF, threshold=3)
## model fitting
out <- HMMpanelFE(subject.id = id, y, X=X, m = estimated.breaks,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, delta0=0, Delta0=1/100)
## print out the slope estimate
## true values are 1 and 1
summary(out)
## compare them with the result from the constant fixed-effects
summary(FEols)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.