A Test for the Subject-level Break using a Unitivariate Linear Regression Model with Breaks
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, ... )
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. |
resid |
A vector of panel residuals. |
max.break |
An upper bound of break numbers for the test. |
minimum |
A minimum length of time series for the test. The test will skip a subject with a time series shorter than this. |
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 the residual mean. |
B0 |
The prior precision of the residual variance |
c0 |
c_0/2 is the shape parameter for the inverse Gamma prior on σ^2. The amount of information in the inverse Gamma prior is something like that from c_0 pseudo-observations. |
d0 |
d_0/2 is the scale parameter for the inverse Gamma prior on σ^2. |
a |
a is the shape1 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states. |
b |
b is the shape2 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states. |
seed |
The seed for the random number generator. If NA, current R system seed is used. |
Time |
Times of the observations. This will be used to find the time of the first observations in panel residuals. |
ps.out |
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 |
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:
e_{it} = α_{im} + \varepsilon_{it}\;\; m = 1, …, M
The errors are assumed to be time-varying at the subject level:
\varepsilon_{it} \sim \mathcal{N}(0, σ^2_{im})
We assume standard, semi-conjugate priors:
β \sim \mathcal{N}(b_0,B_0^{-1})
And:
σ^{-2} \sim \mathcal{G}amma(c_0/2, d_0/2)
Where β and σ^{-2} are assumed a priori independent.
And:
p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, …, M
Where M is the number of states.
OLS estimates are used for starting values.
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.
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.1080/01621459.1995.10476635>
## 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"))
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.