In jongheepark/BridgeChange: Bayesian Change-point models with Bridge prior for analyzing time-series and panel data.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

BridgeChange

R package BridgeChange constains functions useful to analyze time-series data and panel data with possibly large number of covariates and change-points. This package offers a Bayesian inference on the linear regression mode under high-dimensinal covariates whose effects on the outcome are allowed to be heterogeneous across time. This package will be useful to discover a set of predictive variables under dynamic setting where time-varying effect is expected to exist in many cases.

Installation

You can install the most recent version of BridgeChange from Gitub using the devtools package.

# install BridgeChange from Github
# you might need to install "devtools"
devtools::install_github("jongheepark/BridgeChange")

Alvarez et al. data

# load packages 
library(BridgeChange)
library("pcse")
set.seed(1999)
data("agl")
data = agl

Fit the one-way (time) fixed effect model using plm

The following two model arguments are availabel for the fixed-effects HMBB. - the fixed effects model ("within"), the default, - the pooling model ("pooling"),

We will use the one-way (time) fixed effect model because the original model has a time-invariant covariate (central).

model = "within"
index = c('country', 'year')
effect = 'time'
formula <- growth ~ lagg1 + opengdp + openex + openimp + leftc * central
pdata   <- pdata.frame(data, index)
pm <- plm(formula, data = pdata, model = model, effect = effect)
summary(pm)

Inspect the model (response and residuals)

We take a look at the response data and panel residuals to check the sign of misfit due to time-varying effects.

## response
plot(pdata$growth)

## panel residuals
coplot(pm$residuals ~ pdata[,index[2]]|pdata[,index[1]], data=pdata, ## number=length(unique(pdata[,index[1]])),
       overlap=.1, col="brown", type="l", 
       panel = panel.smooth, xlab="panel residuals by group and time")

Fitting HMBB

We fit three HMBBs with no break, one break, and two breaks to see whether effects of covariates change over time. To save time, we set mcmc = 100 here. BridgeFixedPanel transforms the panel data using plm arguments of model, effect and index first. Then, it fits HMBB on the transformed data. n.break sets the number of break to be estimated.

mcmc = 100; burn = 100; verbose = 100; thin = 1;
formula <- growth ~ lagg1 + opengdp + openex + openimp + leftc + central + inter
agl.cp0 <- BridgeFixedPanel(formula=formula, data = data, 
                            model = model, index = index, effect = effect,
                            mcmc=mcmc, verbose=verbose, Waic = TRUE, 
                            n.break = 0)
agl.cp1 <- BridgeFixedPanel(formula=formula, data = data, 
                            model = model, index = index, effect = effect,
                            mcmc=mcmc, verbose=verbose, Waic = TRUE, 
                            n.break = 1)
agl.cp2 <- BridgeFixedPanel(formula=formula, data = data, 
                            model = model, index = index, effect = effect,
                            mcmc=mcmc, verbose=verbose, Waic = TRUE, 
                            n.break = 2)

After fitting multiple models, we can compare their model-fits using WAIC. WaicCompare shows WAIC scores for a list of models. plotWaic() draws a plot of WAIC scores.

## model selection by WAIC
waic <- WaicCompare(list(agl.cp0, agl.cp1, agl.cp2), print = TRUE)
plotWaic(waic)

In addition to WAIC, we compare transitions of hidden states. plotState in MCMCpack can be used to draw hidden state transitions.

## state changes
par(mfrow=c(1, 2))
plotState(agl.cp1, start=1970, legend.control =c(1970, 0.85), main="One break")
plotState(agl.cp2, start=1970, legend.control =c(1970, 0.85), main="Two breaks")

The one break model looks good. We check the time-varying movements of the one break model using dotplotRegime(). Colors in dotplotRegime() are determined by the size of coefficients in the first regime. Red means positive, blue means negative, and grey means close to 0 in the first regime.

## all covariates
dotplotRegime(agl.cp1, hybrid=FALSE, start = 1970, location.bar=12, x.location="default",
              text.cex=0.8, main="Time-varying Movements of All Covariates")
## label as a legend
dotplotRegime(agl.cp1, hybrid=FALSE, start = 1970, location.bar=12, x.location="legend",
              text.cex=0.8, main="Time-varying Movements of All Covariates")

We visualize the movement of the selected covariate using select argument in dotplotRegime(). Here we choose the left-party government-related variable.

## leftc only
## select works like grep()
dotplotRegime(agl.cp1, hybrid=FALSE, start = 1970, location.bar=12, x.location="static",
              text.cex=0.8, select="left", main=("Left party-related covariates"))

jongheepark/BridgeChange documentation built on Jan. 12, 2023, 4:48 p.m.