Statistical Modeling I - Bayesian Synthetic Control With Different Latent Factor Models on Germany Reunification
Zhongming Jiang
2023-07-30
1 Pre-processing the Germany Reunification data
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())d <-
read.dta("/Users/apple/Desktop/Path\ Towards\ Quant\ Mkt\ PhD/Collected\ Data/Germany\ Reunification/repgermany.dta")
df_avg <- d |> group_by(index, country) |>
summarize_at(c("gdp", "infrate", "trade", "industry",
"schooling", "invest70", "invest80"),
mean, na.rm = TRUE)
df_avg <- mutate(df_avg, invest80 = (1/100)*invest80)
df_avg <- mutate(df_avg,
invest75 = mean(c_across(c("invest70", "invest80")), na.rm = TRUE))
df_J_P <- df_avg[, c(3:7, 10)]
rownames(df_J_P) <- 1:17## Warning: Setting row names on a tibble is deprecated.colnames(df_J_P) <-
c("avg_gdp", "avg_infrate", "avg_trade", "avg_industry", "avg_schooling", "avg_invest75")
german_unification <-
read.csv("/Users/apple/Desktop/Path\ Towards\ Quant\ Mkt\ PhD/Collected\ Data/Germany\ Reunification/Numpyro\ BGSC/german_unification.csv")
df_J_T <- german_unification[, 2:45]
colnames(df_J_T) <- 1960:2003
rownames(df_J_T) <- 1:172 Pinkney (2021)’s Data Augmentation Approach
Bayesian Synthetic Control Method with Interactive Fixed Effect Model
T <- length(unique(d$year))
J <- length(unique(d$country))
L <- 8
P <- 6
X <- data.matrix(t(df_J_P))
Y <- data.matrix(df_J_T)
trt_times <- max(d$year) - 1990
data_list <- list(T = T, J = J, L = L, P = P, X = X, Y = Y, trt_times = trt_times)
control_list <- list(max_treedepth = 14, adapt_delta = 0.95)
fit <- stan(file = "/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ I/Germany\ Reunification\ BSCM\ Derivation\ of\ ATT\ Estimator/BSCM_data_augmentation.stan",
data = data_list, init = "0.1",
control = control_list,
chains = 4, warmup = 500, iter = 1000)scm_bg <- df_final[, -(7:8)]
diff_GSCM <- scm_bg$real_gdp_WG - scm_bg$synth_gdp_WG_GSCM
diff_BSCM <- scm_bg$real_gdp_WG - scm_bg$synth_gdp_WG_BSCM
diff_BSCM_lwr <- scm_bg$real_gdp_WG - scm_bg$synth_gdp_WG_cred_int_2.5
diff_BSCM_upr <- scm_bg$real_gdp_WG - scm_bg$synth_gdp_WG_cred_int_97.5
time <- 1960:2003
scm_bg_diff <- data.frame(time, diff_GSCM, diff_BSCM,
diff_BSCM_lwr, diff_BSCM_upr)
SCM <- full_join(scm_bg, scm_bg_diff, by = "time")3 Pang, Liu, Xu (2021)’s Bayesian Causal Panel Analysis Approach
Bayesian Dynamic Multilevel Linear Factor Model
d <- mutate(d,
D_i = ifelse(index == 7 & year > 1990, 1, 0),
trt = ifelse(index == 7, 1, 0))
cols_to_convert <- c("invest60", "invest70")
for (col in cols_to_convert){
d[[col]] <- d[[col]] * 100
}
mice_mod <- mice(d, m = 5, maxit = 50, method = "pmm", seed = 500)##
## iter imp variable
## 1 1 infrate trade schooling invest60 invest70 invest80 industry
## 1 2 infrate trade schooling invest60 invest70 invest80 industry
## 1 3 infrate trade schooling invest60 invest70 invest80 industry
## 1 4 infrate trade schooling invest60 invest70 invest80 industry
## 1 5 infrate trade schooling invest60 invest70 invest80 industry
## 2 1 infrate trade schooling invest60 invest70 invest80 industry
## 2 2 infrate trade schooling invest60 invest70 invest80 industry
## 2 3 infrate trade schooling invest60 invest70 invest80 industry
## 2 4 infrate trade schooling invest60 invest70 invest80 industry
## 2 5 infrate trade schooling invest60 invest70 invest80 industry
## 3 1 infrate trade schooling invest60 invest70 invest80 industry
## 3 2 infrate trade schooling invest60 invest70 invest80 industry
## 3 3 infrate trade schooling invest60 invest70 invest80 industry
## 3 4 infrate trade schooling invest60 invest70 invest80 industry
## 3 5 infrate trade schooling invest60 invest70 invest80 industry
## 4 1 infrate trade schooling invest60 invest70 invest80 industry
## 4 2 infrate trade schooling invest60 invest70 invest80 industry
## 4 3 infrate trade schooling invest60 invest70 invest80 industry
## 4 4 infrate trade schooling invest60 invest70 invest80 industry
## 4 5 infrate trade schooling invest60 invest70 invest80 industry
## 5 1 infrate trade schooling invest60 invest70 invest80 industry
## 5 2 infrate trade schooling invest60 invest70 invest80 industry
## 5 3 infrate trade schooling invest60 invest70 invest80 industry
## 5 4 infrate trade schooling invest60 invest70 invest80 industry
## 5 5 infrate trade schooling invest60 invest70 invest80 industry
## 6 1 infrate trade schooling invest60 invest70 invest80 industry
## 6 2 infrate trade schooling invest60 invest70 invest80 industry
## 6 3 infrate trade schooling invest60 invest70 invest80 industry
## 6 4 infrate trade schooling invest60 invest70 invest80 industry
## 6 5 infrate trade schooling invest60 invest70 invest80 industry
## 7 1 infrate trade schooling invest60 invest70 invest80 industry
## 7 2 infrate trade schooling invest60 invest70 invest80 industry
## 7 3 infrate trade schooling invest60 invest70 invest80 industry
## 7 4 infrate trade schooling invest60 invest70 invest80 industry
## 7 5 infrate trade schooling invest60 invest70 invest80 industry
## 8 1 infrate trade schooling invest60 invest70 invest80 industry
## 8 2 infrate trade schooling invest60 invest70 invest80 industry
## 8 3 infrate trade schooling invest60 invest70 invest80 industry
## 8 4 infrate trade schooling invest60 invest70 invest80 industry
## 8 5 infrate trade schooling invest60 invest70 invest80 industry
## 9 1 infrate trade schooling invest60 invest70 invest80 industry
## 9 2 infrate trade schooling invest60 invest70 invest80 industry
## 9 3 infrate trade schooling invest60 invest70 invest80 industry
## 9 4 infrate trade schooling invest60 invest70 invest80 industry
## 9 5 infrate trade schooling invest60 invest70 invest80 industry
## 10 1 infrate trade schooling invest60 invest70 invest80 industry
## 10 2 infrate trade schooling invest60 invest70 invest80 industry
## 10 3 infrate trade schooling invest60 invest70 invest80 industry
## 10 4 infrate trade schooling invest60 invest70 invest80 industry
## 10 5 infrate trade schooling invest60 invest70 invest80 industry
## 11 1 infrate trade schooling invest60 invest70 invest80 industry
## 11 2 infrate trade schooling invest60 invest70 invest80 industry
## 11 3 infrate trade schooling invest60 invest70 invest80 industry
## 11 4 infrate trade schooling invest60 invest70 invest80 industry
## 11 5 infrate trade schooling invest60 invest70 invest80 industry
## 12 1 infrate trade schooling invest60 invest70 invest80 industry
## 12 2 infrate trade schooling invest60 invest70 invest80 industry
## 12 3 infrate trade schooling invest60 invest70 invest80 industry
## 12 4 infrate trade schooling invest60 invest70 invest80 industry
## 12 5 infrate trade schooling invest60 invest70 invest80 industry
## 13 1 infrate trade schooling invest60 invest70 invest80 industry
## 13 2 infrate trade schooling invest60 invest70 invest80 industry
## 13 3 infrate trade schooling invest60 invest70 invest80 industry
## 13 4 infrate trade schooling invest60 invest70 invest80 industry
## 13 5 infrate trade schooling invest60 invest70 invest80 industry
## 14 1 infrate trade schooling invest60 invest70 invest80 industry
## 14 2 infrate trade schooling invest60 invest70 invest80 industry
## 14 3 infrate trade schooling invest60 invest70 invest80 industry
## 14 4 infrate trade schooling invest60 invest70 invest80 industry
## 14 5 infrate trade schooling invest60 invest70 invest80 industry
## 15 1 infrate trade schooling invest60 invest70 invest80 industry
## 15 2 infrate trade schooling invest60 invest70 invest80 industry
## 15 3 infrate trade schooling invest60 invest70 invest80 industry
## 15 4 infrate trade schooling invest60 invest70 invest80 industry
## 15 5 infrate trade schooling invest60 invest70 invest80 industry
## 16 1 infrate trade schooling invest60 invest70 invest80 industry
## 16 2 infrate trade schooling invest60 invest70 invest80 industry
## 16 3 infrate trade schooling invest60 invest70 invest80 industry
## 16 4 infrate trade schooling invest60 invest70 invest80 industry
## 16 5 infrate trade schooling invest60 invest70 invest80 industry
## 17 1 infrate trade schooling invest60 invest70 invest80 industry
## 17 2 infrate trade schooling invest60 invest70 invest80 industry
## 17 3 infrate trade schooling invest60 invest70 invest80 industry
## 17 4 infrate trade schooling invest60 invest70 invest80 industry
## 17 5 infrate trade schooling invest60 invest70 invest80 industry
## 18 1 infrate trade schooling invest60 invest70 invest80 industry
## 18 2 infrate trade schooling invest60 invest70 invest80 industry
## 18 3 infrate trade schooling invest60 invest70 invest80 industry
## 18 4 infrate trade schooling invest60 invest70 invest80 industry
## 18 5 infrate trade schooling invest60 invest70 invest80 industry
## 19 1 infrate trade schooling invest60 invest70 invest80 industry
## 19 2 infrate trade schooling invest60 invest70 invest80 industry
## 19 3 infrate trade schooling invest60 invest70 invest80 industry
## 19 4 infrate trade schooling invest60 invest70 invest80 industry
## 19 5 infrate trade schooling invest60 invest70 invest80 industry
## 20 1 infrate trade schooling invest60 invest70 invest80 industry
## 20 2 infrate trade schooling invest60 invest70 invest80 industry
## 20 3 infrate trade schooling invest60 invest70 invest80 industry
## 20 4 infrate trade schooling invest60 invest70 invest80 industry
## 20 5 infrate trade schooling invest60 invest70 invest80 industry
## 21 1 infrate trade schooling invest60 invest70 invest80 industry
## 21 2 infrate trade schooling invest60 invest70 invest80 industry
## 21 3 infrate trade schooling invest60 invest70 invest80 industry
## 21 4 infrate trade schooling invest60 invest70 invest80 industry
## 21 5 infrate trade schooling invest60 invest70 invest80 industry
## 22 1 infrate trade schooling invest60 invest70 invest80 industry
## 22 2 infrate trade schooling invest60 invest70 invest80 industry
## 22 3 infrate trade schooling invest60 invest70 invest80 industry
## 22 4 infrate trade schooling invest60 invest70 invest80 industry
## 22 5 infrate trade schooling invest60 invest70 invest80 industry
## 23 1 infrate trade schooling invest60 invest70 invest80 industry
## 23 2 infrate trade schooling invest60 invest70 invest80 industry
## 23 3 infrate trade schooling invest60 invest70 invest80 industry
## 23 4 infrate trade schooling invest60 invest70 invest80 industry
## 23 5 infrate trade schooling invest60 invest70 invest80 industry
## 24 1 infrate trade schooling invest60 invest70 invest80 industry
## 24 2 infrate trade schooling invest60 invest70 invest80 industry
## 24 3 infrate trade schooling invest60 invest70 invest80 industry
## 24 4 infrate trade schooling invest60 invest70 invest80 industry
## 24 5 infrate trade schooling invest60 invest70 invest80 industry
## 25 1 infrate trade schooling invest60 invest70 invest80 industry
## 25 2 infrate trade schooling invest60 invest70 invest80 industry
## 25 3 infrate trade schooling invest60 invest70 invest80 industry
## 25 4 infrate trade schooling invest60 invest70 invest80 industry
## 25 5 infrate trade schooling invest60 invest70 invest80 industry
## 26 1 infrate trade schooling invest60 invest70 invest80 industry
## 26 2 infrate trade schooling invest60 invest70 invest80 industry
## 26 3 infrate trade schooling invest60 invest70 invest80 industry
## 26 4 infrate trade schooling invest60 invest70 invest80 industry
## 26 5 infrate trade schooling invest60 invest70 invest80 industry
## 27 1 infrate trade schooling invest60 invest70 invest80 industry
## 27 2 infrate trade schooling invest60 invest70 invest80 industry
## 27 3 infrate trade schooling invest60 invest70 invest80 industry
## 27 4 infrate trade schooling invest60 invest70 invest80 industry
## 27 5 infrate trade schooling invest60 invest70 invest80 industry
## 28 1 infrate trade schooling invest60 invest70 invest80 industry
## 28 2 infrate trade schooling invest60 invest70 invest80 industry
## 28 3 infrate trade schooling invest60 invest70 invest80 industry
## 28 4 infrate trade schooling invest60 invest70 invest80 industry
## 28 5 infrate trade schooling invest60 invest70 invest80 industry
## 29 1 infrate trade schooling invest60 invest70 invest80 industry
## 29 2 infrate trade schooling invest60 invest70 invest80 industry
## 29 3 infrate trade schooling invest60 invest70 invest80 industry
## 29 4 infrate trade schooling invest60 invest70 invest80 industry
## 29 5 infrate trade schooling invest60 invest70 invest80 industry
## 30 1 infrate trade schooling invest60 invest70 invest80 industry
## 30 2 infrate trade schooling invest60 invest70 invest80 industry
## 30 3 infrate trade schooling invest60 invest70 invest80 industry
## 30 4 infrate trade schooling invest60 invest70 invest80 industry
## 30 5 infrate trade schooling invest60 invest70 invest80 industry
## 31 1 infrate trade schooling invest60 invest70 invest80 industry
## 31 2 infrate trade schooling invest60 invest70 invest80 industry
## 31 3 infrate trade schooling invest60 invest70 invest80 industry
## 31 4 infrate trade schooling invest60 invest70 invest80 industry
## 31 5 infrate trade schooling invest60 invest70 invest80 industry
## 32 1 infrate trade schooling invest60 invest70 invest80 industry
## 32 2 infrate trade schooling invest60 invest70 invest80 industry
## 32 3 infrate trade schooling invest60 invest70 invest80 industry
## 32 4 infrate trade schooling invest60 invest70 invest80 industry
## 32 5 infrate trade schooling invest60 invest70 invest80 industry
## 33 1 infrate trade schooling invest60 invest70 invest80 industry
## 33 2 infrate trade schooling invest60 invest70 invest80 industry
## 33 3 infrate trade schooling invest60 invest70 invest80 industry
## 33 4 infrate trade schooling invest60 invest70 invest80 industry
## 33 5 infrate trade schooling invest60 invest70 invest80 industry
## 34 1 infrate trade schooling invest60 invest70 invest80 industry
## 34 2 infrate trade schooling invest60 invest70 invest80 industry
## 34 3 infrate trade schooling invest60 invest70 invest80 industry
## 34 4 infrate trade schooling invest60 invest70 invest80 industry
## 34 5 infrate trade schooling invest60 invest70 invest80 industry
## 35 1 infrate trade schooling invest60 invest70 invest80 industry
## 35 2 infrate trade schooling invest60 invest70 invest80 industry
## 35 3 infrate trade schooling invest60 invest70 invest80 industry
## 35 4 infrate trade schooling invest60 invest70 invest80 industry
## 35 5 infrate trade schooling invest60 invest70 invest80 industry
## 36 1 infrate trade schooling invest60 invest70 invest80 industry
## 36 2 infrate trade schooling invest60 invest70 invest80 industry
## 36 3 infrate trade schooling invest60 invest70 invest80 industry
## 36 4 infrate trade schooling invest60 invest70 invest80 industry
## 36 5 infrate trade schooling invest60 invest70 invest80 industry
## 37 1 infrate trade schooling invest60 invest70 invest80 industry
## 37 2 infrate trade schooling invest60 invest70 invest80 industry
## 37 3 infrate trade schooling invest60 invest70 invest80 industry
## 37 4 infrate trade schooling invest60 invest70 invest80 industry
## 37 5 infrate trade schooling invest60 invest70 invest80 industry
## 38 1 infrate trade schooling invest60 invest70 invest80 industry
## 38 2 infrate trade schooling invest60 invest70 invest80 industry
## 38 3 infrate trade schooling invest60 invest70 invest80 industry
## 38 4 infrate trade schooling invest60 invest70 invest80 industry
## 38 5 infrate trade schooling invest60 invest70 invest80 industry
## 39 1 infrate trade schooling invest60 invest70 invest80 industry
## 39 2 infrate trade schooling invest60 invest70 invest80 industry
## 39 3 infrate trade schooling invest60 invest70 invest80 industry
## 39 4 infrate trade schooling invest60 invest70 invest80 industry
## 39 5 infrate trade schooling invest60 invest70 invest80 industry
## 40 1 infrate trade schooling invest60 invest70 invest80 industry
## 40 2 infrate trade schooling invest60 invest70 invest80 industry
## 40 3 infrate trade schooling invest60 invest70 invest80 industry
## 40 4 infrate trade schooling invest60 invest70 invest80 industry
## 40 5 infrate trade schooling invest60 invest70 invest80 industry
## 41 1 infrate trade schooling invest60 invest70 invest80 industry
## 41 2 infrate trade schooling invest60 invest70 invest80 industry
## 41 3 infrate trade schooling invest60 invest70 invest80 industry
## 41 4 infrate trade schooling invest60 invest70 invest80 industry
## 41 5 infrate trade schooling invest60 invest70 invest80 industry
## 42 1 infrate trade schooling invest60 invest70 invest80 industry
## 42 2 infrate trade schooling invest60 invest70 invest80 industry
## 42 3 infrate trade schooling invest60 invest70 invest80 industry
## 42 4 infrate trade schooling invest60 invest70 invest80 industry
## 42 5 infrate trade schooling invest60 invest70 invest80 industry
## 43 1 infrate trade schooling invest60 invest70 invest80 industry
## 43 2 infrate trade schooling invest60 invest70 invest80 industry
## 43 3 infrate trade schooling invest60 invest70 invest80 industry
## 43 4 infrate trade schooling invest60 invest70 invest80 industry
## 43 5 infrate trade schooling invest60 invest70 invest80 industry
## 44 1 infrate trade schooling invest60 invest70 invest80 industry
## 44 2 infrate trade schooling invest60 invest70 invest80 industry
## 44 3 infrate trade schooling invest60 invest70 invest80 industry
## 44 4 infrate trade schooling invest60 invest70 invest80 industry
## 44 5 infrate trade schooling invest60 invest70 invest80 industry
## 45 1 infrate trade schooling invest60 invest70 invest80 industry
## 45 2 infrate trade schooling invest60 invest70 invest80 industry
## 45 3 infrate trade schooling invest60 invest70 invest80 industry
## 45 4 infrate trade schooling invest60 invest70 invest80 industry
## 45 5 infrate trade schooling invest60 invest70 invest80 industry
## 46 1 infrate trade schooling invest60 invest70 invest80 industry
## 46 2 infrate trade schooling invest60 invest70 invest80 industry
## 46 3 infrate trade schooling invest60 invest70 invest80 industry
## 46 4 infrate trade schooling invest60 invest70 invest80 industry
## 46 5 infrate trade schooling invest60 invest70 invest80 industry
## 47 1 infrate trade schooling invest60 invest70 invest80 industry
## 47 2 infrate trade schooling invest60 invest70 invest80 industry
## 47 3 infrate trade schooling invest60 invest70 invest80 industry
## 47 4 infrate trade schooling invest60 invest70 invest80 industry
## 47 5 infrate trade schooling invest60 invest70 invest80 industry
## 48 1 infrate trade schooling invest60 invest70 invest80 industry
## 48 2 infrate trade schooling invest60 invest70 invest80 industry
## 48 3 infrate trade schooling invest60 invest70 invest80 industry
## 48 4 infrate trade schooling invest60 invest70 invest80 industry
## 48 5 infrate trade schooling invest60 invest70 invest80 industry
## 49 1 infrate trade schooling invest60 invest70 invest80 industry
## 49 2 infrate trade schooling invest60 invest70 invest80 industry
## 49 3 infrate trade schooling invest60 invest70 invest80 industry
## 49 4 infrate trade schooling invest60 invest70 invest80 industry
## 49 5 infrate trade schooling invest60 invest70 invest80 industry
## 50 1 infrate trade schooling invest60 invest70 invest80 industry
## 50 2 infrate trade schooling invest60 invest70 invest80 industry
## 50 3 infrate trade schooling invest60 invest70 invest80 industry
## 50 4 infrate trade schooling invest60 invest70 invest80 industry
## 50 5 infrate trade schooling invest60 invest70 invest80 industry## Warning: Number of logged events: 1501d_imputed <- complete(mice_mod, 1)# Fit the bpCausal model
out <- bpCausal(data = d_imputed, ## dataset
index = c("country", "year"), ## names for unit and time index
Yname = "gdp", ## outcome variable
Dname = "D_i", ## treatment indicator
Xname = c(), # covariates that have constant (fixed) effect
Zname = c(), # covariates that have unit-level random effect
Aname = c("infrate", "trade", "schooling", "invest60", "invest70", "invest80", "industry"), # covariates that have time-level random effect
re = "both", # two-way random effect: choose from ("unit", "time", "none", "both")
ar1 = TRUE, # whether the time-level random effects is ar1 process or just multilevel (independent)
r = 8, # factor numbers
niter = 15000, # number of mcmc draws
burn = 5000, # burn-in draws
xlasso = 1, ## whether to shrink constant coefs (1 = TRUE, 0 = FALSE)
zlasso = 1, ## whether to shrink unit-level random coefs (1 = TRUE, 0 = FALSE)
alasso = 1, ## whether to shrink time-level coefs (1 = TRUE, 0 = FALSE)
flasso = 1, ## whether to shrink factor loadings (1 = TRUE, 0 = FALSE)
a1 = 0.001, a2 = 0.001, ## parameters for hyper prior shrink on beta (diffuse hyper priors)
b1 = 0.001, b2 = 0.001, ## parameters for hyper prior shrink on alpha_i
c1 = 0.001, c2 = 0.001, ## parameters for hyper prior shrink on xi_t
p1 = 0.001, p2 = 0.001) ## parameters for hyper prior shrink on factor terms## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000sout <- coefSummary(out) ## summary estimated parameters
eout <- effSummary(out, ## summary treatment effects
usr.id = NULL, ## treatment effect for individual treated units, if input NULL, calculate average TT
cumu = FALSE, ## whether to calculate culmulative treatment effects
rela.period = TRUE) ## whether to use time relative to the occurence of treatment (1 is the first post-treatment period) or real period (like year 1998, 1999, ...)
est_att <- eout$est.eff
est_att$time <- c(1960:2003)4 Tabular Aggregation
SCM_ATT_3 <- full_join(est_att, SCM, by = c("time"))
SCM_ATT_3## observed estimated_counterfactual counterfactual_ci_l counterfactual_ci_u
## 1 2284 2241.130 2575.420 1910.236
## 2 2388 2376.520 2703.343 2045.095
## 3 2527 2469.111 2794.070 2134.376
## 4 2610 2584.522 2915.667 2261.693
## 5 2806 2780.498 3107.255 2459.746
## 6 3005 2952.279 3279.529 2623.875
## 7 3168 3118.801 3450.095 2776.276
## 8 3241 3334.092 3663.743 3006.224
## 9 3571 3651.171 3989.872 3307.506
## 10 3998 4046.340 4385.983 3702.604
## 11 4367 4408.281 4733.618 4081.521
## 12 4686 4649.851 5023.309 4277.951
## 13 5055 5010.213 5383.489 4643.856
## 14 5553 5506.717 5869.146 5146.532
## 15 6074 6174.670 6542.595 5805.801
## 16 6603 6707.385 7080.147 6330.092
## 17 7367 7445.730 7805.102 7091.286
## 18 8090 8117.009 8460.303 7761.763
## 19 8928 8902.026 9248.897 8548.204
## 20 10067 10012.212 10370.379 9656.641
## 21 11083 11032.916 11384.794 10690.030
## 22 12115 12091.428 12425.664 11766.698
## 23 12761 12736.086 13079.842 12397.788
## 24 13519 13507.098 13834.219 13171.796
## 25 14481 14438.717 14779.151 14100.484
## 26 15291 15298.593 15640.451 14965.635
## 27 15998 16005.296 16346.896 15660.036
## 28 16679 16704.963 17046.886 16361.634
## 29 17786 17765.700 18112.748 17423.219
## 30 18994 18962.508 19316.589 18605.019
## 31 20465 20526.404 20896.748 20145.736
## 32 21602 21412.769 22043.987 20777.702
## 33 22154 21946.413 22627.548 21230.517
## 34 21878 21693.048 22390.881 20989.574
## 35 22371 23585.675 24313.644 22846.708
## 36 23035 24370.746 25155.826 23583.714
## 37 23742 26284.047 27401.737 25112.294
## 38 24156 27503.694 28601.199 26325.541
## 39 24931 27143.020 28192.753 26075.984
## 40 25755 28593.844 29916.417 27270.554
## 41 26943 29765.502 31406.962 28180.548
## 42 27449 30625.491 32662.531 28633.451
## 43 28348 31491.083 33413.690 29637.719
## 44 28855 32964.288 34864.982 31080.809
## estimated_ATT estimated_ATT_ci_l estimated_ATT_ci_u time count real_gdp_WG
## 1 42.869723 -291.4201 373.7640 1960 1 2284
## 2 11.480457 -315.3431 342.9047 1961 1 2388
## 3 57.888864 -267.0699 392.6239 1962 1 2527
## 4 25.478223 -305.6671 348.3074 1963 1 2610
## 5 25.502110 -301.2553 346.2539 1964 1 2806
## 6 52.720930 -274.5289 381.1253 1965 1 3005
## 7 49.198598 -282.0951 391.7241 1966 1 3168
## 8 -93.091596 -422.7427 234.7764 1967 1 3241
## 9 -80.170822 -418.8716 263.4944 1968 1 3571
## 10 -48.340140 -387.9835 295.3962 1969 1 3998
## 11 -41.280712 -366.6184 285.4791 1970 1 4367
## 12 36.149197 -337.3087 408.0493 1971 1 4686
## 13 44.786808 -328.4889 411.1439 1972 1 5055
## 14 46.283217 -316.1463 406.4679 1973 1 5553
## 15 -100.670189 -468.5948 268.1986 1974 1 6074
## 16 -104.385074 -477.1473 272.9079 1975 1 6603
## 17 -78.729832 -438.1024 275.7137 1976 1 7367
## 18 -27.008619 -370.3035 328.2366 1977 1 8090
## 19 25.974071 -320.8970 379.7959 1978 1 8928
## 20 54.788045 -303.3786 410.3595 1979 1 10067
## 21 50.084000 -301.7938 392.9703 1980 1 11083
## 22 23.571677 -310.6642 348.3022 1981 1 12115
## 23 24.913824 -318.8423 363.2121 1982 1 12761
## 24 11.901940 -315.2194 347.2040 1983 1 13519
## 25 42.283027 -298.1509 380.5164 1984 1 14481
## 26 -7.592729 -349.4513 325.3646 1985 1 15291
## 27 -7.295983 -348.8963 337.9642 1986 1 15998
## 28 -25.962928 -367.8860 317.3657 1987 1 16679
## 29 20.299969 -326.7484 362.7808 1988 1 17786
## 30 31.491853 -322.5887 388.9815 1989 1 18994
## 31 -61.404173 -431.7485 319.2640 1990 1 20465
## 32 189.231154 -441.9873 824.2984 1991 1 21602
## 33 207.587249 -473.5477 923.4835 1992 1 22154
## 34 184.952235 -512.8810 888.4262 1993 1 21878
## 35 -1214.674896 -1942.6440 -475.7080 1994 1 22371
## 36 -1335.746405 -2120.8256 -548.7136 1995 1 23035
## 37 -2542.046551 -3659.7373 -1370.2936 1996 1 23742
## 38 -3347.693677 -4445.1990 -2169.5411 1997 1 24156
## 39 -2212.020274 -3261.7528 -1144.9839 1998 1 24931
## 40 -2838.843998 -4161.4165 -1515.5536 1999 1 25755
## 41 -2822.502359 -4463.9618 -1237.5478 2000 1 26943
## 42 -3176.490997 -5213.5305 -1184.4509 2001 1 27449
## 43 -3143.082892 -5065.6898 -1289.7188 2002 1 28348
## 44 -4109.287574 -6009.9816 -2225.8093 2003 1 28855
## synth_gdp_WG_GSCM synth_gdp_WG_BSCM synth_gdp_WG_cred_int_2.5
## 1 2116.440 2067.584 1935.338
## 2 2218.093 2200.744 2076.437
## 3 2340.969 2335.087 2195.265
## 4 2451.913 2463.256 2321.915
## 5 2626.648 2640.120 2487.879
## 6 2782.204 2807.034 2648.566
## 7 3008.168 3020.359 2855.984
## 8 3183.610 3203.061 3028.299
## 9 3483.878 3519.100 3348.152
## 10 3863.808 3925.046 3755.645
## 11 4289.694 4323.007 4150.061
## 12 4663.123 4665.820 4495.901
## 13 5083.388 5068.631 4884.165
## 14 5602.182 5612.788 5423.441
## 15 6181.158 6180.164 5978.878
## 16 6672.119 6683.389 6495.726
## 17 7303.438 7342.677 7159.881
## 18 8058.963 8008.943 7834.500
## 19 8785.813 8795.937 8612.853
## 20 9869.466 9841.252 9652.084
## 21 10907.325 10911.498 10722.171
## 22 12001.454 11990.890 11803.546
## 23 12715.357 12753.963 12552.612
## 24 13541.532 13479.628 13298.234
## 25 14393.761 14364.420 14177.112
## 26 15275.415 15235.698 15050.070
## 27 15952.404 15943.414 15760.038
## 28 16687.287 16728.066 16536.065
## 29 17853.879 17928.190 17732.491
## 30 19159.964 19190.810 18997.033
## 31 20460.668 20458.084 20250.582
## 32 21338.468 21457.246 21150.762
## 33 22070.885 22127.132 21791.670
## 34 22525.785 22512.562 22189.743
## 35 23440.665 23455.112 23112.115
## 36 24258.348 24399.413 23982.564
## 37 25223.701 25377.303 24829.251
## 38 26126.344 26323.808 25717.377
## 39 26961.147 27071.445 26524.995
## 40 27944.793 27991.598 27282.949
## 41 29598.453 29863.002 28924.840
## 42 30277.618 30844.303 29856.898
## 43 31310.616 31873.340 30900.362
## 44 32237.649 32591.531 31592.277
## synth_gdp_WG_cred_int_97.5 diff_GSCM diff_BSCM diff_BSCM_lwr diff_BSCM_upr
## 1 2197.587 167.560 216.416 348.662 86.413
## 2 2330.747 169.907 187.256 311.563 57.253
## 3 2469.141 186.031 191.913 331.735 57.859
## 4 2607.192 158.087 146.744 288.085 2.808
## 5 2787.297 179.352 165.880 318.121 18.703
## 6 2965.729 222.796 197.966 356.434 39.271
## 7 3185.213 159.832 147.641 312.016 -17.213
## 8 3366.026 57.390 37.939 212.701 -125.026
## 9 3685.524 87.122 51.900 222.848 -114.524
## 10 4090.105 134.192 72.954 242.355 -92.105
## 11 4496.915 77.306 43.993 216.939 -129.915
## 12 4833.099 22.877 20.180 190.099 -147.099
## 13 5242.322 -28.388 -13.631 170.835 -187.322
## 14 5798.232 -49.182 -59.788 129.559 -245.232
## 15 6372.722 -107.158 -106.164 95.122 -298.722
## 16 6870.503 -69.119 -80.389 107.274 -267.503
## 17 7523.740 63.562 24.323 207.119 -156.740
## 18 8189.684 31.037 81.057 255.500 -99.684
## 19 8984.436 142.187 132.063 315.147 -56.436
## 20 10033.432 197.534 225.748 414.916 33.568
## 21 11107.706 175.675 171.502 360.829 -24.706
## 22 12178.423 113.546 124.110 311.454 -63.423
## 23 12942.676 45.643 7.037 208.388 -181.676
## 24 13660.958 -22.532 39.372 220.766 -141.958
## 25 14559.386 87.239 116.580 303.888 -78.386
## 26 15426.274 15.585 55.302 240.930 -135.274
## 27 16126.356 45.596 54.586 237.962 -128.356
## 28 16915.060 -8.287 -49.066 142.935 -236.060
## 29 18116.902 -67.879 -142.190 53.509 -330.902
## 30 19377.836 -165.964 -196.810 -3.033 -383.836
## 31 20664.513 4.332 6.916 214.418 -199.513
## 32 21756.504 263.532 144.754 451.238 -154.504
## 33 22462.489 83.115 26.868 362.330 -308.489
## 34 22832.299 -647.785 -634.562 -311.743 -954.299
## 35 23792.998 -1069.665 -1084.112 -741.115 -1421.998
## 36 24819.090 -1223.348 -1364.413 -947.564 -1784.090
## 37 25921.594 -1481.701 -1635.303 -1087.251 -2179.594
## 38 26898.328 -1970.344 -2167.808 -1561.377 -2742.328
## 39 27601.415 -2030.147 -2140.445 -1593.995 -2670.415
## 40 28677.952 -2189.793 -2236.598 -1527.949 -2922.952
## 41 30767.938 -2655.453 -2920.002 -1981.840 -3824.938
## 42 31770.437 -2828.618 -3395.303 -2407.898 -4321.437
## 43 32813.055 -2962.616 -3525.340 -2552.362 -4465.055
## 44 33581.982 -3382.649 -3736.531 -2737.277 -4726.9825 Trends in Per Capita GDP: West Germany under 3 Synthetic Control Methods
ggplot(SCM_ATT_3) +
geom_ribbon(aes(x = time, y = estimated_counterfactual,
ymin = counterfactual_ci_l,
ymax = counterfactual_ci_u),
alpha = 0.8, fill = "pink") +
geom_ribbon(aes(x = time, y = synth_gdp_WG_BSCM,
ymin = synth_gdp_WG_cred_int_2.5,
ymax = synth_gdp_WG_cred_int_97.5),
alpha = 0.9, fill = "steelblue1") +
geom_line(aes(x = time, y = observed, color = "darkorange4"), size = 0.7) +
geom_line(aes(x = time, y = synth_gdp_WG_GSCM, color = "yellow2"),
linetype = "solid", size = 0.7) +
geom_line(aes(x = time, y = estimated_counterfactual, color = "red2"),
linetype = "longdash", size = 0.7) +
geom_line(aes(x = time, y = synth_gdp_WG_BSCM, color = "darkblue"),
linetype = "dashed", size = 0.7) +
geom_vline(xintercept = 1990, linetype = "dotted", color = "darkgrey", size = 2) +
xlab("year") + ylab("Per Capita GDP (PPP, 2002 USD)") +
scale_color_identity(name = "Per Capita GDP",
breaks = c("darkorange4", "yellow2", "darkblue", "red2"),
labels = c("Observed West Germany",
"SCM (Standard) - West Germany",
"BSCM (IFX) - West Germany",
"BSCM (DM-LFM) - West Germany"),
guide = "legend") +
theme(legend.position = c(0.8, 0.2)) +
annotate("text", x = 1984, y = 31000, label = "German Reunification") +
annotate("text", x = 1992.5, y = 33000,
label = "BSCM (DM-LFM) - 95% CI", color = "pink") +
annotate("text", x = 1991.5, y = 28500,
label = "BSCM (IFX) - 95% CI", color = "steelblue1") +
ggtitle("Trends in Per Capita GDP: West Germany under 3 Synthetic Control Methods")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
6 Esimated Treatment Effects in Per Capita GDP for West Germany under 3 SCMs
ggplot(SCM_ATT_3) +
geom_ribbon(aes(x = time, y = estimated_ATT,
ymin = estimated_ATT_ci_l,
ymax = estimated_ATT_ci_u),
alpha = 0.8, fill = "pink") +
geom_ribbon(aes(x = time, y = diff_BSCM,
ymin = diff_BSCM_lwr,
ymax = diff_BSCM_upr),
alpha = 0.9, fill = "steelblue1") +
geom_line(aes(x = time, y = estimated_ATT, color = "red2"),
linetype = "longdash", size = 1.1) +
geom_line(aes(x = time, y = diff_BSCM, color = "darkblue"),
linetype = "dashed", size = 1.1) +
geom_line(aes(x = time, y = diff_GSCM, color = "yellow2"),
linetype = "solid", size = 1.1) +
geom_vline(xintercept = 1990, linetype = "dashed", size = 0.5, col = "darkgrey") +
geom_hline(yintercept = 0, linetype = "dashed", size = 0.5, col = "darkgrey") +
scale_color_identity(name = "Per Capita GDP",
breaks = c("yellow2", "darkblue", "red2"),
labels = c("SCM (Standard) - West Germany",
"BSCM (IFX) - West Germany",
"BSCM (DM-LFM) - West Germany"),
guide = "legend") +
theme(legend.position = c(0.3, 0.2)) +
annotate("text", x = 1984, y = -700, label = "German Reunification") +
annotate("text", x = 1987, y = -2000,
label = "BSCM (DM-LFM) - 95% CI", color = "pink") +
annotate("text", x = 2000, y = -750,
label = "BSCM (IFX) - 95% CI", color = "steelblue1") +
ylab("Gap in Per Capita GDP (PPP, 2002 USD)") + xlab("year") +
ggtitle("Esimated Treatment Effects in Per Capita GDP for West Germany under 3 SCMs") 
7 ATT Estimator
TT_df <- filter(SCM_ATT_3[, c("time", "estimated_ATT", "estimated_ATT_ci_l", "estimated_ATT_ci_u", "diff_BSCM", "diff_BSCM_lwr", "diff_BSCM_upr", "diff_GSCM")], time > 1990)[, -1]
ATT_list <- lapply(TT_df, mean)
ATT_df <- as.data.frame(ATT_list)
colnames(ATT_df) <- c("ATT_BSCM(DM-LFM)", "ATT_BSCM(DM-LFM)_lwr", "ATT_BSCM(DM-LFM)_upr", "ATT_BSCM(IFX)", "ATT_BSCM(IFX)_upr", "ATT_BSCM(IFX)_lwr", "ATT_SCM")
ATT_df <- ATT_df[, c("ATT_BSCM(DM-LFM)_lwr", "ATT_BSCM(DM-LFM)", "ATT_BSCM(DM-LFM)_upr", "ATT_BSCM(IFX)_lwr", "ATT_BSCM(IFX)", "ATT_BSCM(IFX)_upr", "ATT_SCM")]
t(ATT_df)## [,1]
## ATT_BSCM(DM-LFM)_lwr -3213.320
## ATT_BSCM(DM-LFM) -2012.355
## ATT_BSCM(DM-LFM)_upr -809.701
## ATT_BSCM(IFX)_lwr -2498.237
## ATT_BSCM(IFX) -1897.600
## ATT_BSCM(IFX)_upr -1279.754
## ATT_SCM -1699.6528 References
[1] S. Pinkney, “An Improved and Extended Bayesian Synthetic Control,” arXiv.org, Mar. 30, 2021. https://arxiv.org/abs/2103.16244.
[2] X. Pang, L. Liu, and Y. Xu, “A Bayesian Alternative to Synthetic Control for Comparative Case Studies,” SSRN Electronic Journal, 2020, doi: https://doi.org/10.2139/ssrn.3649226.
9 Appendix - STAN codes
functions {
matrix make_F (int T,
vector diagonal_loadings,
vector lower_tri_loadings) {
int L = num_elements(diagonal_loadings);
int M = num_elements(lower_tri_loadings);
matrix[T, L] F;
int idx = 0; // Index for the lower diagonal
for (j in 1:L) {
F[j, j] = diagonal_loadings[j];
for(i in (j + 1):T) {
idx += 1;
F[i, j] = lower_tri_loadings[idx];
}
}
for (j in 1:(L - 1)) {
for (i in (j + 1):L) F[j, i] = 0;
}
return F;
}
matrix make_beta (int J, matrix off,
vector lambda,
real eta,
vector tau) {
int L = cols(off);
vector[L] cache = ( tan(0.5 * pi() * lambda) *
tan(0.5 * pi() * eta) );
vector[J] tau_ = tan(0.5 * pi() * tau);
matrix[J, L] out;
for (j in 1:J)
out[j] = off[j] * tau_[j];
return diag_pre_multiply(cache, out');
}
}
data {
int T; // times
int J; // countries
int L; // number of factors
int P;
matrix[P, J] X; // predictors
row_vector[T] Y[J]; // data matrix of order [J,T]
int trt_times;
}
transformed data {
int<lower=1> M = L * (T - L) + L * (L - 1) / 2;
row_vector[J] j_ones = rep_row_vector(1, J);
vector[T] t_ones = rep_vector(1.0, T);
matrix[J, P] X_std;
vector[J] y_mu;
vector[J] y_sd;
row_vector[T] Y_scaled[J];
row_vector[T - trt_times] Y_pre_target;
vector[P] x_mu;
vector[P] x_sd;
for (j in 1:J) {
y_mu[j] = Y[j, T - trt_times];
y_sd[j] = sd(Y[j, 1:T - trt_times]);
Y_scaled[j] = ( Y[j] - y_mu[j] ) / y_sd[j];
}
for (p in 1:P) {
x_mu[p] = mean(X[p]);
x_sd[p] = sd(X[p]);
X_std[, p] = ( X[p]' - mean(X[p]) ) / sd(X[p]);
}
Y_pre_target = Y_scaled[1, 1:T - trt_times];
}
parameters{
vector[P] chi; // non-time varying predictors
vector[T] delta; // year fixed effects
row_vector[J] kappa; // country fixed effects
matrix[J, L] beta_off;
vector<lower=0, upper=1>[L] lambda;
real<lower=0, upper=1> eta;
vector<lower=0, upper=1>[J] tau;
row_vector[trt_times] Y_post_target;
real<lower=0> sigma;
vector<lower=0>[L] F_diag;
vector[M] F_lower;
}
transformed parameters {
matrix[L, J] beta = make_beta(J,
beta_off,
lambda,
eta,
tau);
}
model {
chi ~ std_normal();
to_vector(beta_off) ~ std_normal();
F_diag ~ std_normal();
F_lower ~ normal(0, 2);
delta ~ normal(0, 2);
kappa ~ std_normal();
sigma ~ std_normal();
{
vector[J] predictors = X_std * chi;
matrix[T, L] F = make_F(T, F_diag, F_lower);
row_vector[T] Y_target[1];
row_vector[T] Y_temp[J];
Y_target[1] = append_col(Y_pre_target, Y_post_target);
Y_temp = append_array(Y_target, Y_scaled[2:J]);
for (j in 1:J)
Y_temp[j]' ~ normal_id_glm(F, delta + kappa[j] + predictors[j],
beta[ , j], sigma);
}
}
generated quantities {
vector[T] synth_out[J];
matrix[T, L] F_ = make_F(T, F_diag, F_lower);
matrix[T, J] Synth_ = F_ * beta + delta * j_ones +
t_ones * (kappa + (X_std * chi)');
for (j in 1:J)
synth_out[j] = Synth_[, j] * y_sd[j] + y_mu[j];
}