Structural Breaks ala Jong Hee Park

panel data
code
analysis
Author

Robert W. Walker

Published

July 25, 2026

Introduction

Overview

library(MCMCpack)
library(zoo)
## The original data set is from 
## "Bank Rate Policy Under the Interwar Gold Standard: A Dynamic
## Probit Model" by Barry Eichengreen, Mark Watson and Richard Grossman, The
## Economic Journal, Vol. 95, Spetember 1985, pp 725-745.
## Available at http://www.princeton.edu/~mwatson/publi.html
##############################################################################
load("./data/EWG.RData")
EWG$BRold <- BR[-c(1:2),1][-c(329:330)] ## Make it comparable to EWG time frame: May 28, 1925 - Sep 13, 1931
BRraw <- cumsum(EWG$dBR) + 5 ## restore raw data
BRord <- rep(NA, nrow(EWG)) ## ordinal data
for (i in 1:nrow(EWG)){
  if (BRraw[i]<3){
    BRord[i] <- 1
  }
  else if (BRraw[i]>=6){
    BRord[i] <- 5
  }
  else if (BRraw[i]>=5&BRraw[i]<=5.5){
    BRord[i] <- 4
  }
  else if (BRraw[i]>=4&BRraw[i]<=4.5){
    BRord[i] <- 3
  }
  else if (BRraw[i]>=3&BRraw[i]<=3.5){
    BRord[i] <- 2
  }
  else{
    BRord[i] <- NA
  }
}
EWG$BRord <- BRord

## centering
EWG$dGfall.cen <- (EWG$dGfall - mean(EWG$dGfall))/sd(EWG$dGfall)
EWG$dGrise.cen <- (EWG$dGrise - mean(EWG$dGrise))/sd(EWG$dGrise)
EWG$goldstock <- cumsum(EWG$dG)
EWG$gold.cen <- (EWG$goldstock - mean(EWG$goldstock))/sd(EWG$goldstock)
EWG$dummyBRhigh <- ifelse(EWG$BRhigh==0, 0, 1)
indicator <- ifelse(EWG$dBRMRplus==0, TRUE, FALSE)
dBRMRminus <- rep(0, length(BRord))
for (i in 1: (length(BRord)-1)){dBRMRminus[i+1] <- ifelse(diff(EWG$dBRMR)[i]<0, EWG$dBRMR[i+1], 0)}
EWG$dBRMRminus <- dBRMRminus
EWG$MR <- EWG$BRold - EWG$dBRMR ## Market rate
EWG$dMR <- c(NA, diff(EWG$MR))
EWG$time.week <- c(13:52, rep(1:52, 5), 1:28)
EWG$time.year <- c(rep(1925, length(13:52)), rep(1926:1930, each=52), rep(1931, length(1:28)))

## Parameter estimation
attach(EWG)
formula <- BRord ~ dGfall + dGrise + dBRMRminus + dBRMRplus + LDNY + dTB  +  Ex
K = 8; b0 <- rep(0, K); B0 <- diag(10, K)
N <- dim(EWG)[1]
m=1; ns = m + 1
detach(EWG)
change <- MCMCoprobitChange(formula, data = EWG, m=m,
                            b0=0, B0=10, a=floor(N/ns), b=1, tune=c(.1, .5), 
                            mcmc=10000, burnin=10000, thin=1, verbose=2000)


MCMCoprobitChange iteration 2001 of 20000 


 Acceptance rate for state 1 is    0.53400 


 Acceptance rate for state 2 is    0.41000 

 The number of observations in state 1 is      00249
 The number of observations in state 2 is      00079
 beta 0 =   -0.06818  -0.01613   0.22216   0.20938   0.46054   0.63366  -0.00943   0.49317
 beta 1 =    0.02029  -0.03101  -0.31909  -0.21273  -0.35007   0.15412   0.01074   0.26392
 gamma 0 =    0.21381   2.74269   4.75029
 gamma 1 =    2.38844  15.39177  28.40335

MCMCoprobitChange iteration 4001 of 20000 


 Acceptance rate for state 1 is    0.54125 


 Acceptance rate for state 2 is    0.43675 

 The number of observations in state 1 is      00250
 The number of observations in state 2 is      00078
 beta 0 =    0.25474  -0.20074   0.12962   0.54483   0.69213   0.74164  -0.01047   0.66294
 beta 1 =    0.32348  -0.01505  -0.11264   0.22983   0.27722   0.51752  -0.03006   0.08621
 gamma 0 =    0.62972   3.85288   6.05743
 gamma 1 =    3.13803  20.92688  29.06703

MCMCoprobitChange iteration 6001 of 20000 


 Acceptance rate for state 1 is    0.53283 


 Acceptance rate for state 2 is    0.44567 

 The number of observations in state 1 is      00250
 The number of observations in state 2 is      00078
 beta 0 =   -0.22012  -0.16413   0.14503   0.45529   0.41054   0.67066   0.00946   0.55512
 beta 1 =    0.43296   0.05311  -0.10654  -0.37047  -0.02500   0.56132   0.00075   0.22004
 gamma 0 =    0.22586   2.85715   4.87699
 gamma 1 =    3.34402  31.40030  42.94427

MCMCoprobitChange iteration 8001 of 20000 


 Acceptance rate for state 1 is    0.53250 


 Acceptance rate for state 2 is    0.45075 

 The number of observations in state 1 is      00249
 The number of observations in state 2 is      00079
 beta 0 =   -0.15679  -0.03670   0.12881   0.18843   0.44715   0.61289  -0.02784   0.61680
 beta 1 =    0.63001  -0.05512  -0.55185   0.10660   0.00383   0.60817  -0.00900   0.15347
 gamma 0 =    0.34370   3.05740   5.09322
 gamma 1 =    2.58508  14.42322  44.02439

MCMCoprobitChange iteration 10001 of 20000 


 Acceptance rate for state 1 is    0.53730 


 Acceptance rate for state 2 is    0.45470 

 The number of observations in state 1 is      00244
 The number of observations in state 2 is      00084
 beta 0 =    0.27029  -0.12249   0.14590   0.43006   0.65605   1.05815  -0.00261   0.47406
 beta 1 =    0.10488  -0.17043  -0.33159   0.22112   0.05060   0.59058  -0.04146   0.24535
 gamma 0 =    0.60558   3.06734   5.17529
 gamma 1 =    3.00555   8.50475  70.83479

MCMCoprobitChange iteration 12001 of 20000 


 Acceptance rate for state 1 is    0.53767 


 Acceptance rate for state 2 is    0.45733 

 The number of observations in state 1 is      00249
 The number of observations in state 2 is      00079
 beta 0 =    0.10047  -0.15432   0.19107   0.01510   0.31761   0.45437  -0.00700   0.53657
 beta 1 =    0.33612  -0.09438  -0.23707   0.32611   0.25640   0.08145  -0.02559   0.18483
 gamma 0 =    0.51287   2.94322   4.90706
 gamma 1 =    2.74068  25.74449  84.49184

MCMCoprobitChange iteration 14001 of 20000 


 Acceptance rate for state 1 is    0.53507 


 Acceptance rate for state 2 is    0.45750 

 The number of observations in state 1 is      00249
 The number of observations in state 2 is      00079
 beta 0 =    0.14809  -0.05322   0.04532  -0.05231   0.47895   0.67765  -0.01162   0.56080
 beta 1 =   -0.09841  -0.09466  -0.15649   0.14347   0.75264   0.46024  -0.06510   0.12063
 gamma 0 =    0.59130   3.19554   5.04141
 gamma 1 =    1.90567  17.18872 104.93629

MCMCoprobitChange iteration 16001 of 20000 


 Acceptance rate for state 1 is    0.53856 


 Acceptance rate for state 2 is    0.46006 

 The number of observations in state 1 is      00250
 The number of observations in state 2 is      00078
 beta 0 =    0.67525  -0.11530   0.04808   0.53810   0.35696   0.54432  -0.02480   0.49587
 beta 1 =   -0.11790  -0.20828  -0.29578   0.01350  -0.16170   0.25763  -0.05923   0.29803
 gamma 0 =    0.43011   3.33289   4.97938
 gamma 1 =    2.96491   7.65293  90.00330

MCMCoprobitChange iteration 18001 of 20000 


 Acceptance rate for state 1 is    0.53967 


 Acceptance rate for state 2 is    0.46044 

 The number of observations in state 1 is      00251
 The number of observations in state 2 is      00077
 beta 0 =   -0.03218  -0.20265   0.23399   0.71228   0.73318   1.00530  -0.04380   0.90458
 beta 1 =    0.57432  -0.08535  -0.22417  -0.35572  -0.44800   0.54517  -0.05727   0.10745
 gamma 0 =    1.19450   4.99270   7.58066
 gamma 1 =    3.80293  17.52985  93.18674
summary(change)

Iterations = 1:10000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                         Mean      SD  Naive SE Time-series SE
(Intercept)_regime1  0.150005 0.31142 0.0031142      0.0031142
dGfall_regime1      -0.122353 0.07842 0.0007842      0.0010708
dGrise_regime1       0.151983 0.09775 0.0009775      0.0013327
dBRMRminus_regime1   0.217226 0.23086 0.0023086      0.0028175
dBRMRplus_regime1    0.362069 0.21397 0.0021397      0.0026784
LDNY_regime1         0.685722 0.14389 0.0014389      0.0024178
dTB_regime1         -0.006561 0.02003 0.0002003      0.0002798
Ex_regime1           0.669490 0.14143 0.0014143      0.0288761
(Intercept)_regime2  0.048746 0.31014 0.0031014      0.0032567
dGfall_regime2      -0.067696 0.06066 0.0006066      0.0009786
dGrise_regime2      -0.237594 0.13278 0.0013278      0.0019201
dBRMRminus_regime2   0.071535 0.27227 0.0027227      0.0031498
dBRMRplus_regime2   -0.074425 0.26303 0.0026303      0.0032207
LDNY_regime2         0.422526 0.15105 0.0015105      0.0026108
dTB_regime2         -0.018966 0.03655 0.0003655      0.0005594
Ex_regime2           0.222082 0.08346 0.0008346      0.0012338
gamma3_regime1       0.674291 0.52808 0.0052808      0.0832459
gamma4_regime1       3.758298 0.63233 0.0063233      0.1617985
gamma5_regime1       5.890259 0.66608 0.0066608      0.1749004
gamma3_regime2       2.831611 0.33459 0.0033459      0.0104207
gamma4_regime2      16.393400 7.85131 0.0785131      3.4990751
gamma5_regime2      88.552898 9.00288 0.0900288      4.1647557

2. Quantiles for each variable:

                        2.5%      25%       50%       75%     97.5%
(Intercept)_regime1 -0.47745 -0.05504  0.147445  0.357830   0.75526
dGfall_regime1      -0.27589 -0.17496 -0.121287 -0.069609   0.03043
dGrise_regime1      -0.03867  0.08670  0.151193  0.218952   0.34429
dBRMRminus_regime1  -0.23296  0.06339  0.217104  0.373655   0.66689
dBRMRplus_regime1   -0.05071  0.21777  0.359418  0.504370   0.78443
LDNY_regime1         0.40464  0.58832  0.684372  0.781449   0.96698
dTB_regime1         -0.04548 -0.02015 -0.006487  0.007206   0.03296
Ex_regime1           0.42341  0.57202  0.656125  0.753577   0.97646
(Intercept)_regime2 -0.56913 -0.15807  0.048319  0.258299   0.65847
dGfall_regime2      -0.18917 -0.10737 -0.066272 -0.026188   0.04642
dGrise_regime2      -0.50094 -0.32665 -0.235467 -0.148030   0.01987
dBRMRminus_regime2  -0.45695 -0.11445  0.070466  0.252524   0.60669
dBRMRplus_regime2   -0.58493 -0.25091 -0.077081  0.106277   0.44058
LDNY_regime2         0.12699  0.31947  0.421561  0.523459   0.71583
dTB_regime2         -0.09020 -0.04379 -0.018942  0.006109   0.05255
Ex_regime2           0.05593  0.16613  0.222129  0.278746   0.38731
gamma3_regime1       0.05408  0.27768  0.524918  0.917466   2.03387
gamma4_regime1       2.71199  3.33171  3.678765  4.072113   5.18808
gamma5_regime1       4.76899  5.45195  5.800262  6.233119   7.41657
gamma3_regime2       2.21162  2.60493  2.818689  3.047037   3.53269
gamma4_regime2       4.73454  9.87092 15.363716 23.104061  30.58271
gamma5_regime2      72.15228 84.10519 89.105297 92.642352 107.03834

Plot the regimes

draw picture of the posterior regime probability

prob.state <- attr(change, "prob.state")    
require(zoo)
x.Date <- as.Date("1925-05-28") + 7*(0:328)
Pchange <- zoo(c(0,prob.state[,2]), x.Date)
newUNEMP <- zoo(unemp, x.Date)
y <- cbind(Pchange, newUNEMP)[-329,]
colnames(y) <- c("prob","unemp")
plot(Pchange)

plot(newUNEMP)

par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=.02)
plot(y$prob,  y$unemp, xlab="Probability of Regime 2", ylab="Unemployment", cex=2)