A Test for the Group-level Break using a Multivariate Linear Regression Model with Breaks
testpanelGroupBreak fits a multivariate linear regression model with parametric breaks using panel residuals to test the existence of group-level breaks in panel residuals. The details are discussed in Park (2011).
testpanelGroupBreak(
subject.id,
time.id,
resid,
m = 1,
mcmc = 1000,
burnin = 1000,
thin = 1,
verbose = 0,
b0,
B0,
c0,
d0,
a = NULL,
b = NULL,
seed = 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. |
resid |
A vector of panel residuals |
m |
The number of changepoints. |
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. |
marginal.likelihood |
How should the marginal likelihood be calculated?
Options are: |
... |
further arguments to be passed |
testpanelGroupBreak fits a multivariate linear regression model with
parametric breaks using panel residuals to detect the existence of
system-level breaks in unobserved factors as discussed in Park (2011).
The model takes the following form:
e_{i} \sim \mathcal{N}(β_{m}, σ^2_m I)\;\; m = 1, …, M
We assume standard, semi-conjugate priors:
β \sim \mathcal{N}(b0, B0)
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.
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 log-marginal
likelihood of the model (logmarglike).
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:
## data generating
set.seed(1977)
Q <- 3
true.beta1 <- c(1, 1, 1) ; true.beta2 <- c(1, -1, -1)
true.sigma2 <- c(1, 3); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
N=20; T=100;
NT <- N*T
x1 <- rnorm(NT)
x2 <- runif(NT, 5, 10)
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((ruler-break.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+T-1), ]
Wi <- W[j:(j+T-1), ]
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+T-1)]
y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*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))
resid <- rstandard(lm(y ~X-1 + as.factor(subject.id)))
G <- 100
out0 <- testpanelGroupBreak(subject.id, time.id, resid, m=0,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
out1 <- testpanelGroupBreak(subject.id, time.id, resid, m=1,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
out2 <- testpanelGroupBreak(subject.id, time.id, resid, m=2,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
out3 <- testpanelGroupBreak(subject.id, time.id, resid, m=3,
mcmc=G, burnin=G, thin=1, verbose=G,
b0=0, B0=1/100, c0=2, d0=2, marginal.likelihood = "Chib95")
## Note that the code is for a hypothesis test of no break in panel residuals.
## When breaks exist, the estimated number of break in the mean and variance of panel residuals
## tends to be larger than the number of break in the data generating process.
## This is due to the difference in parameter space, not an error of the code.
BayesFactor(out0, out1, out2, out3)
## In order to identify the number of breaks in panel parameters,
## use HMMpanelRE() instead.
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.