testpanelSubjectBreak fits a unitivariate linear regression model with parametric breaks using panel residuals to test the existence of subject-level breaks in panel residuals. The details are discussed in Park (2011).
testpanelSubjectBreak(subject.id, time.id, resid, max.break = 2,
minimum = 10, mcmc = 1000, burnin = 1000, thin = 1, verbose = 0, b0,
B0, c0, d0, a = NULL, b = NULL, seed = NA, Time = NULL,
ps.out = FALSE)
A numeric vector indicating the group number. It should start from 1.
A numeric vector indicating the time unit. It should start from 1.
A vector of panel residuals.
An upper bound of break numbers for the test.
A minimum length of time series for the test. The test will skip a subject with a time series shorter than this.
The number of MCMC iterations after burn-in.
The number of burn-in iterations for the sampler.
The thinning interval used in the simulation. The number of MCMC iterations must be divisible by this value.
A switch which determines whether or not the progress of the
sampler is printed to the screen. If verbose
is greater than 0, the
iteration number and the posterior density samples are printed to the screen
every verbose
th iteration.
The prior mean of the residual mean.
The prior precision of the residual variance
The seed for the random number generator. If NA, current R system seed is used.
Times of the observations. This will be used to find the time of the first observations in panel residuals.
If ps.out == TRUE, state probabilities are exported. If the number of panel subjects is huge, users can turn it off to save memory.
further arguments to be passed
The returned object is a matrix containing log marginal likelihoods
for all HMMs. The dimension of the returned object is the number of panel
subjects by max.break + 1. If psout == TRUE, the returned object has an
array attribute psout
containing state probabilities for all HMMs.
testpanelSubjectBreak
fits a univariate linear regression model for
subject-level residuals from a panel model. The details are discussed in
Park (2011).
The model takes the following form:
The errors are assumed to be time-varying at the subject level:
We assume standard, semi-conjugate priors:
And:
Where
And:
Where
OLS estimates are used for starting values.
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.
Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.'' Journal of Econometrics. 86: 221-241.
# NOT RUN {
# }
# NOT RUN {
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 <- 1000
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)
## estimated break numbers
## thresho
estimated.breaks <- make.breaklist(BF, threshold=3)
## print all posterior model probabilities
print(attr(BF, "model.prob"))
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab